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 . 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.

Path-independent (conserved) integrals are important mathematical tools that are often employed in mechanics as a method of solution to the equations or as a supplemental constraint. Günther , and independently Knowles and Sternberg , introduced a path-independent integral for linear elastic solid mechanics, which Budiansky and Rice  called the M integral. Using the Noether  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α′=xα+γxα and uα′(x¯′)=uα(x¯)+(1−ℓ2)γuα(x¯)$
That is, coordinates xα and displacements uα are self-similarly scaled from the reference configuration . In the framework of linearized kinematics, the strain is equal to the symmetric part of the displacement ui gradient tensor as $εij=sym(∂ui/∂xj)$. We can then write a strain energy density Ws as a product of stresses σij and strains εij by Ws = σijεij/2. The general conserved integral resulting from the Noether procedure, with $yα=xα′−xα$ and $fα(x¯)=uα′(x¯′)−uα(x¯)$, is given by
where this integral is a line integral for two-dimensional problems and a surface integral for three-dimensional problems. In this way, the M integral for a linear elastic material in two dimensions with linearized kinematics is given as

as written by Knowles and Sternberg .

Budiansky and Rice  extended the earlier definitions of M to a generalized elastic material with a strain energy density Ws that is homogeneous of degree m in the strains εij, and therefore, Ws = σijεij/m. Unfortunately, the expression of Budiansky and Rice  for the generalized M integral contains an error. Whether it is typographical, conceptual, or due to the printing process is unknown. He and Hutchinson  give a correct expression (although without derivation) for the three-dimensional generalized M integral in a different geometry than used here or in Ref. : a closed 3D void or flat crack of axially symmetric shape, such that the stresses vary with z and $x2+y2$. Rice  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
(1)

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.

The generalized definition of M for a power-law nonlinear elastic solid is exactly equivalent to the definition for a power-law nonlinear viscous fluid, where the displacement ui is replaced by the velocity vi and the strain εij by the strain rate Dij . The strain rate is the symmetric part of the velocity gradient tensor, as $Dij=sym(∂vi/∂xj)$. Under the definition for a power-law viscous fluid, the elastic strain energy density Ws is replaced by a function called the flow potential W. Just as $dWs=σijdDij$ is an exact differential in generalized elasticity, dW is also an exact differential, satisfying $dW=σ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
(2)

We include a proof of the path independence of Eq. (2) in Appendix  A.

The flow potential W for an incompressible viscous fluid is given as
where $p (=−σkk/3)$ is the isotropic pressure (understood here as a Lagrange multiplier to enforce mass conservation, εkk = 0), and dW is an exact differential. The strain rates are defined in terms of derivatives of the velocity vi as
$Dij=12(∂vi∂xj+∂vj∂xi) with ∂vk∂xk=0$

where the first condition shows that Dij is symmetric, and by the second condition, the mass conservation for an incompressible substance Dkk = 0.

### Ice Rheology.

Here, we describe how the rheology sets the value for the parameter m. In glaciology, ice is commonly modeled as an incompressible viscous fluid with a nonlinear rheology called Glen's law
$DE=AτEn$
(3)
The effective strain rate DE 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 $σ′ij=σij+pδij$, and the subscript 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 . Now, relating the deviatoric stress and strain rate tensors using Eq. (3) implies that
$Dij=AτEn−1(σij+pδij)$
(4)
Using this rheology, the flow potential W can be determined as
$W=∫σijdDij=nn+1σijDij=1mσijDij=2nn+1A−1nDE1+1/n$

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.

In Appendix  B, we show the general class of incompressible viscous fluids for which a flow potential 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
$σij=−pδij+ϕ1(I2,I3)Dij+ϕ2(I2,I3)DinDnj$
where Ik is the kth invariant of the tensor Dij . For a three-dimensional flow, the three invariants are
For a fluid that is incompressible
$I1=0$
Truesdell and Noll  assert that there is little experimental evidence for fluids with $ϕ2≠0$. This assertion is based on Markovitz and Williamson , who found the polymeric data collected by Padden and Dewitt  to be incompatible with $ϕ2≠0$. In glaciology, it is not fully resolved whether Glen's law should be expanded to include a dependence on I3 or $ϕ2≠0$. Glen  describes ice as a Reiner–Rivlin fluid and concludes that the experimental data show sufficient scatter to warrant further study. Baker  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 I3. Still there appears to be little evidence that $ϕ2≠0$ in ice. Such a fluid, as Glen  notes, would be susceptible to swelling or contraction in the direction perpendicular to the plane of shear. Schoof and Clarke  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 I3 with $ϕ2=0$ have a flow potential W, unless
$∂ϕ1∂I3=∂ϕ2∂I2$

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

