Drainage channels are essential components of englacial and subglacial hydrologic systems. Here, we use the M integral, a path-independent integral of the equations of continuum mechanics for a class of media, to unify descriptions of creep closure under a variety of stress states surrounding drainage channels. The advantage of this approach is that the M integral around the hydrologic channels is identical to same integral evaluated in the far field. In this way, the creep closure on the channel wall can be determined as a function of the far-field loading, e.g., involving antiplane shear as well as overburden pressure. We start by analyzing the axisymmetric case and show that the Nye solution for the creep closure of the channels is implied by the path independence of the M integral. We then examine the effects of superimposing antiplane shear. We show that the creep closure of the channels acts as a perturbation in the far field, which we explore analytically and numerically. In this way, the creep closure of channels can be succinctly written in terms of the path-independent M integral, and understanding the variation with applied shear is useful for glacial hydrology models.
Glacial melt water from surface ablation, precipitation, and internal deformation drains via englacial conduits to the base of the glacier where it is evacuated through a subglacial hydrologic network of channels melted into the ice or cavities in the sediments [1–3]. In this paper, we focus on channelized drainage through Röthlisberger channels , where these channels are melted into the ice by the heat dissipation of the turbulently flowing melt water and close by viscous creep of the surrounding ice. We model these channels as very long straight conduits that are oriented along the direction of glacier flow, and we use a conserved integral approach to derive the classical solution found by Nye  for the radial, or in-plane, creep closure velocity of the ice into the channel. We then show how the creep closure of the channels increases when we take the downstream shear present within the ice column, referred to as antiplane shear, into account, such as in an ice stream shear margin.
as written by Knowles and Sternberg .
where ni is the unit outer normal to the closed contour Γ, and s is arc length anticlockwise around the path, such that n1ds = dx2 and n2ds = −dx1, where the tensor subscripts correspond as (x, y, z) = (x1, x2, x3) and (ux, uy, uz) = (u1, u2, u3). Figure 1 shows the domain, the coordinate system, and a path of integration about a void.
where the first condition shows that Dij is symmetric, and by the second condition, the mass conservation for an incompressible substance Dkk = 0.
and, thus, m = 1 + 1/n. From here on, we will use the rheological exponent n instead of m; for a Newtonian viscous fluid (or in linear elasticity), n = 1 and, therefore, m = 2, so the term proportional to vi in Eq. (2) will disappear.
Therefore, our analysis primarily applies for ice modeled using a power-law rheology for ice, such as Glen's law.
In this expression, there are two types of terms: in-plane and antiplane. The in-plane terms are those in the r and θ directions, such as vr, vθ, and σrθ. The antiplane terms, quantities with a subscript x, represent motion into and out of page as a function of only the in-plane coordinates (using the standard glaciology coordinate system with z vertical, y across glacier, and x down glacier). We consider a very long channel with constant ice thickness and, in this way, we can reduce a three-dimensional problem to two dimensions where the quantities are homogeneous along x.
The setup for the problem and these conditions can be seen in Fig. 2. What is also evident is that the problem is purely in-plane and, therefore, we disregard the antiplane terms in Eq. (5). Furthermore, the problem is axisymmetric and so the in-plane shear terms are identically zero. Hence, we have the integral
This can be inserted into Eq. (7). Then, using the fact that M is path independent, we can evaluate two contours: first, along r = a and second, around r → ∞. These two contours are chosen because these are the locations where the boundary conditions are applied. Starting with the latter, we can see from the mass conservation that as r → ∞. This means that the flow potential W also decays to zero in the far field. Although vrr is a constant as r → ∞, the stress (boundary condition) and, therefore, M = 0 as r → ∞.
which corroborates  and, also, reduces to the standard result as b → ∞.
A natural extension for the M integral around an englacial or subglacial channel is to include antiplane terms. These are the terms in Eq. (5) that include x dependence. In glaciology, the antiplane terms can represent the shear flow of ice downstream, which is often ignored in the creep closure of channels [23–25]. However, the downstream flow decreases the effective viscosity of the ice, due to the fact that Glen's law is a shear-thinning rheology, and channels close more quickly than in environments free of antiplane stress. Nye  and Glen  compare the Nye solution to tunnel closure measurements in the field and find that some tunnels close much faster than predicted. Thus, the coupling between the in-plane creep closure and the antiplane motion of the glacier may be important in modeling subglacial hydrologic systems.
Ice stream shear margins are also examples of where antiplane effects can affect the size of drainage channels. Perol et al.  give theoretical arguments for the existence of subglacial channels beneath ice stream shear margins, which is backed up by observations of running water at the base of the dormant Kamb ice stream . Figure 3 shows a schematic for an idealized ice stream shear margin, where the velocity transitions from the fast-flowing ice stream to the nearly stagnant ridge. In the margin, we have schematically drawn a Röthlisberger  subglacial channel. Meyer et al.  show that the shear in the margin leads to faster closure velocities of drainage channels than would be predicted by the Nye solution, due to a decrease in effective viscosity from adding the antiplane shear.
which define the antiplane field. Using these boundary conditions, we evaluate the M integral around two contours: the channel wall r = a as well as in the far-field r → ∞.
Evaluation of the M Integral at the Channel.
Evaluation of M in the Far Field.
using the half-angle formula for . In order for M to be finite, this integral must be zero because of the r2 term in the integrand. However, this integral does not represent the full far-field contributions to M. The boundary conditions represent a constant strain rate, a situation where M = 0, and this evaluation of M in the far field disregards the coupling between the antiplane and in-plane motion. Therefore, we must retain a small perturbation away from a constant background strain rate.
This integral should be negative because in , the antiplane velocity should be less than the boundary value as it approaches the edge.
Thus, the perturbation amplitude ε is related to the unknown strain rates at the edge of the channel. Furthermore, we can see that if vx = 0 (or if vx is a function of r only), then this integral reduces to the integral for the Nye solution and ε = 0.
Nondimensional Equations and Numerics.
and represents the importance of antiplane shearing to in-plane creep closure . Using S we can rewrite the antiplane velocity as .
Although we cannot solve for the constant κ analytically, we can determine its value by numerical simulations. We implement the numerical method described in Ref.  in the existing finite element software comsol . These simulations allow us to compute the value of M as a function of the strain rate ratio S, which is shown in Fig. 4. The two regimes, where M scales as M ∼ S1∕n for S ≪ 1 and M ∼ S1∕n for S ≫ 1, are clearly visible. The best-fit value of κ determined from the simulations is given by
These results show that as the amount of antiplane shear with respect to in-plane shear is increased, as measured by an increase in the strain rate ratio S, the M integral also increases. This is due to a simultaneous increase in both the creep closure velocity Vr as well as the antiplane straining along the edge of the channel. The increase in channel closure velocity is due to a decrease in effective viscosity, as described in Ref. .
where more details are provided in Ref. . In Fig. 5, we show the two limits of the creep closure velocity: for S ≪ 1, the simulations approach the Nye solution, and for S ≫ 1, we verify the scaling given in Eq. (16).
This increase in creep closure velocity due to the addition of antiplane shear is consistent with the increase in tunnel closure velocities observed by Nye  and Glen  due to changes in the glacier stress state. Furthermore, the increase in creep closure velocity due to antiplane straining is analogous to the Rice and Tracey  effect where the growth of voids is strongly enhanced by triaxiality. Intuitively, the in-plane creep closure velocity for large S grows less than linearly with S as the antiplane field only influences the in-plane motion through the viscosity. The consequence is that the dominant balance for large S in Eq. (15) is still between the perturbation in the far field and the antiplane shear at the channel wall.
In this paper, we apply the path-independent M integral to the creep of ice around subglacial and englacial channels. We correct a longstanding error in the implementation of the M integral to problems in generalized elasticity and non-Newtonian power-law fluids. We then use this integral to derive the Nye solution for the creep closure of an englacial or subglacial channel. Building on this solution, we consider applications where the flow of ice is not entirely in-plane and axisymmetric but includes components of flow down glacier (i.e., parallel to the channel axis). Using a simple far-field shear as a representation for ice stream shear margins, we find that an in-plane perturbation exists in the far field. We solve for the perturbation explicitly, up to a constant factor ε. Then, using the M integral, we are able to relate this perturbation back to the in-plane strain rates at the channel. We show that the factor ε exhibits two scaling regimes based on the size of the strain rate ratio : For small antiplane velocity relative to in-plane creep closure, ε is a constant and the M integral approaches zero as M ∼ S1∕n for vanishingly small S. In the other limit, where antiplane straining dominates, M grows as . These two scaling regimes are also present in creep closure velocity where for small S, we retrieve the Nye solution, and for S ≫ 1, the closure velocity increases as due to a decrease in ice viscosity. Thus, M provides a succinct description of the processes affecting channel closure with superimposed antiplane shear.
We would like to thank M. C. Fernandes, T. Perol, and Z. Suo for insightful conversations. We acknowledge the support of NSF Graduate Research Fellowship (Grant No. DGE1144152) to C.R.M. and NSF Office of Polar Programs Grant No. PP1341499.
- a =
- A =
- DE =
effective strain rate,
- Dij =
strain rate field
- g =
acceleration due to gravity
- H =
- Ii =
strain rate tensor invariants
invariant of strain rate tensor
- ℓ =
dimension, i.e., two or three
- m =
homogeneity degree, Ws = σijεij/m
- M =
- n =
rheological exponent, i.e., 3 for Glen's law
- p =
isotropic ice pressure
- pf =
fluid pressure within channel
- S =
strain rate ratio,
- ui =
- vi =
- W =
- Ws =
strain energy density
- xi =
far field antiplane strain rate
- Δp =
- εij =
- ρi =
- σij =
- σ0 =
ice overburden pressure, ρigH
deviatoric stress field, σij + pδij
- τE =
Appendix A: Proof of the Path Independence of M
is path independent.
Appendix B: Flow Potential for a Reiner–Rivlin Fluid
where σij is the stress tensor, Dij is the symmetric part of the velocity gradient tensor, p is the isotropic pressure, and Ik is the kth invariant of the tensor Dij . By writing the stress as an isotropic matrix function of Dij, assuming a symmetric dependence on Dij and Dji, and expanding this function of as a power series, we can use the Cayley–Hamilton theorem to show that the stress is a quadratic polynomial in Dij with coefficients that are functions of the invariants of Dij. For an incompressible fluid with isotropic pressure, this reduces to Eq. (B1).
where the isotropic pressure p is independent of the strain rate. This shows that the deviatoric stress also satisfies the Maxwell reciprocity.
just as was found in Eq. (B3).
A class of Reiner–Rivlin fluids that always satisfy the condition in Eq. (B3) are those for which is a function of I2 solely and . This is the rheological structure of Glen's law in glaciology, and therefore, a flow potential W will always exist.