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 (Reference Robin1967), although in a different context (see also Budd and Radok, Reference Budd and Radok1971). However, here we put this result in a deductive framework.

Dimensionless forms of the governing equations are well established (e. g., Fowler and Larson, Reference Fowler and Larson1978, 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, Reference Fowler2011, 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, Reference Dussan and Davis1974), 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 (Reference Fowler1977) (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 (Reference Nye1967, Reference Nye2015) 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 trajectoryFootnote ¹ 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$.Footnote ² 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 (Reference Nye2015) 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, Reference Post and LaChapelle2000, p. 38).

Fig. 3.

Three congruent glaciers in the St. Elias Mountains, Yukon, Canada. Figure 38 of Post and LaChapelle (Reference Post and LaChapelle2000), 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, Reference Fowler and Larson1980). However, this implies $u\to 0$ at the snout during retreat, in contradiction to Glen's (Reference Glen1961) observation. So in either case, a further examination is necessary.

3. Longitudinal stresses

Fowler (Reference Fowler1977) 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 (Reference Fowler1977) analysed steady solutions of (32) by using the method of strained coordinates (Van Dyke, Reference Van Dyke1975), 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 (Reference Glen1961), 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 }$.