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.

## Introduction

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 [2], 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 [4] 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.

*M*integral. Using the Noether [8] procedure, the

*ℓ*-dimensional, linear elastic

*M*integral is the conservation integral that results from a self-similar scale change by the infinitesimal factor

*γ*, i.e.,

*x*and displacements

_{α}*u*are self-similarly scaled from the reference configuration [9]. In the framework of linearized kinematics, the strain is equal to the symmetric part of the displacement

_{α}*u*gradient tensor as $\epsilon ij=sym(\u2202ui/\u2202xj)$. We can then write a strain energy density

_{i}*W*as a product of stresses

_{s}*σ*and strains

_{ij}*ε*by

_{ij}*W*=

_{s}*σ*/2. The general conserved integral resulting from the Noether procedure, with $y\alpha =x\alpha \u2032\u2212x\alpha $ and $f\alpha (x\xaf)=u\alpha \u2032(x\xaf\u2032)\u2212u\alpha (x\xaf)$, is given by

_{ij}ε_{ij}*M*integral for a linear elastic material in two dimensions with linearized kinematics is given as

as written by Knowles and Sternberg [6].

*M*to a generalized elastic material with a strain energy density

*W*that is homogeneous of degree

_{s}*m*in the strains

*ε*, and therefore,

_{ij}*W*=

_{s}*σ*/

_{ij}ε_{ij}*m*. Unfortunately, the expression of Budiansky and Rice [7] for the generalized

*M*integral contains an error. Whether it is typographical, conceptual, or due to the printing process is unknown. He and Hutchinson [10] give a correct expression (although without derivation) for the three-dimensional generalized

*M*integral in a different geometry than used here or in Ref. [7]: a closed 3D void or flat crack of axially symmetric shape, such that the stresses vary with

*z*and $x2+y2$. Rice [9] gives the correct Noether invariant transformation to generate the

*ℓ*-dimensional

*M*integral for a power-law solid, although, subsequently, only writes the expression of the

*M*integral for the linear (

*m*= 2) material in two dimensions. To set the record straight, the generalized

*M*integral in two dimensions with combined in-plane and antiplane straining in the

*y–z*plane, and with void aligned in the

*x*direction, is

where *n _{i}* is the unit outer normal to the closed contour Γ, and

*s*is arc length anticlockwise around the path, such that

*n*

_{1}

*ds*=

*dx*

_{2}and

*n*

_{2}

*ds*= −

*dx*

_{1}, where the tensor subscripts correspond as (

*x*,

*y*,

*z*) = (

*x*

_{1},

*x*

_{2},

*x*

_{3}) and (

*u*,

_{x}*u*,

_{y}*u*) = (

_{z}*u*

_{1},

*u*

_{2},

*u*

_{3}). Figure 1 shows the domain, the coordinate system, and a path of integration about a void.

*M*for a power-law nonlinear elastic solid is exactly equivalent to the definition for a power-law nonlinear viscous fluid, where the displacement

*u*is replaced by the velocity

_{i}*v*and the strain

_{i}*ε*by the strain rate

_{ij}*D*[11]. The strain rate is the symmetric part of the velocity gradient tensor, as $Dij=sym(\u2202vi/\u2202xj)$. Under the definition for a power-law viscous fluid, the elastic strain energy density

_{ij}*W*is replaced by a function called the flow potential

_{s}*W*. Just as $dWs=\sigma ijdDij$ is an exact differential in generalized elasticity,

*dW*is also an exact differential, satisfying $dW=\sigma ijdDij$. From here on, we will use the notation related to the flow of a viscous fluid. This change of notation and extension to viscous fluids allow us to apply the

*M*integral to the flow of ice in glaciers. In the notation of viscous fluids, the two-dimensional

*M*integral is written as

*W*for an incompressible viscous fluid is given as

*ε*= 0), and

_{kk}*dW*is an exact differential. The strain rates are defined in terms of derivatives of the velocity

*v*as

_{i}where the first condition shows that *D _{ij}* is symmetric, and by the second condition, the mass conservation for an incompressible substance

*D*= 0.

_{kk}### Ice Rheology.

*m*. In glaciology, ice is commonly modeled as an incompressible viscous fluid with a nonlinear rheology called Glen's law