## Analysis

To analyze the creep closure of an drainage channel, it is convenient to write M in polar coordinates and adopt a circular path of radius r (ds = rdθ) as
(5)

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 σ. 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.

### Nye Solution.

Nye  derived the rate of closure of a circular channel in a Glen rheology subject to a stress $σrr(r=a)=−pw$ (water pressure) applied at the channel and the stress $σrr(r→∞)=−σ0=−ρigH$ (overburden ice pressure for a glacier of height H and density ρi) far away. By adding a uniform tensile stress σ0 to the mass of ice, we transform our problem and apply a tensile stress $σrr(r=a)=σ0−pw=Δ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
$σrr(r=a)=σ0−pw=Δp and σrr(r→∞)=0$

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

(6)
The mass conservation gives that
$dvrdr+vrr=0$
and, therefore, we can simplify the terms in Eq. (6) as
(7)
For in-plane, axisymmetric motion, the flow potential W can be written as
$W=2nn+1A−1nDE1+1/n=2nn+1A−1n|dvrdr|1+1/n$

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→0$ 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 $σrr(r→∞)=0$ (boundary condition) and, therefore, M = 0 as r.

Thus, the M integral around the channel must also be zero. Using the boundary condition $σrr(r=a)=Δp$, the mass conservation, and the expression for the flow potential, we can write Eq. (7) as
Now, due to axisymmetry, the integrand must be independent of θ and is therefore equal to zero. Taking care with the absolute value term, we can rearrange the integrand to find
$vr(a)=−Aa(Δpn)n$
(8)

which is the Nye solution for the creep closure rate at the edge of the channel [4,21].

This analysis can be easily extended to the case where the outer boundary is finite, i.e., where r = b on the outer edge in Fig. 2. Following the same method, and using $Dθθ=vr/r$ from the geometry and $Dθθ=−Drr$ from the mass conservation, we have that
(9)
Now, the 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πrvr=constant$, we can related the radial velocity at $vr(r=a)$ to $vr(r=a)$ as
$avr(a)=bvr(b)$
Solving for the creep closure rate at the edge of the channel $vr(r=a)$ using Eq. (9), we find that the finite domain Nye solution is given as
$vr(a)=−Aa[1−(a/b)2/n]n(Δpn)n$

which corroborates  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 . 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.

The problem setup for including antiplane shear follows Fig. 2 with the additional boundary conditions
$σrx|r=a=0 and vx|r→∞=γ˙farr cos(θ)$
(10)

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.

We start by evaluating M at the channel. We cancel the in-plane shear terms (e.g., σ) and antiplane terms in Eq. (5), due to the boundary conditions in Eq. (10). Thus, we arrive at
Using the mass conservation, we can write the radial velocity derivative as
$∂vr∂r=−vrr−1r∂vθ∂θ$
The second term in this expression can be omitted because it will always be zero when integrated around a closed loop for a constant boundary stress. If we insert the effective pressure Δp for the radial stress, we find that
(11)
which is nearly identical to Eq. (7) in symbols. The difference is in the value of W. Along r = a with superimposed antiplane shear, we have that
$Wr=a=2nn+1A−1/nDE1+1/n|r=a=2nn+1A−1n[(∂vr∂r)2+(12a∂vx∂θ)2]1+1/n$
where the shear term $∂vx/∂r$ was not included due to the boundary condition specified in Eq. (10). Inserting this expression for Wr=a into Eq. (11), we find that
(12)

#### Evaluation of M in the Far Field.

Following the derivation of the Nye solution, we seek to determine the value of 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
We start by assuming that the antiplane fields are given by the boundary conditions, we have that
$vx=γ˙farr cos(θ) and σrx=A−1/n(γ˙far/2)1/n cos(θ)$
which is a state of constant strain rate. Inserting these expressions into M gives
Evaluating W as r gives
$W|r→∞=2nn+1A−1/n(γ˙far/2)1+1/n$
Thus, the integral for M as r reduces to

using the half-angle formula for $cos(θ)$. 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.

#### Far-Field Perturbation.

