A note on the snout

Abstract

The shallow ice approximation for glaciers and ice sheets is a degenerate model in which the ice surface slope at the margin may be infinite. This result is due to the neglect of the otherwise small longitudinal stress terms. Here we derive a corrected approximation for the basal shear stress, and show that the resulting model provides an explanation for the observed finite slope margins.

1. Introduction

It is well known that the shallow ice approximation for glacier and ice-sheet flow, based on the assumption that $\delta \ll 1$, where

(1)$$\delta = {d\over l},\; $$

and $d$ is a typical depth scale and $l$ is a suitable length scale, leads to solutions which are singular at the margin, where the depth $h\to 0$. The singularity leads to an infinite slope. The same occurs if the ice is assumed to be perfectly plastic (Nye, 1967), when the ice thickness has a square root singularity at the margin. Nye nicely enunciated the problem as follows:

The theoretical treatments of glacier flow are all based on the approximation, which is valid far from the end, that the top and bottom surfaces of the glacier are almost parallel. If a model based on this approximation is extended towards the end, the top surface becomes more and more steeply inclined to the bed, and the approximation becomes useless before the end is reached.

Lliboutry (1956) addressed this issue, in particular concerning the observed finite slope at the snout of a glacier, but his methods were criticised by Nye (1957). Later, Nye (1967) simply stated that they were wrong (‘unsuccessful’), but his earlier comments were more muted in tone. Lliboutry (1958) responded to the criticism, and again Nye (1958) provides an answering (slightly exasperated) comment. The essence of the disagreement lies in the contrast between Nye's precision and Lliboutry's more heuristic approach.

Nye (1967) provided a theory to address the singularity at the snout based on a plastic flow law. Because his paper is based on a perfectly plastic flow law, the terminal snout angle is always $45^{\circ }$, which is much larger than normal values $\sim 15^{\circ }$ (e. g., Paterson, 1964). Nye adduces various reasons why his theory may nevertheless be reasonable, but the simplest reason may be that a Glen's law fluid is a better model to employ. Almost 50 years later, Nye (2015) returned to the subject, and provided an explanation for finite slope angle which indeed uses Glen's flow law. He states that Glen's (1961) observations of finite slope and compressive stress at the snout of Austerdalsbreen ‘are fully explained by a model based on the nonlinear ($n\approx 3$) Glen flow law that superposes longitudinal strain rate and simple shearing’.

Given this statement, one might wonder what the purpose of the present note is, since I have the same object in view. In elucidating Nye's note, one needs to understand what he means by laminar flow. This refers not to the smooth flow distinct from turbulent flow as a fluid dynamicist would suppose, but to an axial flow in which the velocity is parallel to the bed (Nye, 1952). Nye states that a suitable model is to allow a downstream shear flow with a compressive term (the compression causes an upwards velocity component, and is thus not laminar in Nye's sense). This observation is consistent with the work presented here, but rather than assuming a linear surface slope and compressive strain rate, I aim to deduce these from the governing equations of the motion. The present exposition is based on chapter 6 of my thesis (Fowler, 1977), which was motivated by the altercations of Nye and Lliboutry in the 1950s mentioned above.

2. Glacier flow model

The basic result of the analysis below is that near the snout the classical approximation for the basal shear stress must be amended by inclusion of a longitudinal stress term. This in itself is not new, having been established by Robin (1967), although in a different context (see also Budd and Radok, 1971). However, here we put this result in a deductive framework.

Dimensionless forms of the governing equations are well established (e. g., Fowler and Larson, 1978, equation (3.16)), so we will cut straight to these. We consider an isothermal slow two-dimensional flow of a glacier with axes $x$ along the mean bedrock slope, and $z$ upwards and normal to the $x$ axis. The momentum conservation equations take the form

(2)$$\eqalign{\tau_{3z}& = -1 + \mu s_x + \delta^{2}( p_x-\tau_{1x}) ,\; \cr p_z& = \tau_{3x}-\tau_{1z},\; }$$

subject to surface boundary conditions (of no stress)