*D*is a function of the second invariant of the strain rate tensor, i.e., $DE=DijDij/2$,

_{E}*τ*is defined in the same way for the deviatoric stress tensor $\sigma \u2032ij=\sigma ij+p\delta ij$, and the subscript

_{E}*E*stands for effective. The two parameters are

*A*, the ice softness, and

*n*, the rheological exponent. The standard value used in glaciology is

*n*= 3 and is called Glen's law, which is appropriate for the typical values of stress and strain rate encountered in the field [12]. Now, relating the deviatoric stress and strain rate tensors using Eq. (3) implies that

*W*can be determined as

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 *v _{i}* in Eq. (2) will disappear.

*W*exists, which is required for the

*M*integral. All purely viscous incompressible fluids fall into the general class of Reiner–Rivlin fluids, which have a rheology given by

*I*is the

_{k}*k*th invariant of the tensor

*D*[13]. For a three-dimensional flow, the three invariants are

_{ij}*I*

_{3}or $\varphi 2\u22600$. Glen [17] describes ice as a Reiner–Rivlin fluid and concludes that the experimental data show sufficient scatter to warrant further study. Baker [18] reviews the subsequent experiments in determining the effects of the third invariant on the flow of ice and compares the results with his own experimental setup, which show that there is a significant dependence on

*I*

_{3}. Still there appears to be little evidence that $\varphi 2\u22600$ in ice. Such a fluid, as Glen [17] notes, would be susceptible to swelling or contraction in the direction perpendicular to the plane of shear. Schoof and Clarke [19] exploit this generation of deviatoric normal stress and use a Reiner–Rivlin fluid to model subglacial flutes by way of a secondary transverse flow. Here, we show that only Reiner–Rivlin fluids that are independent of

*I*

_{3}with $\varphi 2=0$ have a flow potential

*W*, unless

Therefore, our analysis primarily applies for ice modeled using a power-law rheology for ice, such as Glen's law.

## Analysis

*M*in polar coordinates and adopt a circular path of radius

*r*(

*ds*=

*rdθ*) as

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 *v _{r}*,

*v*, and

_{θ}*σ*. The antiplane terms, quantities with a subscript

_{rθ}*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*.

### Nye Solution.

*H*and density

*ρ*) far away. By adding a uniform tensile stress

_{i}*σ*

_{0}to the mass of ice, we transform our problem and apply a tensile stress $\sigma rr(r=a)=\sigma 0\u2212pw=\Delta p$ at the channel wall and a traction free condition at infinity. We are able to do this without changing the problem because of incompressibility and the pressure independence of Glen's law. Thus, the boundary conditions are

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

*W*can be written as

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 $dvr/dr\u21920$ as *r* → *∞*. This means that the flow potential *W* also decays to zero in the far field. Although *v _{r}r* is a constant as

*r*→

*∞*, the stress $\sigma rr(r\u2192\u221e)=0$ (boundary condition) and, therefore,

*M*= 0 as

*r*→

*∞*.

*M*integral around the channel must also be zero. Using the boundary condition $\sigma rr(r=a)=\Delta p$, the mass conservation, and the expression for the flow potential, we can write Eq. (7) as

*θ*and is therefore equal to zero. Taking care with the absolute value term, we can rearrange the integrand to find

*r*=

*b*on the outer edge in Fig. 2. Following the same method, and using $D\theta \theta =vr/r$ from the geometry and $D\theta \theta =\u2212Drr$ from the mass conservation, we have that

*M*integral around the outer edge of the domain is no longer zero. Since the rate of volume change of any ring is zero, i.e., $2\pi rvr=constant$, we can related the radial velocity at $vr(r=a)$ to $vr(r=a)$ as

which corroborates [22] and, also, reduces to the standard result as *b* → *∞*.

### Antiplane Shear.

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 [4] and Glen [26] 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. [27] 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 [28]. 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 [2] subglacial channel. Meyer et al. [29] 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.

*p*for the radial stress, we find that

*W*. Along

*r*=

*a*with superimposed antiplane shear, we have that

*W*

_{r}_{=}

*into Eq. (11), we find that*

_{a}#### Evaluation of M in the Far Field.

*M*by evaluating the contour as

