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

Meyer, Colin R.; Hutchinson, John W.; Rice, James R.

doi:10.1115/1.4034828

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.

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_{α} and u_{α}^{'} (\underline{x}') = u_{α} (\underline{x}) + (1 - \frac{ℓ}{2}) γ u_{α} (\underline{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 u_i gradient tensor as

ε_{i j} = sym (\partial u_{i} / \partial x_{j})

⁠. We can then write a strain energy density W_s as a product of stresses σ_ij and strains ε_ij by W_s = σ_ijε_ij/2. The general conserved integral resulting from the Noether procedure, with

y_{α} = x_{α}^{'} - x_{α}

and

f_{α} (\underline{x}) = u_{α}^{'} (\underline{x}') - u_{α} (\underline{x})

⁠, is given by

\oint_{Γ} {W_{s} y_{α} n_{α} + σ_{α k} n_{k} [f_{α} (\underline{x}) - y_{β} \frac{\partial u_{α}}{\partial x_{β}}]} d s

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 = \oint_{Γ} {W_{s} x_{α} n_{α} - σ_{α k} n_{k} x_{β} \frac{\partial u_{α}}{\partial x_{β}}} d s

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 W_s that is homogeneous of degree m in the strains ε_ij, and therefore, 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

\sqrt{x^{2} + y^{2}}

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

M = \oint_{Γ} {W_{s} x_{k} n_{k} + σ_{i k} n_{k} [(\frac{m - 2}{m}) u_{i} - x_{j} \frac{\partial u_{i}}{\partial x_{j}}]} d s

(1)

where n_i is the unit outer normal to the closed contour Γ, and s is arc length anticlockwise around the path, such that n₁ds = dx₂ and n₂ds = −dx₁, where the tensor subscripts correspond as (x, y, z) = (x₁, x₂, x₃) and (u_x, u_y, u_z) = (u₁, u₂, u₃). Figure 1 shows the domain, the coordinate system, and a path of integration about a void.

Fig. 1

View large Download slide

Schematic of the path-independent M integral around a void. The path of integration is Γ and the vector normal to this integration path is n_i. 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).

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 u_i is replaced by the velocity v_i and the strain ε_ij by the strain rate D_ij [11]. The strain rate is the symmetric part of the velocity gradient tensor, as

D_{i j} = sym (\partial v_{i} / \partial x_{j})

⁠. Under the definition for a power-law viscous fluid, the elastic strain energy density W_s is replaced by a function called the flow potential W. Just as

d W_{s} = σ_{i j} d D_{i j}

is an exact differential in generalized elasticity, dW is also an exact differential, satisfying

d W = σ_{i j} d D_{i j}

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

M = \oint_{Γ} {W x_{k} n_{k} + σ_{i k} n_{k} [(\frac{m - 2}{m}) v_{i} - x_{j} \frac{\partial v_{i}}{\partial x_{j}}]} d s

(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

\frac{\partial W}{\partial D_{i j}} = σ_{i j} + p δ_{i j}, W = \frac{1}{m} σ_{i j} D_{i j}, and d W = σ_{i j} d D_{i j}

where

p (= - σ_{k k} / 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 v_i as

D_{i j} = \frac{1}{2} (\frac{\partial v_{i}}{\partial x_{j}} + \frac{\partial v_{j}}{\partial x_{i}}) with \frac{\partial v_{k}}{\partial x_{k}} = 0

where the first condition shows that D_ij is symmetric, and by the second condition, the mass conservation for an incompressible substance D_kk = 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

D_{E} = A τ_{E}^{n}

(3)

The effective strain rate D_E is a function of the second invariant of the strain rate tensor, i.e.,

D_{E} = \sqrt{D_{i j} D_{i j} / 2}

⁠, τ_E is defined in the same way for the deviatoric stress tensor

σ'_{i j} = σ_{i j} + p δ_{i j}

⁠, 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. implies that

D_{i j} = A τ_{E}^{n - 1} (σ_{i j} + p δ_{i j})

(4)

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

W = \int σ_{i j} d D_{i j} = \frac{n}{n + 1} σ_{i j} D_{i j} = \frac{1}{m} σ_{i j} D_{i j} = \frac{2 n}{n + 1} A^{- \frac{1}{n}} D_{E}^{1 + 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 v_i 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

σ_{i j} = - p δ_{i j} + ϕ_{1} (I_{2}, I_{3}) D_{i j} + ϕ_{2} (I_{2}, I_{3}) D_{i n} D_{n j}

where I_k is the kth invariant of the tensor D_ij [13]. For a three-dimensional flow, the three invariants are

I_{1} = D_{k k}, I_{2} = \frac{1}{2} D_{l m} D_{m l}, I_{3} = \frac{1}{3} D_{l m} D_{m n} D_{n l}

For a fluid that is incompressible

I_{1} = 0

Truesdell and Noll [14] assert that there is little experimental evidence for fluids with

ϕ_{2} \neq 0

⁠. This assertion is based on Markovitz and Williamson [15], who found the polymeric data collected by Padden and Dewitt [16] to be incompatible with

ϕ_{2} \neq 0

⁠. In glaciology, it is not fully resolved whether Glen's law should be expanded to include a dependence on I₃ or

ϕ_{2} \neq 0

⁠. 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₃. Still there appears to be little evidence that

ϕ_{2} \neq 0

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₃ with

ϕ_{2} = 0

have a flow potential W, unless

\frac{\partial ϕ_{1}}{\partial I_{3}} = \frac{\partial ϕ_{2}}{\partial I_{2}}

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

\begin{matrix} M = \int_{0}^{2 π} {W r - (\frac{n - 1}{n + 1}) [σ_{r r} v_{r} + σ_{r θ} v_{θ} + σ_{r x} v_{x}] \\ - [σ_{r r} \frac{\partial v_{r}}{\partial r} + σ_{r θ} \frac{\partial v_{θ}}{\partial r} + σ_{r x} \frac{\partial v_{x}}{\partial r}] r} r d θ \end{matrix}

(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 v_r, v_θ, and σ_rθ. The antiplane terms, quantities with a subscript x, represent motion into and out of page as a function of only the in-plane coordinates (using the standard glaciology coordinate system with z vertical, y across glacier, and x down glacier). We consider a very long channel with constant ice thickness and, in this way, we can reduce a three-dimensional problem to two dimensions where the quantities are homogeneous along x.

Nye Solution.

Nye [4] derived the rate of closure of a circular channel in a Glen rheology subject to a stress

σ_{r r} (r = a) = - p_{w}

(water pressure) applied at the channel and the stress

σ_{r r} (r \to \infty) = - σ_{0} = - ρ_{i} g H

(overburden ice pressure for a glacier of height H and density ρ_i) far away. By adding a uniform tensile stress σ₀ to the mass of ice, we transform our problem and apply a tensile stress

σ_{r r} (r = a) = σ_{0} - p_{w} = Δ 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

σ_{r r} (r = a) = σ_{0} - p_{w} = Δ p and σ_{r r} (r \to \infty) = 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

Fig. 2

View large Download slide

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

M = \int_{0}^{2 π} {W r^{2} - (\frac{n - 1}{n + 1}) σ_{r r} v_{r} r - σ_{r r} \frac{d v_{r}}{d r} r^{2}} d θ

(6)

The mass conservation gives that

\frac{d v_{r}}{d r} + \frac{v_{r}}{r} = 0

and, therefore, we can simplify the terms in Eq. as

M = \int_{0}^{2 π} {W r^{2} + \frac{2 σ_{r r} v_{r}}{n + 1} r} d θ

(7)

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

W = \frac{2 n}{n + 1} A^{- \frac{1}{n}} D_{E}^{1 + 1 / n} = \frac{2 n}{n + 1} A^{- \frac{1}{n}} {| \frac{d v_{r}}{d r} |}^{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 $d v_{r} / d r \to 0$ as r → ∞. This means that the flow potential W also decays to zero in the far field. Although v_rr is a constant as r → ∞, the stress $σ_{r r} (r \to \infty) = 0$ (boundary condition) and, therefore, M = 0 as r → ∞.

Thus, the M integral around the channel must also be zero. Using the boundary condition

σ_{r r} (r = a) = Δ p

⁠, the mass conservation, and the expression for the flow potential, we can write Eq. as

M = \frac{2 a}{n + 1} \int_{0}^{2 π} {n A^{- \frac{1}{n}} {| - v_{r} |}^{1 + 1 / n} a^{- 1 / n} + Δ p v_{r}} 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

v_{r} (a) = - A a {(\frac{Δ p}{n})}^{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_{θ θ} = v_{r} / r

from the geometry and

D_{θ θ} = - D_{r r}

from the mass conservation, we have that

\begin{matrix} M = \frac{2 n a^{2}}{n + 1} \int_{0}^{2 π} {A^{- \frac{1}{n}} {| - \frac{v_{r} (a)}{a} |}^{1 + \frac{1}{n}} + \frac{Δ p v_{r} (a)}{n a}} d θ \\ = \frac{2 n b^{2}}{n + 1} \int_{0}^{2 π} A^{- \frac{1}{n}} {| - \frac{v_{r} (b)}{b} |}^{1 + \frac{1}{n}} d θ \end{matrix}

(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 π r v_{r} = constant

⁠, we can related the radial velocity at

v_{r} (r = a)

to

v_{r} (r = a)

as

a v_{r} (a) = b v_{r} (b)

Solving for the creep closure rate at the edge of the channel

v_{r} (r = a)

using Eq. , we find that the finite domain Nye solution is given as

v_{r} (a) = - \frac{A a}{{[1 - {(a / b)}^{2 / n}]}^{n}} {(\frac{Δ p}{n})}^{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 [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.

Fig. 3

View large Download slide

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

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

{σ_{r x} |}_{r = a} = 0 and {v_{x} |}_{r \to \infty} = {\dot{γ}}_{far} r 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., σ_rθ) and antiplane terms in Eq. , due to the boundary conditions in Eq. . Thus, we arrive at

M = \int_{0}^{2 π} {W r^{2} - (\frac{n - 1}{n + 1}) σ_{r r} v_{r} r - σ_{r r} \frac{\partial v_{r}}{\partial r} r^{2}} d θ

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

\frac{\partial v_{r}}{\partial r} = - \frac{v_{r}}{r} - \frac{1}{r} \frac{\partial v_{θ}}{\partial θ}

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

M = \int_{0}^{2 π} {W a^{2} + \frac{2 Δ p}{n + 1} v_{r} a} d θ

(11)

which is nearly identical to Eq. in symbols. The difference is in the value of W. Along r = a with superimposed antiplane shear, we have that

\begin{matrix} W_{r = a} = {\frac{2 n}{n + 1} A^{- 1 / n} D_{E}^{1 + 1 / n} |}_{r = a} \\ = \frac{2 n}{n + 1} A^{- \frac{1}{n}} {[\sqrt{{(\frac{\partial v_{r}}{\partial r})}^{2} + {(\frac{1}{2 a} \frac{\partial v_{x}}{\partial θ})}^{2}}]}^{1 + 1 / n} \end{matrix}

where the shear term

\partial v_{x} / \partial r

was not included due to the boundary condition specified in Eq. . Inserting this expression for W_r₌_a into Eq. , we find that

\begin{matrix} M = \frac{2 a^{2} Δ p}{n + 1} \int_{0}^{2 π} {\frac{n A^{- \frac{1}{n}}}{Δ p} [{(\frac{\partial v_{r}}{\partial r})}^{2} \\ + {(\frac{1}{2 a} \frac{\partial v_{x}}{\partial θ})}^{2}]^{\frac{n + 1}{2 n}} + \frac{v_{r}}{a}} d θ \end{matrix}

(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. , we now cancel all the in-plane terms, as they must decay as r → ∞. This leaves us with

M = \int_{0}^{2 π} {W r^{2} - (\frac{n - 1}{n + 1}) σ_{r x} v_{x} r - σ_{r x} \frac{\partial v_{x}}{\partial r} r^{2}} d θ

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

v_{x} = {\dot{γ}}_{far} r cos (θ) and σ_{r x} = A^{- 1 / n} {({\dot{γ}}_{far} / 2)}^{1 / n} cos (θ)

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

M = \int_{0}^{2 π} {W - \frac{2 n}{n + 1} {(2 A)}^{- 1 / n} {\dot{γ}}_{far}^{1 + 1 / n} {cos}^{2} (θ)} r^{2} d θ

Evaluating W as r → ∞ gives

{W |}_{r \to \infty} = \frac{2 n}{n + 1} A^{- 1 / n} {({\dot{γ}}_{far} / 2)}^{1 + 1 / n}

Thus, the integral for M as r → ∞ reduces to

M = \frac{n}{n + 1} {(2 A)}^{- 1 / n} {\dot{γ}}_{far}^{1 + 1 / n} \int_{0}^{2 π} {1 - 2 {cos}^{2} (θ)} r^{2} d θ = 0

using the half-angle formula for $cos (θ)$ ⁠. In order for M to be finite, this integral must be zero because of the r² 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

v_{x} = {\dot{γ}}_{far} y + ε h (y, z), v_{y} = ε g (y, z), and v_{z} = ε 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

D_{E} = \frac{{\dot{γ}}_{far}}{2} \sqrt{1 + \frac{2 ε}{{\dot{γ}}_{far}} \frac{\partial h}{\partial 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

\begin{matrix} σ_{y x} = A^{- 1 / n} D_{E}^{(1 - n) / n} (\frac{{\dot{γ}}_{far}}{2} + \frac{ε}{2} \frac{\partial h}{\partial y}) \\ σ_{z x} = A^{- 1 / n} D_{E}^{(1 - n) / n} \frac{ε}{2} \frac{\partial h}{\partial z} \end{matrix}

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

\begin{matrix} σ_{y x} = {[\frac{A}{2 {\dot{γ}}_{far}}]}^{1 / n} (1 + \frac{ε}{n {\dot{γ}}_{far}} \frac{\partial h}{\partial y}) \\ σ_{z x} = {[\frac{A}{2 {\dot{γ}}_{far}}]}^{1 / n} \frac{ε}{{\dot{γ}}_{far}} \frac{\partial h}{\partial z} \end{matrix}

We then insert the linearized stress into the force balance equation

\frac{\partial σ_{x y}}{\partial y} + \frac{\partial σ_{y z}}{\partial z} = 0

and we find that

\frac{1}{n} \frac{\partial^{2} h}{\partial y^{2}} + \frac{\partial^{2} h}{\partial z^{2}} = 0

By making the transformation,

η = \sqrt{n} y

⁠, we can write this equation as

{\hat{\nabla}}^{2} h = 0

i.e., Laplace's equation with

\hat{\nabla} = (\partial_{η}, \partial_{z})

⁠. If we return to polar coordinates, now with

\hat{r} = \sqrt{η^{2} + z^{2}} and \hat{θ} = \arctan (\frac{z}{η})

we can write

h = {\hat{r}}^{λ} f (\hat{θ})

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

f = B cos (k \hat{θ}) and λ = \pm k

where k is some unknown integer. The term proportional to

sin (k \hat{θ})

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 = \frac{B}{\hat{r}} cos (\hat{θ}) = \frac{B}{r} (\frac{cos (\arctan (\frac{tan (θ)}{\sqrt{n}}))}{\sqrt{1 + (n - 1) {cos}^{2} (θ)}})

(13)

which we refer to as

h (r, θ) = \frac{B}{r} χ (θ)

Now inserting the perturbation ansatz

h = ε r^{λ} χ (θ)

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

\begin{matrix} M = ε {[\frac{{\dot{γ}}_{far}}{2 A}]}^{1 / n} r^{λ + 1} \int_{0}^{2 π} {λ cos (θ) χ - sin (θ) χ'} \\ - \frac{2 (n - 1)}{n + 1} [λ {cos}^{3} (θ) χ - \frac{1}{2} cos (θ) sin (2 θ) χ'] \\ - \frac{n - 1}{n + 1} (λ + 1) cos (θ) χ + 2 λ cos (θ) χ d θ \end{matrix}

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

\begin{matrix} M = ε {[\frac{{\dot{γ}}_{far}}{2 A}]}^{1 / n} \int_{0}^{2 π} \frac{2 (n - 1)}{n + 1} [{cos}^{3} (θ) χ \\ + sin (θ) {cos}^{2} (θ) χ'] + 2 cos (θ) χ d θ \end{matrix}

which is of the form

M = ε {[\frac{{\dot{γ}}_{far}}{2 A}]}^{1 / n} I_{M} (n)

where I_M(n) is given by

\int_{0}^{2 π} \frac{2 (n - 1)}{n + 1} [{cos}^{3} (θ) χ + sin (θ) {cos}^{2} (θ) χ'] + 2 cos (θ) χ d θ

Inserting our expression for χ(θ) from Eq. , we can integrate over θ numerically to find I_M(n). For n = 3, i.e., Glen's law for ice, we find that

I_{M} (3) = - 0.113 \dots

This integral should be negative because in $v_{x} = {\dot{γ}}_{far} y + ε 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 v_z 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. . This gives

\begin{matrix} M = \frac{2 a^{2} Δ p}{n + 1} \int_{0}^{2 π} {\frac{n A^{- \frac{1}{n}}}{Δ p} {[{(\frac{\partial v_{r}}{\partial r})}^{2} + {(\frac{1}{2 a} \frac{\partial v_{x}}{\partial θ})}^{2}]}^{\frac{n + 1}{2 n}} + \frac{v_{r}}{a}} d θ \\ = ε {[\frac{{\dot{γ}}_{far}}{2 A}]}^{1 / n} I_{M} (n) \end{matrix}

(14)

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_x 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. ), i.e.,

v_{r} = A a Δ p^{n} V_{r}

⁠. From the antiplane boundary conditions, we can see that

v_{x} = {\dot{γ}}_{far} a V_{x}

⁠. 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 = \frac{{\dot{γ}}_{far}}{A Δ p^{n}}

and represents the importance of antiplane shearing to in-plane creep closure [29]. Using S we can rewrite the antiplane velocity as $v_{x} = A a Δ p^{n} S V_{x}$ ⁠.

Using these scalings, we can rewrite M as

M = a^{2} A Δ p^{n + 1} \hat{M}

where we immediately drop the variable hats. Thus, Eq. gives

\begin{matrix} M = \frac{2}{n + 1} \int_{0}^{2 π} {n {[{(\frac{\partial V_{r}}{\partial R})}^{2} + \frac{S^{2}}{4} {(\frac{\partial V_{x}}{\partial θ})}^{2}]}^{\frac{n + 1}{2 n}} + V_{r}} d θ \\ = \frac{ε}{a^{2} A Δ p^{n}} {(\frac{S}{2})}^{1 / n} I_{M} (n) \end{matrix}

(15)

When the far field of the domain is dominated by antiplane shear, we have that S ≫ 1 and the appropriate scaling for ε is

ε = κ {\dot{γ}}_{far} a^{2}

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 Δ p^{n} a^{2}

⁠. From this definition of ε, we can see that M in the far field is either

M = \frac{κ I_{M} (n)}{2^{1 / n}} S^{1 / n} (S ≪ 1) or \frac{κ I_{M} (n)}{2^{1 / n}} S^{(n + 1) / n} (S ≫ 1)

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∕ⁿ for S ≪ 1 and M ∼ S^1∕ⁿ for S ≫ 1, are clearly visible. The best-fit value of κ determined from the simulations is given by

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.

View large Download slide

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

κ = - 2.45 (S ≪ 1) or κ = - 86.97 (S ≫ 1)

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

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. ), which is written nondimensionally as

V_{r} = - \frac{n^{- n}}{R}

When the deformation is dominated by antiplane motion, i.e., 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

V_{r} \sim - S^{(n - 1) / n}

(16)

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

Fig. 5

View large Download slide

Numerical simulations showing the nondimensional channel creep closure velocity $| V_{r} |$ 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}$ ⁠.

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 \equiv {\dot{γ}}_{far} / (A Δ p^{n})$ ⁠: For small antiplane velocity relative to in-plane creep closure, ε is a constant and the M integral approaches zero as M ∼ S^1∕ⁿ for vanishingly small S. In the other limit, where antiplane straining dominates, M grows as $M \sim 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 $V_{r} \sim - 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 =: channel radius
A =: ice softness
D_E =: effective strain rate, $\sqrt{D_{i j} D_{i j} / 2}$
D_ij =: strain rate field
g =: acceleration due to gravity
H =: ice thickness
I_i =: strain rate tensor invariants
$I_{k}$ =: 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, ${\dot{γ}}_{far} / (A Δ p^{n})$
u_i =: displacement field
v_i =: velocity field
W =: flow potential
W_s =: strain energy density
x_i =: position field
${\dot{γ}}_{far}$ =: far field antiplane strain rate
Δp =: effective pressure
ε_ij =: strain field
ρ_i =: ice density
σ_ij =: stress field
σ₀ =: ice overburden pressure, ρ_igH
$σ'_{i j}$ =: deviatoric stress field, σ_ij + pδ_ij
τ_E =: effective stress, $\sqrt{σ'_{i j} σ'_{i j} / 2}$

Appendix A: Proof of the Path Independence of M

Here, we start with the two-dimensional generalized M integral as written in Eq. with m = 1 + 1/n, i.e.,

M = \oint_{Γ} {W x_{k} n_{k} - σ_{i k} n_{k} [(\frac{n - 1}{n + 1}) v_{i} + x_{j} \frac{\partial v_{i}}{\partial x_{j}}]} d s

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

M = \oint_{Ω} \frac{\partial}{\partial x_{k}} (W x_{k} - σ_{i k} [(\frac{n - 1}{n + 1}) v_{i} + x_{j} \frac{\partial v_{i}}{\partial x_{j}}]) d A

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

\frac{\partial}{\partial x_{k}} (W x_{k} - σ_{i k} [(\frac{n - 1}{n + 1}) v_{i} + x_{j} \frac{\partial v_{i}}{\partial x_{j}}]) = 0

Taking the derivatives, we find that

\begin{matrix} \frac{\partial W}{\partial x_{k}} x_{k} + W \frac{\partial x_{k}}{\partial x_{k}} - \frac{\partial σ_{i k}}{\partial x_{k}} [(\frac{n - 1}{n + 1}) v_{i} + x_{j} \frac{\partial v_{i}}{\partial x_{j}}] \\ - σ_{i k} [(\frac{n - 1}{n + 1}) \frac{\partial v_{i}}{\partial x_{k}} + \frac{\partial x_{j}}{\partial x_{k}} \frac{\partial v_{i}}{\partial x_{j}} + x_{j} \frac{\partial^{2} v_{i}}{\partial x_{j} \partial x_{k}}] = 0 \end{matrix}

(A1)

The equilibrium equations are

\frac{\partial σ_{i k}}{\partial x_{k}} = 0

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

\frac{\partial W}{\partial x_{k}} x_{k} + W \frac{\partial x_{k}}{\partial x_{k}} - σ_{i k} [(\frac{n - 1}{n + 1}) \frac{\partial v_{i}}{\partial x_{k}} + \frac{\partial x_{j}}{\partial x_{k}} \frac{\partial v_{i}}{\partial x_{j}} + x_{j} \frac{\partial^{2} v_{i}}{\partial x_{j} \partial x_{k}}] = 0

The derivatives of coordinates lead to Kronecker delta functions as

\frac{\partial x_{i}}{\partial x_{j}} = δ_{i j}

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

\frac{\partial W}{\partial x_{k}} x_{k} + 2 W - σ_{i k} [(\frac{n - 1}{n + 1}) \frac{\partial v_{i}}{\partial x_{k}} + δ_{j k} \frac{\partial v_{i}}{\partial x_{j}} + x_{j} \frac{\partial^{2} v_{i}}{\partial x_{j} \partial x_{k}}] = 0

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

\frac{\partial W}{\partial x_{k}} x_{k} + 2 W - \frac{2 n}{n + 1} σ_{i k} \frac{\partial v_{i}}{\partial x_{k}} - σ_{i k} x_{j} \frac{\partial^{2} v_{i}}{\partial x_{j} \partial x_{k}} = 0

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

W = \frac{n}{n + 1} σ_{i j} D_{i j} = \frac{n}{n + 1} σ_{i j} \frac{\partial v_{i}}{\partial x_{j}}

This allows us to write

\frac{\partial W}{\partial x_{k}} x_{k} + 2 W - 2 W - σ_{i j} \frac{\partial^{2} v_{i}}{\partial x_{k} \partial x_{j}} x_{k} = 0

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

\frac{\partial W}{\partial D_{i j}} \frac{\partial D_{i j}}{\partial x_{k}} x_{k} - σ_{i j} \frac{\partial^{2} v_{i}}{\partial x_{k} \partial x_{j}} x_{k} = 0

If W is written to depend symmetrically on D_ij, with constraint D_kk = 0, then

σ_{i j} + p δ_{i j} = \frac{\partial W}{\partial D_{i j}}

and thus, the integrand is

σ_{i j} \frac{\partial D_{i j}}{\partial x_{k}} x_{k} - σ_{i j} \frac{\partial^{2} v_{i}}{\partial x_{k} \partial x_{j}} x_{k} = 0

where the isotropic components of the stress tensor cancel when multiplied by ∂D_ij/∂x_k because ∂D_ll/∂x_k = 0. Now, using the symmetry of the stress tensor, we have that

σ_{i j} \frac{\partial^{2} v_{i}}{\partial x_{k} \partial x_{j}} x_{k} - σ_{i j} \frac{\partial^{2} v_{i}}{\partial x_{k} \partial x_{j}} x_{k} = 0

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

M = \oint_{Γ} {W x_{k} n_{k} - σ_{i k} n_{k} [(\frac{n - 1}{n + 1}) v_{i} + x_{j} \frac{\partial v_{i}}{\partial x_{j}}]} d s

is path independent.

Appendix B: Flow Potential for a Reiner–Rivlin Fluid

A Reiner–Rivlin fluid is an incompressible fluid (I₁ = D_kk = 0) for which

σ_{i j} = - p δ_{i j} + ϕ_{1} (I_{2}, I_{3}) D_{i j} + ϕ_{2} (I_{2}, I_{3}) D_{i n} D_{n j}

(B1)

where σ_ij is the stress tensor, D_ij is the symmetric part of the velocity gradient tensor, p is the isotropic pressure, and I_k is the kth invariant of the tensor D_ij [13]. By writing the stress as an isotropic matrix function of D_ij, assuming a symmetric dependence on D_ij and D_ji, 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 D_ij with coefficients that are functions of the invariants of D_ij. 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 σ_ijdD_ij. This requires that

\frac{\partial σ_{i j}}{\partial D_{k l}} = \frac{\partial σ_{k l}}{\partial D_{i j}}

(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

\frac{\partial (σ'_{i j} + p δ_{i j})}{\partial D_{k l}} = \frac{\partial (σ'_{k l} + p δ_{i j})}{\partial D_{i j}}

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. , we find that

(\frac{\partial ϕ_{1}}{\partial I_{3}} - \frac{\partial ϕ_{2}}{\partial I_{2}}) (D_{i j} D_{k n} D_{n l} - D_{k l} D_{i n} D_{n j}) = 0

Now since

D_{i j} D_{k n} D_{n l} \neq D_{k l} D_{i n} D_{n j}

for all i, j, k, and l, the only way for this condition to be satisfied for all flows is for

\frac{\partial ϕ_{1}}{\partial I_{3}} = \frac{\partial ϕ_{2}}{\partial I_{2}}

(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 D_ij, given as

I_{1} = D_{k k}, I_{2} = \frac{1}{2} D_{l m} D_{m l}, I_{3} = \frac{1}{3} D_{l m} D_{m n} D_{n l}

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

σ_{i j} + p δ_{i j} = \frac{\partial W}{\partial D_{i j}} = \frac{\partial W}{\partial I_{2}} \frac{\partial I_{2}}{\partial D_{i j}} + \frac{\partial W}{\partial I_{3}} \frac{\partial I_{3}}{\partial D_{i j}}

(B4)

Now using the facts

\frac{\partial D_{i j}}{\partial D_{k l}} = δ_{i k} δ_{j l}, \frac{\partial I_{2}}{\partial D_{k l}} = D_{k l}, and \frac{\partial I_{3}}{\partial D_{k l}} = D_{k n} D_{n l}

we have that

σ_{i j} + p δ_{i j} = \frac{\partial W}{\partial D_{i j}} = \frac{\partial W}{\partial I_{2}} D_{i j} + \frac{\partial W}{\partial I_{3}} D_{i n} D_{n j}

(B5)

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

ϕ_{1} = \frac{\partial W}{\partial I_{2}} and ϕ_{2} = \frac{\partial W}{\partial I_{3}}

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

\frac{\partial^{2} W}{\partial I_{2} \partial I_{3}} = \frac{\partial^{2} W}{\partial I_{3} \partial I_{2}}

which implies that

\frac{\partial ϕ_{1}}{\partial I_{3}} = \frac{\partial ϕ_{2}}{\partial I_{2}}

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 I₂ 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

1.

Lliboutry

,

L.

,

1968

, “

General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers

,”

J. Glaciol.

,

7

, pp.

21

–

58

.

Google Scholar

Crossref

2.

Röthlisberger

,

H.

,

1972

, “

Water Pressure in Intra- and Subglacial Channels

,”

J. Glaciol.

,

11

(

62

), pp.

177

–

203

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

PubMed

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

7.

Budiansky

,

B.

, and

Rice

,

J. R.

,

1973

, “

Conservation Laws and Energy-Release Rates

,”

ASME J. Appl. Mech.

,

40

(

1

), pp.

201

–

203

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

12.

Glen

,

J. W.

,

1955

, “

The Creep of Polycrystalline Ice

,”

Proc. R. Soc. London, Ser. A

,

228

(

1175

), pp.

519

–

538

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

16.

Padden

,

F. J.

, and

Dewitt

,

T. W.

,

1954

, “

Some Rheological Properties of Concentrated Polyisobutylene Solutions

,”

J. Appl. Phys.

,

25

(

9

), pp.

1086

–

1091

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

23.

Schoof

,

C.

,

2010

, “

Ice-Sheet Acceleration Driven by Melt Supply Variability

,”

Nature

,

468

(

7325

), pp.

803

–

806

.

Google Scholar

Crossref

PubMed

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

.

Google Scholar

Crossref

25.

Hewitt

,

I. J.

,

2013

, “

Seasonal Changes in Ice Sheet Motion Due to Melt Water Lubrication

,”

Earth Planet. Sci. Lett.

,

371

, pp.

16

–

25

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

Introduction

Ice Rheology.

Analysis

Nye Solution.

Antiplane Shear.

Evaluation of the M Integral at the Channel.

Evaluation of M in the Far Field.

Far-Field Perturbation.

Nondimensional Equations and Numerics.

Conclusions

Acknowledgment

Nomenclature

Appendix A: Proof of the Path Independence of M

Appendix B: Flow Potential for a Reiner–Rivlin Fluid

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

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

Introduction

Ice Rheology.

Analysis

Nye Solution.

Antiplane Shear.

Evaluation of the M Integral at the Channel.

Evaluation of M in the Far Field.

Far-Field Perturbation.

Nondimensional Equations and Numerics.

Conclusions

Acknowledgment

Nomenclature

Appendix A: Proof of the Path Independence of M

Appendix B: Flow Potential for a Reiner–Rivlin Fluid

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only