(3)$$\eqalign{& \tau_3 = -\delta^{2}( p-\tau_1) s_x,\; \cr & p + \tau_1 + \tau_3s_x = 0\hbox{ on } z = s,\; }$$

and basal conditions of prescribed sliding velocity $u_b$ and zero normal velocity (melt-induced velocity being very small):

(4)$$w = ub_x,\; \quad {u + \delta^{2}wb_x\over \left(1 + \delta^{2}b_x^{2}\right)^{1/2}} = u_b\hbox{ on } z = b( x) .$$

Here, $\tau _1$ and $\tau _3$ are (dimensionless) tensile and shear components of the deviatoric stress tensor, which are related to the velocity components $u$ and $w$ by the flow law

(5)$$\tau_3 = \eta( u_z + \delta^{2}w_x) ,\; \quad \tau_1 = 2\eta u_x,\; $$

in which $\eta$ is the scaled viscosity; assuming Glen's law, this takes the form

(6)$$\eta = {1\over \tau^{n-1}},\; \quad \tau = \left(\tau_3^{2} + \delta^{2}\tau_1^{2}\right)^{1/2}.$$

Further, $p$ is the (dimensionless) pressure deviation from cryostatic overburden, $z = s$ is the ice upper surface, $z = b$ is the glacier bed, and subscripts $x$ and $z$ denote partial derivatives. The dimensionless parameter $\delta$ was defined in (1), and measures the shallowness of the flow: a typical value is $\delta \sim 10^{-2}$. The parameter $\mu$ is defined by

(7)$$\mu = {\delta\over S},\; $$

where $S$ is the mean bed slope, with a typical value of $\sim 0.1$. A typical value of $\mu$ may thus also be small (e. g., $\sim 0.1$), but not as small as $\delta$. It should be emphasised that the present discussion is restricted to valley glaciers with non-zero bed slope. Ice sheets with zero mean bed slope can be treated similarly (e. g., Fowler, 2011, pp. 631 ff.); essentially the downslope term $-1$ in (2) $_1$ is not present, and the (slightly different) non-dimensionalisation of the model leads to the replacement of $\mu$ by $1$.

The model is completed by the conservation of mass equation, which together with appropriate boundary conditions leads to the vertically integrated form

(8)$${{\islant20%\partial} {h}\over {\islant20%\partial} {t}} + {{\islant20%\partial} {}\over {\islant20%\partial} {x}}\left[hu_b + \int_b^{s}u\, dz\right] = a,\; $$

in which $h = s-b$ is the depth and $a$ is the accumulation rate (which is thus negative near the snout).

2.1 The shallow ice approximation

For simplicity we will take the bed to be flat, $b = 0$ (but still inclined to the horizontal). The shallow ice approximation simply ignores terms of $O( \delta )$, and the resulting integration leads to (8) in the form of a single equation for $h$,

(9)$${{\islant20%\partial} {h}\over {\islant20%\partial} {t}} + {{\islant20%\partial} {}\over {\islant20%\partial} {x}}\left[{( 1-\mu h_x) ^{n}h^{n + 2}\over n + 2} + hu_b\right] = a.$$

The sliding velocity needs to be prescribed, and we will assume a Weertman law of the form

(10)$$u_b = C\tau_b^{m},\; \quad m\leq n,\; $$

where $\tau _b$ is the basal shear stress, here equal to $\tau _3$ since the bed is flat. It is more realistic to suppose a dependence also on the effective pressure $N$, but this would simply distract from our focus. It is important to include sliding, otherwise we would run into the awkwardness of the contact line problem (Dussan V. and Davis, 1974), which refers to the fact that for an advancing fluid front (the snout) with a finite slope and no slip at the base, the local basal shear stress becomes infinite (and worse, non-integrable, implying an infinite force at the snout).

2.2 Marginal behaviour

In the shallow ice approximation, the basal shear stress is

(11)$$\tau_b = h( 1-\mu h_x) \, \colon \, $$