*r*→

*∞*. Starting from Eq. (5), we now cancel all the in-plane terms, as they must decay as

*r*→

*∞*. This leaves us with

*M*gives

*W*as

*r*→

*∞*gives

*M*as

*r*→

*∞*reduces to

using the half-angle formula for $cos(\theta )$. In order for *M* to be finite, this integral must be zero because of the *r*^{2} 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.

#### Far-Field Perturbation.

*ε*is an unknown small parameter and

*f*(

*y*,

*z*),

*g*(

*y*,

*z*), and

*h*(

*y*,

*z*) are unknown functions. To first order in

*ε*, the effective strain rate is given as

*ε*, and we are able to ignore the velocities

*v*and

*w*. Thus, we concentrate on writing an equation for

*h*(

*y*,

*z*). The antiplane stresses are given by

*ε*give

*h*, we find that

*k*is some unknown integer. The term proportional to $sin(k\theta \u0302)$ can be ignored due to symmetry, and the positive solutions for

*λ*can be disregarded as they are singular as

*r*→

*∞*. The term that decays the slowest is

*k*=−1, and therefore, we have

*M*integral in the far field, we find that

*B*into

*ε*as the unknown far-field amplitude. We know that

*λ*= −1, and therefore, we have that

*I*(

_{M}*n*) is given by

*χ*(

*θ*) from Eq. (13), we can integrate over

*θ*numerically to find

*I*(

_{M}*n*). For

*n*= 3, i.e., Glen's law for ice, we find that

This integral should be negative because in $vx=\gamma \u02d9fary+\epsilon h(y,z)$, the antiplane velocity should be less than the boundary value as it approaches the edge.

*ε*, we need to know the behavior of

*v*for small

_{z}*r*, which requires solving the fully coupled (in-plane and antiplane) partial differential equation. Simultaneously, the path independence of the

*M*integral relates the perturbation integral in the far field to the integral around the channel, i.e., Eq. (12). This gives

Thus, the perturbation amplitude *ε* is related to the unknown strain rates at the edge of the channel. Furthermore, we can see that if *v _{x}* = 0 (or if

*v*is a function of

_{x}*r*only), then this integral reduces to the integral for the Nye solution and

*ε*= 0.

#### Nondimensional Equations and Numerics.

*M*in far-field nondimensionally. The natural length scale is the channel radius, and therefore, we write

*r*=

*aR*, where

*R*is the nondimensional radial coordinate (using capital letters to denote nondimensional variables). We proceed by using the boundary conditions to scale the stress and velocity. For the in-plane components of stress, Δ

*p*is a sensible scaling. This leads to a scaling for the in-plane velocity, by dimensional arguments alone, that is reminiscent of the Nye solution (Eq. (8)), i.e., $vr=Aa\Delta pnVr$. From the antiplane boundary conditions, we can see that $vx=\gamma \u02d9faraVx$. It is evident immediately that the in-plane and antiplane velocities do not scale in the same manner, and thus, a logical control parameter for the system is their ratio

*S*, which is

and represents the importance of antiplane shearing to in-plane creep closure [29]. Using *S* we can rewrite the antiplane velocity as $vx=Aa\Delta pnSVx$.

*M*as

*S*≫ 1 and the appropriate scaling for

*ε*is

*κ*is now a constant related to the unknown strain rates at the edge of the channel. When

*S*≪ 1, we scale

*ε*using the Nye strain rate as $\epsilon =\kappa A\Delta pna2$. From this definition of

*ε*, we can see that

*M*in the far field is either

Although we cannot solve for the constant *κ* analytically, we can determine its value by numerical simulations. We implement the numerical method described in Ref. [29] in the existing finite element software comsol [30]. 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* ∼ *S*^{1∕}* ^{n}* for

*S*≪ 1 and

*M*∼

*S*

^{1∕}

*for*

^{n}*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 *V _{r}* 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. [29].

*S*. When there is very little antiplane motion as compared to in-plane straining, the creep closure velocity is given by the Nye solution (Eq. (8)), which is written nondimensionally as

*S*≫ 1, the effective strain rate scales as

*D*∼

_{E}*S*. The radial force balance gives that the averaged creep closure velocity around the channel is given as