In the far field, we consider a perturbation from the constant antiplane strain rate boundary condition. We write the far-field velocities in Cartesian coordinates as
$vx=γ˙fary+εh(y,z), vy=εg(y,z), and vz=εf(y,z)$
Here, ε 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
$DE=γ˙far21+2εγ˙far∂h∂y+O(ε2)$
Consequently, the antiplane problem in the far field is coupled to the in-plane motions through ε, 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
$σyx=A−1/nDE(1−n)/n(γ˙far2+ε2∂h∂y)σzx=A−1/nDE(1−n)/nε2∂h∂z$
Inserting the effective strain rate and linearizing the stress about ε give
$σyx=[A2γ˙far]1/n(1+εnγ˙far∂h∂y)σzx=[A2γ˙far]1/nεγ˙far∂h∂z$
We then insert the linearized stress into the force balance equation
$∂σxy∂y+∂σyz∂z=0$
and we find that
$1n∂2h∂y2+∂2h∂z2=0$
By making the transformation, $η=ny$, we can write this equation as
$∇̂2h=0$
i.e., Laplace's equation with $∇̂=(∂η,∂z)$. If we return to polar coordinates, now with
$r̂=η2+z2 and θ̂=arctan(zη)$
we can write
$h=r̂λf(θ̂)$
Inserting this ansatz into Laplace's equation for h, we find that
$f=B cos(kθ̂) and λ=±k$
where k is some unknown integer. The term proportional to $sin(kθ̂)$ 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
$h=Br̂cos(θ̂)=Br( cos(arctan( tan(θ)n))1+(n−1) cos2(θ))$
(13)
which we refer to as
$h(r,θ)=Brχ(θ)$
Now inserting the perturbation ansatz $h=εrλχ(θ)$ into the M integral in the far field, we find that
where we have incorporated the unknown factor B into ε as the unknown far-field amplitude. We know that λ = −1, and therefore, we have that
which is of the form
$M=ε[γ˙far2A]1/nIM(n)$
where IM(n) is given by
Inserting our expression for χ(θ) from Eq. (13), we can integrate over θ numerically to find IM(n). For n = 3, i.e., Glen's law for ice, we find that
$IM(3)=−0.113…$

This integral should be negative because in $vx=γ˙fary+εh(y,z)$, the antiplane velocity should be less than the boundary value as it approaches the edge.

Now, to determine the far-field amplitude factor ε, we need to know the behavior of vz for small 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
(14)

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.

We now write the expression for 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ΔpnVr$. From the antiplane boundary conditions, we can see that $vx=γ˙faraVx$. 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
$S=γ˙farAΔpn$

and represents the importance of antiplane shearing to in-plane creep closure . Using S we can rewrite the antiplane velocity as $vx=AaΔpnSVx$.

Using these scalings, we can rewrite M as
$M=a2AΔpn+1M̂$
where we immediately drop the variable hats. Thus, Eq. (14) gives
(15)
When the far field of the domain is dominated by antiplane shear, we have that S ≫ 1 and the appropriate scaling for ε is
$ε=κγ˙fara2$
where κ 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 $ε=κAΔ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.  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. .

We now describe the evolution of the creep closure velocity as a function of 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
$Vr=−n−nR$
When the deformation is dominated by antiplane motion, i.e., S ≫ 1, the effective strain rate scales as DE ∼ S. The radial force balance gives that the averaged creep closure velocity around the channel is given as
$Vr∼−S(n−1)/n$
(16)

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.

## 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≡γ˙far/(AΔpn)$: 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 $M∼S(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∼−S(n−1)/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 =

•
• A =

ice softness

•
• DE =

effective strain rate, $DijDij/2$

•
• Dij =

strain rate field

•
• g =

acceleration due to gravity

•
• H =

ice thickness

•
• Ii =

strain rate tensor invariants

•
• $Ik$ =

invariant of strain rate tensor

•
• =

dimension, i.e., two or three

•
• m =

homogeneity degree, Ws = σijεij/m

•
• M =

path-independent integral

•
• n =

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

•
• p =

isotropic ice pressure

•
• pf =

fluid pressure within channel

•
• S =

strain rate ratio, $γ˙far/(AΔpn)$

•
• ui =

displacement field

•
• vi =

velocity field

•
• W =

flow potential

•
• Ws =

strain energy density

•
• xi =

position field

•
• $γ˙far$ =

far field antiplane strain rate

