0
Research Papers

The Path-Independent M Integral Implies the Creep Closure of Englacial and Subglacial Channels OPEN ACCESS

[+] Author and Article Information
Colin R. Meyer

John A. Paulson School of Engineering and Applied Sciences,
Harvard University,
Cambridge, MA 02138
e-mail: colinrmeyer@gmail.com

John W. Hutchinson

Professor
Fellow ASME
John A. Paulson School of Engineering and Applied Sciences,
Harvard University,
Cambridge, MA 02138
email: hutchinson@husm.harvard.edu

James R. Rice

Professor
Fellow ASME
John A. Paulson School of Engineering and Applied Sciences,
Harvard University,
Cambridge, MA 02138;
Department of Earth and Planetary Sciences,
Harvard University,
Cambridge, MA 02138
e-mail: rice@seas.harvard.edu

1Corresponding author.

Manuscript received September 7, 2016; final manuscript received September 21, 2016; published online October 18, 2016. Editor: Yonggang Huang.

J. Appl. Mech 84(1), 011006 (Oct 18, 2016) (9 pages) Paper No: JAM-16-1433; doi: 10.1115/1.4034828 History: Received September 07, 2016; Revised September 21, 2016

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.

FIGURES IN THIS ARTICLE
<>

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 [13]. 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.

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 [5], and independently Knowles and Sternberg [6], introduced a path-independent integral for linear elastic solid mechanics, which Budiansky and Rice [7] called the 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α=xα+γxαanduα(x¯)=uα(x¯)+(12)γuα(x¯)

That is, coordinates 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 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

Γ{Wsyαnα+σαknk[fα(x¯)yβuαxβ]} ds

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

M=Γ{Wsxαnασαknkxβuαxβ} ds

as written by Knowles and Sternberg [6].

Budiansky and Rice [7] 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 [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 Display Formula

(1)M=Γ{Wsxknk+σiknk[(m2m)uixjuixj]} ds

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 [11]. 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 Display Formula

(2)M=Γ{Wxknk+σiknk[(m2m)vixjvixj]} ds

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

WDij=σij+pδij,   W=1mσijDij,anddW=σijdDij

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(vixj+vjxi)withvkxk=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 Display Formula

(3)DE=AτEn

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 [12]. Now, relating the deviatoric stress and strain rate tensors using Eq. (3) implies that Display Formula

(4)Dij=AτEn1(σij+pδij)

Using this rheology, the flow potential W can be determined as

W=σijdDij=nn+1σijDij=1mσijDij=2nn+1A1nDE1+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 [13]. For a three-dimensional flow, the three invariants are

I1=Dkk,   I2=12DlmDml,   I3=13DlmDmnDnl

For a fluid that is incompressible

I1=0

Truesdell and Noll [14] assert that there is little experimental evidence for fluids with ϕ20. This assertion is based on Markovitz and Williamson [15], who found the polymeric data collected by Padden and Dewitt [16] to be incompatible with ϕ20. In glaciology, it is not fully resolved whether Glen's law should be expanded to include a dependence on I3 or ϕ20. 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 I3. Still there appears to be little evidence that ϕ20 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 I3 with ϕ2=0 have a flow potential W, unless

ϕ1I3=ϕ2I2

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

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 Display Formula

(5)M=02π{Wr(n1n+1)[σrrvr+σrθvθ+σrxvx][σrrvrr+σrθvθr+σrxvxr]r} rdθ

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 [4] 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)=σ0pw=Δ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)=σ0pw=Δpandσ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

Display Formula

(6)M=02π{Wr2(n1n+1)σrrvrrσrrdvrdrr2} dθ

The mass conservation gives that

dvrdr+vrr=0

and, therefore, we can simplify the terms in Eq. (6) as Display Formula

(7)M=02π{Wr2+2σrrvrn+1r} dθ

For in-plane, axisymmetric motion, the flow potential W can be written as

W=2nn+1A1nDE1+1/n=2nn+1A1n|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/dr0 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

M=2an+102π{nA1n|vr|1+1/na1/n+Δpvr} dθ=0

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 Display Formula

(8)vr(a)=Aa(Δpn)n

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 Display Formula