where more details are provided in Ref. [29]. 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 [4] and Glen [26] 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 [31] 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.

## Conclusions

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 $S\u2261\gamma \u02d9far/(A\Delta pn)$: For small antiplane velocity relative to in-plane creep closure, *ε* is a constant and the *M* integral approaches zero as *M* ∼ *S*^{1∕}* ^{n}* for vanishingly small

*S*. In the other limit, where antiplane straining dominates,

*M*grows as $M\u223cS(n+1)/n$. 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 $Vr\u223c\u2212S(n\u22121)/n$ due to a decrease in ice viscosity. Thus,

*M*provides a succinct description of the processes affecting channel closure with superimposed antiplane shear.

## Acknowledgment

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.

## Nomenclature

*a*=channel radius

*A*=ice softness

*D*=_{E}effective strain rate, $DijDij/2$

*D*=_{ij}strain rate field

*g*=acceleration due to gravity

*H*=ice thickness

*I*=_{i}strain rate tensor invariants

- $Ik$ =
invariant of strain rate tensor

*ℓ*=dimension, i.e., two or three

*m*=homogeneity degree,

*W*=_{s}*σ*/_{ij}ε_{ij}*m**M*=path-independent integral

*n*=rheological exponent, i.e., 3 for Glen's law

*p*=isotropic ice pressure

*p*=_{f}fluid pressure within channel

*S*=strain rate ratio, $\gamma \u02d9far/(A\Delta pn)$

*u*=_{i}displacement field

*v*=_{i}velocity field

*W*=flow potential

*W*=_{s}strain energy density

*x*=_{i}position field

- $\gamma \u02d9far$ =
far field antiplane strain rate

- Δ
*p*= effective pressure

*ε*=_{ij}strain field

*ρ*=_{i}ice density

*σ*=_{ij}stress field

*σ*_{0}=ice overburden pressure,

*ρ*_{i}gH- $\sigma \u2032ij$ =
deviatoric stress field,

*σ*+_{ij}*pδ*_{ij}*τ*=_{E}effective stress, $\sigma \u2032ij\sigma \u2032ij/2$

### Appendix A: Proof of the Path Independence of *M*

*M*integral as written in Eq. (2) with

*m*= 1 + 1/

*n*, i.e.,

*M*integral is path independent if the terms inside the area integral are zero. That is, if

*δ*= 2, and thus

_{kk}*W*can be related to the stress and strain rates as

*W*terms, we can use the chain rule to relate spatial derivatives on

*W*to derivatives of strain as

*W*is written to depend symmetrically on

*D*, with constraint

_{ij}*D*= 0, then

_{kk}*∂D*/

_{ij}*∂x*because

_{k}*∂D*/

_{ll}*∂x*= 0. Now, using the symmetry of the stress tensor, we have that

_{k}*M*integral, written as

is path independent.

### Appendix B: Flow Potential for a Reiner–Rivlin Fluid

*I*

_{1}=

*D*= 0) for which

_{kk}where *σ _{ij}* is the stress tensor,

*D*is the symmetric part of the velocity gradient tensor,

_{ij}*p*is the isotropic pressure, and

*I*is the

_{k}*k*th invariant of the tensor

*D*[13]. By writing the stress as an isotropic matrix function of

_{ij}*D*, assuming a symmetric dependence on

_{ij}*D*and

_{ij}*D*, 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

_{ji}*D*with coefficients that are functions of the invariants of

_{ij}*D*. For an incompressible fluid with isotropic pressure, this reduces to Eq. (B1).

_{ij}*dW*is a perfect differential of

*σ*. This requires that

_{ij}dD_{ij}where the isotropic pressure *p* is independent of the strain rate. This shows that the deviatoric stress also satisfies the Maxwell reciprocity.

*i*,

*j*,

*k*, and

*l*, the only way for this condition to be satisfied for all flows is for

*W*. If we start with the invariants of

*D*, given as

_{ij}*W*to exist, and to be a perfect differential of

*σ*as is required to write Eq. (B4), is that

_{ij}dD_{ij}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 $\varphi 1$ is a function of *I*_{2} solely and $\varphi 2=0$. This is the rheological structure of Glen's law in glaciology, and therefore, a flow potential *W* will always exist.