this is the dimensionless equivalent of the famous formula $\tau = \rho g h \sin \alpha$. In particular we assume that the ice slopes downhill, so that $1-\mu h_x> 0$, and this has been assumed in writing (9) and (12). We shall in fact assume that sliding is dominant, and neglect shearing in the ice altogether. The reason for this is that Fowler (1977) (pp. 116 ff., leading to equation (6.19)) showed that near the snout, sliding dominates shearing, so that a consideration of sliding only will suffice. In that case, we can take the velocity scaling to be such that $C = 1$ in (10). In the shallow ice approximation, (9) is then just

(12)$${{\islant20%\partial} {h}\over {\islant20%\partial} {t}} + {{\islant20%\partial} {}\over {\islant20%\partial} {x}}[ ( 1-\mu h_x) ^{m}h^{m + 1}] = a.$$

2.2.1 Steady state

It is worth pointing out that Nye (1967, 2015) assumed a steady state, even though the glacier under consideration (Austerdalsbreen) was in retreat: for this, see the following section. If the head of the glacier (where the ice flux is zero) is at $x = 0$, and we define the balance function

(13)$$B( x) = \int_0^{x}a( x') \, dx',\; $$

then the steady state of (12) satisfies

(14)$$( 1-\mu h_x) ^{m}h^{m + 1} = B( x) .$$

We can assume that the accumulation function $a( x)$ is a monotonically decreasing function of $x$, so that $B( x)$ first increases from zero at $x = 0$ to a maximum, and then decreases, passing through zero at a unique positive value $x_s$. Because the ice flux must be zero at the snout, this defines the snout position $x_s$. We can define the accumulation rate at the snout to be $a( x_s) = -1$ (negative because it represents ablation), by choice of the scale for $a$; then near the snout,

(15)$$( -\mu h_x) ^{m}h^{m + 1}\approx x_s-x$$

(because the slope is large and thus $-\mu h_x\gg 1$), and integrating this with the boundary condition that $h = 0$ at $x = x_s$, we find

(16)$$h\approx \left[{2m + 1\over ( m + 1) \mu}\right]^{{m\over 2m + 1}}( x_s-x) ^{{m + 1\over 2m + 1}},\; $$

with an almost square root singularity (it is a square root in the plastic limit $m\to \infty$).

2.2.2 The head of the glacier

Although we are mostly concerned with the snout, it is worth commenting on the behaviour near the head, as there is a difficulty there. We write (14) in the form

(17)$$h_x = {1\over \mu}\left[1-\left\{{h^{\ast }( x) \over h}\right\}^{( m + 1) /m}\right],\; \quad h^{\ast }( x) = \{ B( x) \} ^{1/( m + 1) }.$$

First note that if $\mu \ll 1$, then the solution is apparently just $h\approx h^{\ast }( x)$. With $m> 0$ and finite accumulation rate at the head, so that $B\sim x$ there, $h$ would seem to have an infinite positive slope, and the glacier surface slope would imply upstream flow! This makes no sense. In fact, for any positive value of $\mu$, (17) implies that the actual downstream slope (relative to the true horizontal) which is $\propto\ 1-\mu h_x$ must remain non-negative.

Let us examine the solution of (17) more closely. To do this, we consider the direction of trajectories in the $( x,\; \, h)$ ‘phase plane’, as shown in Figure 1. These cannot intersect, and slope respectively upwards and downwards above and below the curve $\Gamma$ defined by $h = h^{\ast }( x)$. The steady solution is that unique trajectory¹ which reaches $x = x_s$ and has a maximum at $x = x_M$, say, where it crosses $\Gamma$, and $h_x> 0$ for $x< x_M$. As $x$ decreases in $x< x_M$, $h_x$ increases as $x\to 0$ and reaches the value $1/\mu$ at $x = 0$, corresponding to a true horizontal surface.

Fig. 1. The curve $\Gamma$ given by $h = h^{\ast }( x)$ divides the $( x,\; \, h)$ plane into two regions. Above the curve, solutions of (17) point upwards, and below it they point downwards, as indicated by the arrows. The solution which terminates at the snout, where $h^{\ast } = 0$, is shown in red, and cannot reach zero at $x = 0$.

There are then two choices. We suppose the mountain slope has a peak. If we assume ice accumulates all the way to the peak, then $x = 0$ represents the peak, and the head of the glacier is at $x = 0$, with non-zero depth: it is actually an ice cap.

More realistic may be to suppose that the bed slope becomes sufficiently large that snow accumulation can not thicken to form ice, as avalanches remove it to lower altitudes. In that case we might suppose $a> 0$ for $x\geq 0$, but $a = 0$ (and thus also $B = 0$) for $x< 0$ (and the mountain peak is somewhere in $x< 0$). In that case $h$ is positive at $x = 0$, and $1-\mu h_x = 0$ for $x< 0$, as indicated by the dashed line in Figure 1. The ice surface is physically horizontal in $x< 0$ until $h = 0$.² This is a primitive representation of a bergschrund. In reality, snow accumulation on the flat surface would cause ice accumulation, and the model should be adjusted to allow $a = a( x,\; \, s_x)$ ($s$ not $h$ since when $h = 0$ it is the bedslope which is important). The two situations are illustrated in Figure 2.

Fig. 2. Two possible interpretations of the solution indicated in Figure 1. On the top an ice cap; on the bottom, a pair of mountain glaciers. The dashed line on the bottom indicates the position where $x = 0$.

2.3 Moving margins

As mentioned earlier, Nye (2015) was concerned with the finite surface slope of a retreating glacier. It is well known that the signature of an advancing glacier is a steep slope, while the slope of a retreating glacier is lower. A famous example is in the picture by Austin Post in Figure 3 which shows three glaciers in the Yukon: the outer two are advancing, while the central one is retreating (Post and LaChapelle, 2000, p. 38).

Fig. 3. Three congruent glaciers in the St. Elias Mountains, Yukon, Canada. Figure 38 of Post and LaChapelle (2000), reproduced with permission of the University of Washington Press, Seattle.

We can generalise the description of the local behaviour of the ice surface near a moving margin $x_s( t)$ by supposing $h\sim c( x_s-x) ^{\lambda }$, and then finding the leading order term balances in (12). Putting this into (12) gives (taking the accumulation rate $a( x_s) = -1$ as before)

(18)$$\eqalign{\lambda c( x_s-x) ^{\lambda -1}\dot{x}_s + \left[\left\{1 + \nu\lambda c( x_s-x) ^{\lambda -1}\right\}^{m}c^{m + 1}( x_s-x) ^{\lambda( m + 1) }\right]_x\approx -1.}$$

The relative sizes of the terms depend on the value of $\lambda$. If $\lambda \geq 1$, then the second term on the left is of $O[ ( x_s-x) ^{\lambda ( m + 1) -1}\ll 1$, and to balance the other two terms we must choose $\lambda = 1$, and then

(19)$$h\sim c( x_s-x) ,\; \quad \dot{x}_s = -{1\over c},\; $$

which describes retreat since $\dot {x}_s< 0$.

On the other hand, if $\lambda < 1$, then $-\mu h_x\gg 1$, the ablation term is much smaller than $h_t$, and so equating the two terms on the left-hand side of (18) leads to

(20)$$\lambda c( x_s-x) ^{\lambda-1}\dot{x}_s\approx -\left\{\mu^{m}\lambda^{m} c^{2m + 1}( x_s-x) ^{( \lambda -1) m + \lambda( m + 1) }\right\}_x.$$

Equating these requires $\lambda = {1\over 2}$ (thus indeed $< 1$), and then

(21)$$h\sim c( x_s-x) ^{1/2},\; \quad \dot{x}_s = \left({\mu c^{2}\over 2}\right)^{m},\; $$

which describes advance since $\dot {x}_s> 0$. Thus the slope is automatically finite during retreat, but becomes infinite during advance. This is associated with the fact that the diffusionless ($\mu = 0$) model forms shocks during advance, but retreats via a rarefaction wave (Fowler and Larson, 1980). However, this implies $u\to 0$ at the snout during retreat, in contradiction to Glen's (1961) observation. So in either case, a further examination is necessary.

3. Longitudinal stresses

Fowler (1977) used a convoluted sequence of scaling arguments to derive an approximate correction to the expression (11) for the basal shear stress, but we can abbreviate these by simply rescaling the equations (2)–(6) in the following way: we define

(22)$$\eqalign{x-x_s& = \varepsilon X,\; \quad h = \varepsilon^{1/( m + 1) }H,\; \ z = \varepsilon^{1/( m + 1) }Z,\; \cr \tau_3& = \varepsilon^{1/( m + 1) }T_3,\; \quad u = \varepsilon^{m/( m + 1) }U,\; \ \tau = {\delta^{1/n}\over \varepsilon^{1/( m + 1) n}}T,\; \cr \tau_1& = {1\over \delta^{( n-1) /n}\varepsilon^{1/( m + 1) n}}T_1,\; \quad p = {1\over \delta^{( n-1) /n}\varepsilon^{1/( m + 1) n}}P,\; }$$

where we choose

(23)$$\varepsilon = \delta^{{( m + 1) ( n + 1) \over ( m + 1) n + 1}},\; $$

thus $\varepsilon < \delta$. Suppose for example that $m = 2$, $n = 3$. Then $\varepsilon = \delta ^{6/5}$. A dimensionless length scale of $x_s-x\sim \delta$ corresponds, in view of the definition of $\delta$ in (1), to distances from the snout comparable to the depth, thus of the order of hundreds of metres. The extra reduction of $\delta ^{0.2}$ corresponds, if $\delta = 10^{-2}$, to a further reduction by a factor $0.4$. So the rescaling in (22) puts us in the last few hundred metres or so from the snout.

With these definitions, the equations become

(24)$$\eqalign{T_{3Z}& = -1 + MH_X + P_X-T_{1X},\; \cr P_Z + T_{1Z}& = \delta^{2\gamma}T_{3X},\; \cr U_Z& = \delta^{2\gamma}[ T^{n-1}T_3-w_X] ,\; \cr T_1& = {2U_X\over T^{n-1}},\; \cr T& = [ T_1^{2} + \delta^{2\gamma}T_3^{2}] ^{1/2},\; }$$

with

(25)$$\left\{\matrix{ T_3 = -( P-T_1) H_X,\; \cr P + T_1 + \delta^{2\gamma}T_3H_X = 0 }\right.\hbox{ on } Z = H,\; $$

where

(26)$$M = {\mu\over \delta^{1-\gamma}},\; \quad \gamma = {n-m + 1\over mn + n + 1},\; $$

and with an obvious notation, the sliding law (10) (with $C = 1$) is

(27)$$U_b = T_b^{m}.$$

Note that $0< \gamma < 1$. A typical value, with $m = 2$ and $n = 3$, is $\gamma = {1\over 5}$.

There are two awkwardnesses about this scaling. One is that the jump from $\tau \approx \tau _3$ to $T\approx {T_1}$ means that there must be another distinguished limit where both stresses are of equal size, and full Stokes flow may be appropriate; this would be expected to be where $x_s-x\sim \delta$, but in fact this complication is probably over-ridden by the fact that shearing is in any case negligible. The other awkwardness is the extra parameter $M$, which is formally large if $\mu = O( 1)$. This seems not to be an issue in practice, and we can for example formally assume that $M = O( 1)$ (with $m = 2$, $n = 3$, $\mu = 0.1$, $\delta = 0.01$, we have $M = 4$).

Ignoring small terms in (24) and (25), we find $P\approx -T_1$, $U\approx U( X)$, $T\approx \vert T_1\vert$, $T_1\approx ( 2\vert U_X\vert ) ^{1/n}\hbox {sgn}\, U_X$, and then on integrating,

(28)$$T_3 = ( 1-MH_X + 2T_{1X}) ( H-Z) + 2T_1H_X,\; $$

whence the basal shear stress is

(29)$$T_b = H( 1-MH_X) + 2[ H( 2\vert U_X\vert ) ^{1/n}\hbox{sgn} \,U_X] _X.$$

We now rewrite this in terms of the original macroscopic scales; this gives the correction to (11) as

(30)$$\tau_b = h( 1-\mu h_x) + \nu\{ h\vert u_x\vert ^{1/n}\hbox{sgn} \,u_x\} _x,\; $$

where

(31)$$\nu = ( 2\delta) ^{( n + 1) /n}.$$

For $\delta = 0.01$ and $n = 3$, we find $\nu = 0.005$. The point of this is that (30) gives a uniformly valid approximation to the basal shear stress, even away from the snout. For the predominantly sliding flow we assume, $h$ and $u$ are determined by

(32)$$\eqalign{h_t + ( hu) _x& = a,\; \cr u& = [ h( 1-\mu h_x) + \nu\{ h\vert u_x\vert ^{1/n}\hbox{sgn} \,u_x\} _x] ^{m}.}$$

Solutions of this problem are considered in the following section.

4. Local snout solutions

Fowler (1977) analysed steady solutions of (32) by using the method of strained coordinates (Van Dyke, 1975), which provides a uniformly valid approximate solution, that is to say a solution which provides an accurate approximation both in the main trunk of the glacier and also near the snout. The idea is that the singularity at the snout is displaced by straining the coordinate to a position which is not attained by the glacier. In general, one writes

(33)$$\eqalign{h& = h_0( \xi,\; t) + \nu h_1( \xi,\; t) + \ldots,\; \cr u& = u_0( \xi,\; t) + \nu u_1( \xi,\; t) + \ldots,\; \cr x& = \xi + \nu x_1( \xi,\; t) + \ldots,\; }$$

and chooses $x_1$ so that the higher approximations are no more singular than the first. For example, in the case of advance, (21) implies a square root singularity, $h_0\sim ( x_s-\xi ) ^{1/2}$. The choice of $x_1$ must be such that any singularity of $h_1$ must be no worse than this, i. e., $h_1/h_0$ remains bounded as $\xi \to x_s$. The leading order solution has a singularity at $\xi = x_s$, but the straining of $x$ moves this beyond the physical snout position. The procedure is, however, algebraically laborious, and beyond the scope of this note, so I will limit myself to showing that the model (32) does indeed provide for regular solutions at the snout.

We search for local expansions of the form

(34)$$h\approx c( x_s-x) ,\; \quad u\approx u_s + k( x_s-x) ,\; $$

and we allow $x_s$ to be time-dependent. Substitution of this into (32) shows that $k> 0$, as found by Glen (1961), and then

(35)$$u_s = \dot{x}_s + {1\over c} = ( \nu ck^{1/n}) ^{m},\; $$

which give the snout velocity and compressive strain in terms of the snout speed and the snout ice velocity, which can only be determined from global considerations. Smaller values of the snout slope and velocity are associated with retreat, with $\dot {x}_s\approx -{1\over c}$, as before in (19). For advance, $c$ is large, $u_s\approx \dot {x}_s$, as is consistent with shock propagation, and $ck^{1/n}\sim {1\over \nu }$. If for example we take $k = 1$, then this implies the actual slope (in terms of dimensional lengths) corresponds to an angle of $\tan ^{-1}( {1\over ( 2^{n + 1}\delta ) ^{1/n}})$, which for $\delta = 0.01$ and $n = 3$ is $28^{\circ }$.

5. Discussion

The main result of this note has been to show that the longitudinal stress correction to the basal stress introduced by Robin (1967) applies also to the snout of a glacier, and that it provides a perturbation which allows an explanation of both the finite slope and compressive stress at the snout. The treatment differs from that of Nye (2015) both in not assuming a steady state, and in providing a formulation which can be applied to the whole glacier (with some further consideration necessary at the glacier head).

Further insight would require a global solution of the boundary value problem (32). For the steady-state problem, this has been done by Fowler (1977) using the method of strained coordinates. Of more interest is the situation where the snout is advancing or retreating. In retreat, the unperturbed model can be applied all the way to the snout, and the perturbing compressive stress term provides a regular perturbation.

More challenging is the case of an advancing snout, when the perturbation is singular, and an aim of future work should be to combine the local snout description with the propagation of small amplitude surface waves to provide a uniform solution, in combination with linear wave theory (Nye, 1960; Fowler and Larson, 1980).