(9)M=2na2n+102π{A1n|vr(a)a|1+1n+Δpvr(a)na} dθ=2nb2n+102πA1n|vr(b)b|1+1n dθ

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 [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 [2325]. 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.

The problem setup for including antiplane shear follows Fig. 2 with the additional boundary conditions Display Formula

(10)σrx|r=a=0andvx|r=γ˙farrcos(θ)

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

M=02π{Wr2(n1n+1)σrrvrrσrrvrrr2} dθ

Using the mass conservation, we can write the radial velocity derivative as

vrr=vrr1rvθθ

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 Display Formula

(11)M=02π{Wa2+2Δpn+1vra} dθ

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+1A1/nDE1+1/n|r=a=2nn+1A1n[(vrr)2+(12avxθ)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 Display Formula

(12)M=2a2Δpn+102π{nA1nΔp[(vrr)2+(12avxθ)2]n+12n+vra} dθ

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

M=02π{Wr2(n1n+1)σrxvxrσrxvxrr2} dθ

We start by assuming that the antiplane fields are given by the boundary conditions, we have that

vx=γ˙farrcos(θ)andσrx=A1/n(γ˙far/2)1/ncos(θ)

which is a state of constant strain rate. Inserting these expressions into M gives

M=02π{W2nn+1(2A)1/nγ˙far1+1/ncos2(θ)}r2 dθ

Evaluating W as r gives

W|r=2nn+1A1/n(γ˙far/2)1+1/n

Thus, the integral for M as r reduces to

M=nn+1(2A)1/nγ˙far1+1/n02π{12cos2(θ)}r2 dθ=0

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),andvz=ε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εγ˙farhy+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=A1/nDE(1n)/n(γ˙far2+ε2hy)σzx=A1/nDE(1n)/nε2hz

Inserting the effective strain rate and linearizing the stress about ε give

σyx=[A2γ˙far]1/n(1+εnγ˙farhy)σzx=[A2γ˙far]1/nεγ˙farhz

We then insert the linearized stress into the force balance equation

σxyy+σyzz=0

and we find that

1n2hy2+2hz2=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+z2andθ̂=arctan(zη)

we can write

h=r̂λf(θ̂)

Inserting this ansatz into Laplace's equation for h, we find that

f=Bcos(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 Display Formula

(13)h=Br̂cos(θ̂)=Br(cos(arctan(tan(θ)n))1+(n1)cos2(θ))

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

M=ε[γ˙far2A]1/nrλ+102π{λcos(θ)χsin(θ)χ}2(n1)n+1[λcos3(θ)χ12cos(θ)sin(2θ)χ]n1n+1(λ+1)cos(θ)χ+2λcos(θ)χ dθ

where we have incorporated the unknown factor B into ε as the unknown far-field amplitude. We know that λ = −1, and therefore, we have that

M=ε[γ˙far2A]1/n02π2(n1)n+1[cos3(θ)χ+sin(θ)cos2(θ)χ]+2cos(θ)χ dθ

which is of the form

M=ε[γ˙far2A]1/nIM(n)

where IM(n) is given by

02π2(n1)n+1[cos3(θ)χ+sin(θ)cos2(θ)χ]+2cos(θ)χ dθ

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 Display Formula

(14)M=2a2Δpn+102π{nA1nΔp[(vrr)2+(12avxθ)2]n+12n+vra} dθ=ε[γ˙far2A]1/nIM(n)

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 [29]. 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 Display Formula

(15)M=2n+102π{n[(VrR)2+S24(Vxθ)2]n+12n+Vr} dθ=εa2AΔpn(S2)1/nIM(n)

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

M=κIM(n)21/nS1/n   (S1)orκIM(n)21/nS(n+1)/n   (S1)

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

κ=2.45   (S1)orκ=86.97   (S1)

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. [29].

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=nnR

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 Display Formula

(16)VrS(n1)/n

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.

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 MS(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 VrS(n1)/n 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 =

    channel radius

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

M=Γ{Wxknkσiknk[(n1n+1)vi+xjvixj]} ds

where Γ is a closed path, as shown in Fig. 1. We use the divergence theorem to write this integral as

M=Ωxk(Wxkσik[(n1n+1)vi+xjvixj]) dA

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

xk(Wxkσik[(n1n+1)vi+xjvixj])=0

Taking the derivatives, we find that Display Formula

(A1)Wxkxk+Wxkxkσikxk[(n1n+1)vi+xjvixj]σik[(n1n+1)vixk+xjxkvixj+xj2vixjxk]=0

The equilibrium equations are

σikxk=0

which allow us to cancel the third group of terms in Eq. (A1) and, therefore, write

Wxkxk+Wxkxkσik[(n1n+1)vixk+xjxkvixj+xj2vixjxk]=0

The derivatives of coordinates lead to Kronecker delta functions as

xixj=δij

where in two dimensions, δkk = 2, and thus

Wxkxk+2Wσik[(n1n+1)vixk+δjkvixj+xj2vixjxk]=0

We can simplify this expression slightly by multiplying through by the stress and combining like terms as

Wxkxk+2W2nn+1σikvixkσikxj2vixjxk=0

The strain rate energy density function W can be related to the stress and strain rates as

W=nn+1σijDij=nn+1σijvixj

This allows us to write

Wxkxk+2W2Wσij2vixkxjxk=0

Canceling the 2 W terms, we can use the chain rule to relate spatial derivatives on W to derivatives of strain as

WDijDijxkxkσij2vixkxjxk=0

If W is written to depend symmetrically on Dij, with constraint Dkk = 0, then

σij+pδij=WDij

and thus, the integrand is

σijDijxkxkσij2vixkxjxk=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

σij2vixkxjxkσij2vixkxjxk=0

which is zero. Thus, we have shown that the M integral, written as

M=Γ{Wxknkσiknk[(n1n+1)vi+xjvixj]} ds

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 Display Formula

(B1)σij=pδij+ϕ1(I2,I3)Dij+ϕ2(I2,I3)DinDnj

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 [13]. 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 Display Formula

(B2)σijDkl=σklDij

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

(ϕ1I3ϕ2I2)(DijDknDnlDklDinDnj)=0

Now since DijDknDnlDklDinDnj for all i, j, k, and l, the only way for this condition to be satisfied for all flows is for Display Formula

(B3)ϕ1I3=ϕ2I2

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

I1=Dkk,   I2=12DlmDml,   I3=13DlmDmnDnl

we can write the relationship between the flow potential and the stress as Display Formula

(B4)σij+pδij=WDij=WI2I2Dij+WI3I3Dij

Now using the facts

DijDkl=δikδjl,   I2Dkl=Dkl,andI3Dkl=DknDnl

we have that Display Formula

(B5)σij+pδij=WDij=WI2Dij+WI3DinDnj

from which we can see that Eq. (B5) is Reiner–Rivlin fluid with

ϕ1=WI2andϕ2=WI3

Thus, a requirement for W to exist, and to be a perfect differential of σijdDij as is required to write Eq. (B4), is that

2WI2I3=2WI3I2

which implies that

ϕ1I3=ϕ2I2

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.

Lliboutry, L. , 1968, “ General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers,” J. Glaciol., 7, pp. 21–58.
Röthlisberger, H. , 1972, “ Water Pressure in Intra- and Subglacial Channels,” J. Glaciol., 11(62), pp. 177–203.
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. [CrossRef] [PubMed]
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. [CrossRef]
Günther, W. , 1962, “ Über einige randintegrale der elastomechanik,” Abh. Braunschw. Wiss. Ges., 14, pp. 53–72.
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. [CrossRef]
Budiansky, B. , and Rice, J. R. , 1973, “ Conservation Laws and Energy-Release Rates,” ASME J. Appl. Mech., 40(1), pp. 201–203. [CrossRef]
Noether, E. , 1918, “ Invariante variations-probleme,” Nachr. Ges. Gottingen, Math. Phys. Klasse, pp. 235–257.
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.
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. [CrossRef]
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. [CrossRef]
Glen, J. W. , 1955, “ The Creep of Polycrystalline Ice,” Proc. R. Soc. London, Ser. A, 228(1175), pp. 519–538. [CrossRef]
Schowalter, W. R. , 1978, Mechanics of Non-Newtonian Fluids, Pergamon Press, Elmsford, NY.
Truesdell, C. , and Noll, W. , 1965, The Non-Linear Field Theories of Mechanics, 3rd ed., Springer, New York.
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. [CrossRef]
Padden, F. J. , and Dewitt, T. W. , 1954, “ Some Rheological Properties of Concentrated Polyisobutylene Solutions,” J. Appl. Phys., 25(9), pp. 1086–1091. [CrossRef]
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.
Baker, R. W. , 1987, The Physical Basis of Ice Sheet Modelling, IAHS, Washington, DC, pp. 7–16.
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. [CrossRef]
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. [CrossRef]
Cuffey, K. M. , and Paterson, W. S. B. , 2010, The Physics of Glaciers, 4th ed., Elsevier, Burlington, MA.
Evatt, G. W. , 2015, “ Röthlisberger Channels With Finite Ice Depth and Open Channel Flow,” Ann. Glaciol., 50(70), pp. 45–50. [CrossRef]
Schoof, C. , 2010, “ Ice-Sheet Acceleration Driven by Melt Supply Variability,” Nature, 468(7325), pp. 803–806. [CrossRef] [PubMed]
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. [CrossRef]
Hewitt, I. J. , 2013, “ Seasonal Changes in Ice Sheet Motion Due to Melt Water Lubrication,” Earth Planet. Sci. Lett., 371, pp. 16–25. [CrossRef]
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.
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. [CrossRef]
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. [CrossRef]
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. [CrossRef]
COMSOL, 2006, “ COMSOL Multiphysics: Version 4.3,” COMSOL, Inc., Stockholm, Sweden.
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. [CrossRef]
Copyright © 2017 by ASME
View article in PDF format.

References

Lliboutry, L. , 1968, “ General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers,” J. Glaciol., 7, pp. 21–58.
Röthlisberger, H. , 1972, “ Water Pressure in Intra- and Subglacial Channels,” J. Glaciol., 11(62), pp. 177–203.
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. [CrossRef] [PubMed]
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. [CrossRef]
Günther, W. , 1962, “ Über einige randintegrale der elastomechanik,” Abh. Braunschw. Wiss. Ges., 14, pp. 53–72.
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. [CrossRef]
Budiansky, B. , and Rice, J. R. , 1973, “ Conservation Laws and Energy-Release Rates,” ASME J. Appl. Mech., 40(1), pp. 201–203. [CrossRef]
Noether, E. , 1918, “ Invariante variations-probleme,” Nachr. Ges. Gottingen, Math. Phys. Klasse, pp. 235–257.
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.
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. [CrossRef]
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. [CrossRef]
Glen, J. W. , 1955, “ The Creep of Polycrystalline Ice,” Proc. R. Soc. London, Ser. A, 228(1175), pp. 519–538. [CrossRef]
Schowalter, W. R. , 1978, Mechanics of Non-Newtonian Fluids, Pergamon Press, Elmsford, NY.
Truesdell, C. , and Noll, W. , 1965, The Non-Linear Field Theories of Mechanics, 3rd ed., Springer, New York.
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. [CrossRef]
Padden, F. J. , and Dewitt, T. W. , 1954, “ Some Rheological Properties of Concentrated Polyisobutylene Solutions,” J. Appl. Phys., 25(9), pp. 1086–1091. [CrossRef]
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.
Baker, R. W. , 1987, The Physical Basis of Ice Sheet Modelling, IAHS, Washington, DC, pp. 7–16.
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. [CrossRef]
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. [CrossRef]
Cuffey, K. M. , and Paterson, W. S. B. , 2010, The Physics of Glaciers, 4th ed., Elsevier, Burlington, MA.
Evatt, G. W. , 2015, “ Röthlisberger Channels With Finite Ice Depth and Open Channel Flow,” Ann. Glaciol., 50(70), pp. 45–50. [CrossRef]
Schoof, C. , 2010, “ Ice-Sheet Acceleration Driven by Melt Supply Variability,” Nature, 468(7325), pp. 803–806. [CrossRef] [PubMed]
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. [CrossRef]
Hewitt, I. J. , 2013, “ Seasonal Changes in Ice Sheet Motion Due to Melt Water Lubrication,” Earth Planet. Sci. Lett., 371, pp. 16–25. [CrossRef]
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.
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. [CrossRef]
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. [CrossRef]
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. [CrossRef]
COMSOL, 2006, “ COMSOL Multiphysics: Version 4.3,” COMSOL, Inc., Stockholm, Sweden.
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. [CrossRef]

Figures

Grahic Jump Location
Fig. 1

Schematic of the path-independent M integral around a void. The path of integration is Γ and the vector normal to this integration path is ni. The down glacier component is x which points out of the page, i.e., antiplane, the in-plane coordinates are y (across glacier) and z (depth).

Grahic Jump Location
Fig. 2

Nye problem setup and boundary conditions. Inward-pointing arrows indicate the creep closure of the ice (Ur), and the eddy structures in the channel highlight the turbulence (turbulent pipe flow simulation by Feldmann and Wagner [20]; photo credit: Noel Fitzpatrick).

Grahic Jump Location
Fig. 3

Schematic for an idealized ice stream shear margin with a Röthlisberger [2] subglacial channel. Surface ice velocity increases away from the margin and is nearly stagnant in the ridge. Adapted from Ref. [29].

Grahic Jump Location
Fig. 4

Numerical evaluation of M at the edge of the channel Minner and in the far-field Mouter at R = 500. The scalings show that for S ≪ 1, M ∼ S1∕n and for S ≫ 1, M ∼ S(n+1)∕n. Numerical artifacts lead to negative values for Minner when S is small; therefore, these simulations are omitted but yield absolute values close to Mouter.

Grahic Jump Location
Fig. 5

Numerical simulations showing the nondimensional channel creep closure velocity |Vr| at the edge of the channel R = 1. For S ≪ 1, the antiplane straining is negligible and the creep closure velocity approaches the nondimensional Nye solution. When the antiplane motions are dominant, i.e., S ≫ 1, the creep closure velocity increases as S(n−1)/n.

Tables

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In