•
• Δp =

effective pressure

•
• εij =

strain field

•
• ρi =

ice density

•
• σij =

stress field

•
• σ0 =

ice overburden pressure, ρigH

•
• $σ′ij$ =

deviatoric stress field, σij + ij

•
• τE =

effective stress, $σ′ijσ′ij/2$

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

Here, we start with the two-dimensional generalized M integral as written in Eq. (2) with m = 1 + 1/n, i.e.,
where Γ is a closed path, as shown in Fig. 1. We use the divergence theorem to write this integral as
The M integral is path independent if the terms inside the area integral are zero. That is, if
$∂∂xk(Wxk−σik[(n−1n+1)vi+xj∂vi∂xj])=0$
Taking the derivatives, we find that
$∂W∂xkxk+W∂xk∂xk−∂σik∂xk[(n−1n+1)vi+xj∂vi∂xj]−σik[(n−1n+1)∂vi∂xk+∂xj∂xk∂vi∂xj+xj∂2vi∂xj∂xk]=0$
(A1)
The equilibrium equations are
$∂σik∂xk=0$
which allow us to cancel the third group of terms in Eq. (A1) and, therefore, write
$∂W∂xkxk+W∂xk∂xk−σik[(n−1n+1)∂vi∂xk+∂xj∂xk∂vi∂xj+xj∂2vi∂xj∂xk]=0$
The derivatives of coordinates lead to Kronecker delta functions as
$∂xi∂xj=δij$
where in two dimensions, δkk = 2, and thus
$∂W∂xkxk+2W−σik[(n−1n+1)∂vi∂xk+δjk∂vi∂xj+xj∂2vi∂xj∂xk]=0$
We can simplify this expression slightly by multiplying through by the stress and combining like terms as
$∂W∂xkxk+2W−2nn+1σik∂vi∂xk−σikxj∂2vi∂xj∂xk=0$
The strain rate energy density function W can be related to the stress and strain rates as
$W=nn+1σijDij=nn+1σij∂vi∂xj$
This allows us to write
$∂W∂xkxk+2W−2W−σij∂2vi∂xk∂xjxk=0$
Canceling the 2 W terms, we can use the chain rule to relate spatial derivatives on W to derivatives of strain as
$∂W∂Dij∂Dij∂xkxk−σij∂2vi∂xk∂xjxk=0$
If W is written to depend symmetrically on Dij, with constraint Dkk = 0, then
$σij+pδij=∂W∂Dij$
and thus, the integrand is
$σij∂Dij∂xkxk−σij∂2vi∂xk∂xjxk=0$
where the isotropic components of the stress tensor cancel when multiplied by ∂Dij/∂xk because ∂Dll/∂xk = 0. Now, using the symmetry of the stress tensor, we have that
$σij∂2vi∂xk∂xjxk−σij∂2vi∂xk∂xjxk=0$
which is zero. Thus, we have shown that the M integral, written as

is path independent.

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

A Reiner–Rivlin fluid is an incompressible fluid (I1 = Dkk = 0) for which
$σij=−pδij+ϕ1(I2,I3)Dij+ϕ2(I2,I3)DinDnj$
(B1)

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).

Here, we ask what is the most general fluid rheology that still posses a flow potential, i.e., where dW is a perfect differential of σijdDij. This requires that
$∂σij∂Dkl=∂σkl∂Dij$
(B2)
which is called Maxwell reciprocity. For an incompressible fluid, we include pressure as a Lagrange multiplier to enforce the mass conservation and write that
$∂(σ′ij+pδij)∂Dkl=∂(σ′kl+pδij)∂Dij$

where the isotropic pressure p is independent of the strain rate. This shows that the deviatoric stress also satisfies the Maxwell reciprocity.

If we insert the Reiner–Rivlin fluid rheology into Eq. (B2), we find that
$(∂ϕ1∂I3−∂ϕ2∂I2)(DijDknDnl−DklDinDnj)=0$
Now since $DijDknDnl≠DklDinDnj$ for all i, j, k, and l, the only way for this condition to be satisfied for all flows is for
$∂ϕ1∂I3=∂ϕ2∂I2$
(B3)
This condition also arises from an equality requirement of the mixed partial derivatives of the flow potential W. If we start with the invariants of Dij, given as
we can write the relationship between the flow potential and the stress as
$σij+pδij=∂W∂Dij=∂W∂I2∂I2∂Dij+∂W∂I3∂I3∂Dij$
(B4)
Now using the facts
we have that
$σij+pδij=∂W∂Dij=∂W∂I2Dij+∂W∂I3DinDnj$
(B5)
from which we can see that Eq. (B5) is Reiner–Rivlin fluid with
$ϕ1=∂W∂I2 and ϕ2=∂W∂I3$
Thus, a requirement for W to exist, and to be a perfect differential of σijdDij as is required to write Eq. (B4), is that
$∂2W∂I2∂I3=∂2W∂I3∂I2$
which implies that
$∂ϕ1∂I3=∂ϕ2∂I2$

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 $ϕ1$ is a function of I2 solely and $ϕ2=0$. This is the rheological structure of Glen's law in glaciology, and therefore, a flow potential W will always exist.

## References

References
1.
Lliboutry
,
L.
,
1968
, “
General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers
,”
J. Glaciol.
,
7
, pp.
21
58
.
2.
Röthlisberger
,
H.
,
1972
, “
Water Pressure in Intra- and Subglacial Channels
,”
J. Glaciol.
,
11
(
62
), pp.
177
203
.
3.
Kamb
,
B.
,
Raymond
,
C. F.
,
Harrison
,
W. D.
,
Engelhardt
,
H.
,
Echelmeyer
,
K. A.
,
Humphrey
,
N.
,
Brugman
,
M. M.
, and
Pfeffer
,
T.
,
1985
, “
Glacier Surge Mechanism: 1982–1983 Surge of Variegated Glacier, Alaska
,”
Science
,
227
(
4686
), pp.
469
479
.
4.
Nye
,
J. F.
,
1953
, “
The Flow Law of Ice From Measurements in Glacier Tunnels, Laboratory Experiments and the Jungfraufirn Borehole Experiment
,”
Proc. R. Soc. London, Ser. A
,
219
(
1139
), pp.
477
489
.
5.
Günther
,
W.
,
1962
, “
Über einige randintegrale der elastomechanik
,”
Abh. Braunschw. Wiss. Ges.
,
14
, pp.
53
72
.
6.
Knowles
,
J. K.
, and
Sternberg
,
E.
,
1972
, “
On a Class of Conservation Laws in Linearized and Finite Elastostatics
,”
Arch. Ration. Mech. Anal.
,
44
, pp.
187
211
.
7.
Budiansky
,
B.
, and
Rice
,
J. R.
,
1973
, “
Conservation Laws and Energy-Release Rates
,”
ASME J. Appl. Mech.
,
40
(
1
), pp.
201
203
.
8.
Noether
,
E.
,
1918
, “
Invariante variations-probleme
,”
Nachr. Ges. Gottingen, Math. Phys. Klasse
, pp.
235
257
.
9.
Rice
,
J. R.
,
1985
, “
Conserved Integrals and Energetic Forces
,”
Eshelby Memorial Symposium: Fundamentals of Deformation and Fracture
,
B. A.
Bilby
,
K. J.
Miller
, and
J. R.
Willis
, eds., Sheffield, UK, Apr. 2–5, Cambridge University Press, Cambridge, UK, pp.
33
56
.
10.
He
,
M. Y.
, and
Hutchinson
,
J. W.
,
1981
, “
The Penny-Shaped Crack and the Plane Strain Crack in an Infinite Body of Power-Law Material
,”
ASME J. Appl. Mech.
,
48
(
4
), pp.
830
840
.
11.
Ben Amar
,
M.
, and
Rice
,
J. R.
,
2002
, “
Exact Results With the J-Integral Applied to Free-Boundary Flows
,”
J. Fluid Mech.
,
461
, pp.
321
341
.
12.
Glen
,
J. W.
,
1955
, “
The Creep of Polycrystalline Ice
,”
Proc. R. Soc. London, Ser. A
,
228
(
1175
), pp.
519
538
.
13.
Schowalter
,
W. R.
,
1978
,
Mechanics of Non-Newtonian Fluids
,
Pergamon Press
,
Elmsford, NY
.
14.
Truesdell
,
C.
, and
Noll
,
W.
,
1965
,
The Non-Linear Field Theories of Mechanics
, 3rd ed.,
Springer
,
New York
.
15.
Markovitz
,
H.
, and
Williamson
,
R. B.
,
1957
, “
Normal Stress Effect in Polyisobutylene Solutions. I. Measurements in a Cone and Plate Instrument
,”
Trans. Soc. Rheol.
,
1
(
1
), pp.
25
36
.
16.
,
F. J.
, and
Dewitt
,
T. W.
,
1954
, “
Some Rheological Properties of Concentrated Polyisobutylene Solutions
,”
J. Appl. Phys.
,
25
(
9
), pp.
1086
1091
.
17.
Glen
,
J. W.
,
1958
,
The Flow Law of Ice: A Discussion of the Assumptions Made in Glacier Theory, Their Experimental Foundations and Consequences
,
IAHS
,
Washington, DC
, Vol.
47
, pp.
171
183
.
18.
Baker
,
R. W.
,
1987
,
The Physical Basis of Ice Sheet Modelling
,
IAHS
,
Washington, DC
, pp.
7
16
.
19.
Schoof
,
C.
, and
Clarke
,
G. K. C.
,
2008
, “
A Model for Spiral Flows in Basal Ice and the Formation of Subglacial Flutes Based on a Reiner–Rivlin Rheology for Glacial Ice
,”
J. Geophys. Res.
,
113
(
B5
), pp.
1
12
.
20.
Feldmann
,
D.
, and
Wagner
,
C.
,
2012
, “
Direct Numerical Simulation of Fully Developed Turbulent and Oscillatory Pipe Flows at Reτ = 1440
,”
J. Turbul.
,
13
, pp.
1
28
.
21.
Cuffey
,
K. M.
, and
Paterson
,
W. S. B.
,
2010
,
The Physics of Glaciers
, 4th ed.,
Elsevier
,
Burlington, MA
.
22.
Evatt
,
G. W.
,
2015
, “
Röthlisberger Channels With Finite Ice Depth and Open Channel Flow
,”
Ann. Glaciol.
,
50
(
70
), pp.
45
50
.
23.
Schoof
,
C.
,
2010
, “
Ice-Sheet Acceleration Driven by Melt Supply Variability
,”
Nature
,
468
(
7325
), pp.
803
806
.
24.
Bartholomous
,
T. C.
,
Anderson
,
R. S.
, and
Anderson
,
S. P.
,
2011
, “
Growth and Collapse of the Distributed Subglacial Hydrologic System of Kennicott Glacier, Alaska, USA, and Its Effects on Basal Motion
,”
J. Glaciol.
,
57
(
206
), pp.
985
1002
.
25.
Hewitt
,
I. J.
,
2013
, “
Seasonal Changes in Ice Sheet Motion Due to Melt Water Lubrication
,”
Earth Planet. Sci. Lett.
,
371
, pp.
16
25
.
26.
Glen
,
J. W.
,
1956
, “
Measurement of the Deformation of Ice in a Tunnel at the Foot of an Ice Fall
,”
J. Glaciol.
,
2
(
20
), pp.
735
745
.
27.
Perol
,
T.
,
Rice
,
J. R.
,
Platt
,
J. D.
, and
Suckale
,
J.
,
2015
, “
Subglacial Hydrology and Ice Stream Margin Locations
,”
J. Geophys. Res.
,
120
(
7
), pp.
1352
1368
.
28.
Vogel
,
S. W.
,
Tulaczyk
,
S.
,
Kamb
,
B.
,
Engelhardt
,
H.
,
Carsey
,
F. D.
,
Behar
,
A. E.
,
Lane
,
A. L.
, and
Joughin
,
I.
,
2005
, “
Subglacial Conditions During and After Stoppage of an Antarctic Ice Stream: Is Reactivation Imminent?
,”
Geophys. Res. Lett.
,
32
(
14
), pp.
1
4
.
29.
Meyer
,
C. R.
,
Fernandes
,
M. C.
,
Creyts
,
T. T.
, and
Rice
,
J. R.
,
2016
, “
Effects of Ice Deformation on Röthlisberger Channels and Implications for Transitions in Subglacial Hydrology
,”
J. Glaciol.
,
62
(
234
), pp.
750
762
.
30.
COMSOL
,
2006
, “
COMSOL Multiphysics: Version 4.3
,” COMSOL, Inc., Stockholm, Sweden.
31.
Rice
,
J. R.
, and
Tracey
,
D. M.
,
1969
, “
On the Ductile Enlargement of Voids in Triaxial Stress Fields
,”
J. Mech. Phys. Solids
,
17
(
3
), pp.
201
217
.