Asymptotic analysis of flow near a glacier terminus

Abstract

The non-linearly viscous ice flow in the vicinity of a glacier terminus, an observation region, depends crucially on the upstream flow as well as on the local surface and bed conditions. The former requires a likely complex solution of the balance laws and boundary conditions for the complete glacier. However, if the profile and downstream surface tangential velocity in the observation region are measured at an observation time t = 0, and a two-dimensional flow approximation is satisfactory, the complete stress and velocity fields satisfying local reduced model equations in the observation region at time t = 0 can be determined by asymptotic expansions in upstream distance from the (moving) terminus. Thus the full strain-rate and stress tensors are determined without prescribing the basal conditions. The terminus velocity is determined in terms of the net accumulation or melt flux and surface velocity at the terminus, with bounds for advance or retreat. The analysis and illustration are presented for a plane flow approximation.

1. INTRODUCTION

Nye (2015) constructed an approximate plane flow velocity distribution near the glacier terminus consistent with Glen (1961)'s measured profile and surface velocity distributions on Austerdalsbreen in 1958 and 1959. He adopts the conventional viscous power law for the ice, assumed to be close to melting so the temperature-dependent rate factor is constant (unity). It is then assumed that the driving shear stress is that given by the reduced model, an expression valid for small surface and bed undulation slopes relative to a mean bed plane which is close to horizontal. The next key simplification is that the longitudinal strain-rate is constant, independent of depth, defined by the observed, approximately constant, surface velocity gradient. This is incorporated in the viscous relation to construct, independent of any further prescribed surface and bed conditions, a non-laminar flow field from which a realistic compressive longitudinal stress magnitude is obtained.

An immediate question is whether it is necessary to make this longitudinal strain-rate simplification, or any alternative simplification. It is not. A plane flow local reduced model solution is now constructed for the terminus vicinity – an observation region – assuming a known surface profile and known surface tangential velocity at the observation time t = 0; a snapshot of the flow. The local reduced model is a set of leading order relations in a small parameter ε, which neglects terms of order ε compared with unity. ε is a co-ordinate and velocity scaling factor which is determined by the small dimensionless viscosity ν incorporating the observation region thickness magnitude d ₀.The validity of the local reduced model requires that the slopes of the bed undulations relative to the mean bed plane are not greater than ε in magnitude. There is no reference to conditions upstream of the observation region. Distinct balances arise for mean bed planes close to horizontal, slopes of order ε as in the above application and for steep glaciers with bed plane slope of order unity, but only the former relevant to Nye (2015)'s application is presented. The solution is achieved by Taylor series expansions in an upstream dimensionless distance ξ from the (moving) terminus, which covers a few thickness magnitudes. The Taylor series expansions in ξ are truncated when the next term does not exceed ε in magnitude. A bed tangential condition, such as a sliding law, is not required in this solution, but the solution incorporates the tangential traction and velocity at the bed and so imposes restrictions on any such relation. Further, the asymptotic expansion of the unsteady profile equation then determines the terminus velocity and instantaneous rate of change of the profile in terms of the net ablation/accumulation at the terminus and provides the criterion for the retreat/advance change-over. Analogous analyses can be performed for other two-dimensional flows, for example those describing a terminus spreading with radial or elliptic symmetry.

This observation region flow analysis is therefore of wider generality than Nye (2015)'s flow deduction for Austerdalsbreen, but comparison with the latter is used for illustration. In particular, the validity of the simplification and its consequences are noted.

2. TERMINUS VICINITY FLOW

Morland and Johnson (1980) developed the reduced model for steady plane flow of an ice sheet over a flat bed, then Morland and Johnson (1982) incorporated bed topography for near horizontal and steep mean bed planes. Subsequently, Morland (1984) presented the corresponding analysis for general three-dimensional unsteady flow over an inclined mean bed plane with arbitrary slope relative to the horizontal. The trivial solution of a horizontal ice surface – an ice reservoir – can be rejected by co-ordinate and velocity component scalings involving a small parameter ε related to a small dimensionless ice viscosity factor ν. Hutter (1981, 1983) introduced this co-ordinate and velocity stretching directly to reflect the large aspect ratio of an ice sheet, described as the shallow ice approximation (SIA). Here it is only a terminus vicinity flow which is being constructed independent of the likely complex upstream flow over varying bed slopes, but this can be described by a local reduced model derived by the same analysis with a length scale defined by the thickness magnitude of the observation region. The dimensionless reduced model relations derived by Morland (1984) will be adopted as a starting point. A plane flow approximation is adopted, as in Nye (2015), but analogous analyses follow for a terminus spreading with radial or elliptic symmetry, or any one-dimensional spreading.

The plane flow is described in rectangular Cartesian co-ordinates (x,z) with the x-axis pointing down a mean bed plane inclined at an angle χ to the horizontal, as shown in Figure 1. The associated velocity components are (u, w) with the velocity normal to the x, z plane v = 0. The origin represents the start of the glacier terminus vicinity – the observation region – over which the ice thickness reaches a magnitude d ₀; the Austerdalsbreen data d ₀ = 60 m is adopted for calculation. The moving terminus at time t is defined by x = x _M(t), with x _M(0) the observed location. The fixed bed is z = f(x), not necessarily plane, but with only undulations ofa small slope, and the unsteady glacier surface is z = h(x, t), also with small slope imposed by the reduced model equations. The slopes are exaggerated in the figure. Nye (2015) adopts an approximation tan χ = 0.095 based on his later bed profile measurements, and a constant surface slope γ = 0.182 relative to the bed at the observation time t = 0, shown in Figure 1. This corresponds to a constant surface slope with tangent 0.282 relative to the horizontal. Note that the surface slope 0.182 relative to the bed plane is not very small, so a reduced model which neglects terms of this magnitude compared with unity is only a rough approximation. Here a plane bed with small slope 0.095 is again adopted for illustration, but the ice surface is not assumed to be plane. The asymptotic analysis is performed retaining terms with magnitudes greater than or equal to ε so as not to compound the reduced model errors.

Fig. 1. Plane co-ordinate system (x, z) with x-axis down mean bed plane, showing ice surface and the horizontal from terminus location x _M at time t = 0. Slopes are exaggerated.

The ice is assumed to be incompressible with density ρ = 918 kgm⁻³, so the pressure p is a workless constraint not given by a constitutive law, but determined by momentum balance and boundary conditions. The deviatoric response is defined by a non-linearly viscous fluid relation between the strain-rate tensor D and the deviatoric stress tensor $ {\hat {\bi \sigma} } = {\bi \sigma } + p {\bi I} $, where σ is the Cauchy stress and I is the unit tensor. It is assumed that the ice in the terminus vicinity is temperate so the temperature dependent rate factor is unity. The conventional simplifications that D is co-axial with $ {\hat {\bi \sigma }}$, and the viscosity coefficient depends only on the second principal invariant J of $ {\hat {\bi \sigma }}$, are adopted. The viscous law in dimensionless form is

(1)$$ \displaystyle{{\bi D} \over {D_{0}}} = \psi({J}) {\displaystyle{ {\hat{ \bi \sigma}} \over {\sigma_{0}}}}, \quad {J} = {\displaystyle {{1} \over {2}}}{tr} \left(\displaystyle{ {\hat{\bi \sigma}} \over {\sigma_{0}}}\right)^{2}, $$

where the units σ ₀ and D ₀ are chosen so that the constitutive function ψ(J) is of order unity for deviatoric stresses and strain-rates arising in typical ice sheet flows; namely

(2)$$ \sigma_{0} = 10^{5}\,{Nm}^{-2}, \quad D_{0} = 1\,a^{-1}. $$

From Morland and Johnson (1980), the customary Glen's power law, (Glen, 1955, 1958), based on his laboratory experiments on ice close to melting, adopted by Nye (2015), is given by

(3)$$ \psi(J) = 1.5\,k(3J)^{({n-1})/{2}}, \quad k = 0.17, \quad n= 3.17, $$

which has an infinite viscosity at zero stress giving rise to singularities at the free surface. An alternative expression with smaller least squares error to the same data over a shear stress range 0–5 × 10⁵ Nm⁻², 0 ≤ J ≤ 25, is the polynomial (Smith and Morland, 1981) :

(4)$$ \eqalign{& \psi(J) = \psi_0 + \psi_1\,J + \psi_2\,J^{2},\cr & \psi_0 = 0.3336, \quad \psi_1 = 0.32, \quad \psi_2 = 0.02963,} $$

which has the mathematical advantage of a finite viscosity at zero stress. Solutions for both laws will be derived.

3. REDUCED MODEL

First introduce normalised dimensionless stress, strain-rate and time variables, denoted by an overbar:

(5)$$ (\,p, {\hat{\bi \sigma}}) = \rho g d_{0}( {\bar p}, {\hat{\bar{\bi \sigma}}}), \quad {\bi D} = q_{0} {\bar{\bi D}}/d_0, \quad t = d_0 {\bar t}/q_0, $$

where g =  9.81 ms⁻² and q ₀ = 1 ma⁻¹ is an accumulation/melt magnitude. Note that atmospheric pressure p _a = 1.013 × 10⁵ Nm⁻², assumed constant over the terminus vicinity, is not negligible compared with ice overburden pressure for these depths of ice, unlike the situation for large thick ice sheets; pointed out to me by Nye (personal communication). The viscous relation (1) now becomes

(6)$$ \eqalign{& {\bar{\bi D}} = \nu^{-1} \,\psi(J) {\hat{\bar{\bi \sigma}}}, \quad \nu = \displaystyle{ {{{\sigma}_0} q_0} \over {\rho g {d_{0}}^2 D_0}} = 0.003, \cr & J = \displaystyle{{1} \over {2}} \left( \displaystyle{{\rho g d_{0}} \over {\sigma}_{0}}\right)^{2} {tr}( {\hat{\bar{ \bi \sigma}}})^{2},} $$

where ν is a small dimensionless viscosity. Formal series expansions in ν results in an ice reservoir with horizontal surface – no descent to a terminus, and it follows that a non-horizontal surface requires co-ordinate and velocity scalings

(7)$$ \eqalign{& X = {\epsilon} {\bar x}/d_0, \quad Z = {\bar z}/d_0, \quad H = {\bar h}/d_0, \cr & U = {\epsilon} {\bar u}/q_0, \quad W = {\bar w}/q_0, \quad {\epsilon} \ll 1,} $$

where the capital variables and gradients are order unity. The viscous law (6) then becomes

(8)$$ \eqalign{&{\bar{\bi D}} = \nu^{-1}\,\psi( {\vartheta} {\check J}) {\hat{\bar{\bi \sigma}}}, \quad {\check J} = \displaystyle{{1} \over {2}}{\rm tr} \left(\displaystyle{ {\hat{\bar{\bi \sigma}}} \over { {\epsilon}}}\right)^{2}, \cr & {\vartheta} = \left( \displaystyle{\rho g d_{0}} \over {\sigma_{0}}\right)^{2}\, {\epsilon}^{2} = 29.2\, {\epsilon}^{2}.} $$

where for the power law given by (3) and the polynomial law given by (4),

(9)$$ \eqalign{&(3): {\epsilon} = \nu^{{n}/({n+1})} = 0.012, \quad {\vartheta} = 0.004; \cr & (4): {\epsilon} = \nu^{1/2} = 0.055, \quad {\vartheta} = 0.088.} $$

but any small-scale factor ε of similar magnitude provides the magnitude balance. It will be convenient to investigate both laws with a common scale factor ε = 0.04.

Now $ {\bar D}_{xx} = - {\bar D}_{zz} = O(1)$ so $ {\bar \sigma }_{xx} = - {\bar \sigma }_{zz} = O( {\epsilon }^{2})$, and assuming small mean bed inclination χ, of order ε or less, the lead order two-dimensional mass and momentum balances are

(10)$$ \eqalign{& {\displaystyle{\partial {U} \over {\partial {X}}}} + {\displaystyle{\partial {W}} \over {\partial {Z}}} = 0, \quad {\displaystyle{\partial { \bar \sigma}_{xz} \over {\partial {Z}}}} - {\epsilon}{\displaystyle{{\partial {\bar p}} \over {\partial {X}}}} + {\sin}\,\chi = 0, \cr & {\displaystyle{\partial { {\bar p}} \over {\partial {Z}}}} + {\cos}\,\chi = 0.} $$

These are subject to surface and bed conditions, treating the surface slope relative to the inclined bed plane, θ, as order ε,

(11)$$ \eqalign{& Z = H(X, {\bar t}): {\bar p} = {\bar p}_a = 0.187, \quad {\bar \sigma}_{xz} = 0, \cr & \displaystyle{{\partial {H}} \over {\partial {\bar t}}} + U\displaystyle{{\partial {H}} \over {\partial {X}}} - W = Q_n,} $$

(12)$$ Z = 0: \quad -W = B_n, \quad {\bar \sigma}_{xz} = F( {\bar p}_{b}, {\bar p}_a,U_b), $$

where Q _n and B _n are dimensionless fluxes, and $U_{b}(X,\, {\bar t})$ is the dimensionless velocity along the bed, all with unit q ₀. $ {\bar p}_b$ is the dimensionless total pressure on the bed where $ {\bar p}_{b} - {\bar p}_{a}$ is that due to the ice overburden. The actual relation, function F, will not be required, but the solution imposes restrictions on F. The lead order viscous relation becomes

(13)$$ {\displaystyle{\partial {U}} \over {\partial {Z}}} = {\epsilon}^{-1} {\bar \sigma}_{xz}\psi( {\vartheta} {\check J}), \quad {\check J} = ( {\epsilon}^{-1} {\bar \sigma}_{xz})^{2}. $$

With the restriction χ = O(ε), ε ⁻¹sin χ is order unity, while cos χ = 1, neglecting ε ², but is retained. Then integrating (10)₂ and (10)₃ in turn through the depth and applying the surface traction conditions (11) and the bed flux condition (12)₁ gives the lead order stress expressions

(14)$$ {\bar p} - {\bar p}_{\rm a} = {\cos}\,\chi[H(X, {\bar t}) - Z], \quad \Gamma(X, {\bar t}) = {\displaystyle{\partial {H}} \over {\partial {X}}}, $$

(15)$$ \eqalign{&{\bar \sigma}_{xz} = {\epsilon}\tau_{\rm s}(X, {\bar t})[H(X, {\bar t}) - Z], \cr & \tau_{\rm s}(X, {\bar t}) = {\epsilon}^{-1}{\sin}\,\chi - {\cos}\,\chi\, \Gamma(X, {\bar t}).} $$

Different expressions are obtained for a steep glacier bed χ = O(1). From (14) and (15) the tractions on the bed Z = 0 are

(16)$$ {\bar p}_{\rm b} = {\bar p}_{\rm a} + {\cos}\,\chi\,H(X, {\bar t}), \quad [ {\bar \sigma}_{xz}]_{\rm b} = {\epsilon}\tau_{\rm s}(X, {\bar t})\,H(X, {\bar t}), $$

and hence $H(X,\, {\bar t})$ determines the pressure on the bed, and $H(X,\, {\bar t})$ and $ \Gamma (X,\, {\bar t})$ determine the shear traction on the bed, independent of any prescribed sliding law.

Define

(17)$$ \tau(X,Z, {\bar t}) = {\epsilon}^{-1} {\bar \sigma}_{xz} = \tau_{\rm s}(X, {\bar t})[H(X, {\bar t}) - Z], $$

then $ {\check J} = \tau ^2$and the viscous relation (13) is

(18)$$ {\displaystyle{\partial {U}} \over {\partial {Z}}} = g(X,Z, {\bar t}) \; {\rm where} \; g(X,Z, {\bar t}) = \tau\psi( {\vartheta} \tau^{2}). $$

Integrating (18) with respect to Z from the bed Z = 0:

(19)$$ \eqalign{& U(X,Z, {\bar t}) = U_b(X, {\bar t}) + g_{1}(X,Z, {\bar t}), \cr & g_{1}(X,Z, {\bar t}) = \int_{0}^{Z} g(X,Z^{\prime}, {\bar t}) {d}Z^{\prime}.} $$

Incompressibility (10) is satisfied by introducing a stream function $ \omega (X,\,Z,\, {\bar t})$ such that

(20)$$ U = \displaystyle{{\partial { \omega}} \over {\partial {Z}}}, \quad W = -\displaystyle{{\partial { \omega}} \over {\partial {X}}}, $$

then integrating again with respect to Z:

(21)$$ \eqalign{&\omega(X,Z, {\bar t}) = \omega_{\rm b}(X, {\bar t}) + Z\,U_{\rm b}(X, {\bar t}) + g_{2}(X,Z, {\bar t}), \cr & \quad g_{2}(X,Z, {\bar t}) = \int_{0}^{Z} g_{1}(X,Z^{\prime}, {\bar t}) {d}Z^{\prime},} $$

where $ \omega _{\rm b}(X,\, {\bar t}) $ is the stream function value on the bed. Now the surface and bed kinematic conditions (11)₃ and (12)₁ are expressed by

(22)$$ \eqalign{& {\displaystyle{\partial {H}} \over {\partial { {\bar t}}}} + {\displaystyle{\partial} \over {\partial X}}\{[ \omega[X,H(X, {\bar t})]\} = Q_n, \cr & {\displaystyle{\partial} \over {\partial X}}\{ \omega[X,0, {\bar t})]\} = B_n,} $$

then differencing and eliminating ω by (21)₁ gives the profile equation

(23)$$ \eqalign{& {\displaystyle{\partial {H}} \over {\partial { {\bar t}}}} + {\displaystyle{{\partial} \over {\partial X}}}\{H(X, {\bar t})U_{\rm b}(X, {\bar t}) + g_{2}[X,H(X, {\bar t}), {\bar t}]\} \cr & = Q_n- B_n = -\bar{B}(X, {\bar t}),} $$

where $\bar {B}(X,\, {\bar t})$ is the net dimensionless melt flux.

4. TERMINUS ASYMPTOTICS

As near a divide and near a free surface, near the terminus the reduced model approximation is not valid – an inner boundary layer asymptotic expansion is required to be matched with the outer reduced model solution. However, in all these cases, the inner boundary layer expansion is passive. That is, the reduced model relations continue to satisfy all the boundary conditions through the inner boundary layer, so the outer reduced model solution is not influenced by matching with the boundary layer expansion. The inner region approximation of the full equations is, of course ,different from the reduced model solution, but here only in some small vicinity of the terminus not likely containing data points. Following is the reduced model analysis which will apply except in some small vicinity of the terminus.

First express variables as functions of the dimensionless time $ {\bar t}$ and a dimensionless distance ξ from the moving terminus $X = {X_{m}( {\bar t})}$, denoted by an over-check, and set X _m(0) = 0; thus

(24)$$ \eqalign{& \xi = X_{m}( {\bar t}) - X, \quad H(X, {\bar t}) = {\check{H}}(\xi, {\bar t}), \cr & \Gamma(X, {\bar t}) = -{\displaystyle{\partial { {\check{H}}}} \over {\partial {\xi}}} =- {\check{\Gamma}}(\xi, {\bar t}),} $$

(25)$$ \eqalign{&\displaystyle{{\partial {H}} \over {\partial { {\bar t}}}} = \displaystyle{{\partial { {\check H}}} \over {\partial { {\bar t}}}}+ {X_{m}^{\prime}( {\bar t})} {\check{\Gamma}}(\xi, {\bar t}), \cr & U(X,Z, {\bar t}) = {\check U}(\xi,Z, {\bar t}), \quad W(X,Z, {\bar t}) = {\check W}(\xi,Z, {\bar t}).} $$

Glen (1961)'s Austerdalsbreen measurements cover a distance 4 d ₀ from the terminus at time $ {\bar t} = 0$, equivalent to ξ = 4ε  ≈ 0.16 with a mean ε = 0.04 from (9). So consider a Taylor series expansion in small ξ for the profile :

(26)$$ \eqalign{& {\check H}(\xi, {\bar t}) = H_{1}( {\bar t})\xi + H_{2}( {\bar t})\xi ^{2} + H_{3}( {\bar t})\xi ^{3} + O(\xi^{4}), \cr & {\check{\Gamma}}(\xi, {\bar t}) = H_{1}( {\bar t}) + 2H_{2}( {\bar t})\xi + 3H_{3}( {\bar t})\xi ^{2} + O(\xi^{3}),} $$

where the coefficients at $ {\bar t} = 0$, H ₁(0), H ₂(0) and H ₃(0), are determined by correlation with measured profile data; Glen (1961)'s profile data are used in the later illustration for comparison with Nye (2015)'s solution. A profile expansion to order ξ ³, slope to order ξ ² and down-surface velocity expansion to order ξ ², results in a balance to order ξ ² of the profile equation (23). That is, an error of order ξ ³ compared with unity, with maximum value 0.004.

Now from (15), (17) and (18),

(27)$$ \eqalign{& \tau(X, {\bar t}) = {\check \tau} (\xi,Z, {\bar t})= {\check \tau}_{{\rm s}}(\xi, {\bar t}) [ {\check H}(\xi, {\bar t}) - Z], \cr & {\check \tau}_{{\rm s}}(\xi, {\bar t}) = \tau_{{\rm s}}(X, {\bar t}), \quad {\check{g}}(\xi,Z, {\bar t}) = {\check \tau}\psi( {\vartheta} {\check \tau}^{2}),} $$

and by (15) τ _s is order unity, so by (26) and (27), $ {\check \tau }$ is order ξ. For the power law,

(28)$$ \eqalign{&{\bar g}( X,Z, {\bar t}) = {\check g}(\chi,Z, {\bar t}) = a_{0} {\check \tau}^{n} \lt O(\xi^{3}), \cr & a_0 = 3k(3 {\vartheta})^{({n-1})/{2}} = 0.947,} $$

so $ {\check g}_{1}$ = O(ξ ⁿ⁺¹) and $ {\check g}_{2}$ = O(ξ ⁿ⁺²), which are neglected magnitudes. For the polynomial law,

(29)$$ \eqalign{& {\bar g}( X,Z, {\bar t}) = {\check g}(\chi,Z, {\bar t}) = a_{1} {\check \tau} + O(\xi^{2}), \cr & a_1 = 2\psi_{0} = 0.6672,} $$

from which

(30)$$ \eqalign{ {\check g}_{1}( \xi,Z, {\bar t}) & = {{\displaystyle{1} \over {2}}} a_{1} {\check \tau}_{{\rm s}}(\xi, {\bar t}) \{[ {\check H}(\xi, {\bar t})]^2 \cr & - [ {\check H}(\xi, {\bar t}) - Z]^2\} + O(\xi^3),} $$

(31)$$ \eqalign{{\check g}_{2}( \xi,Z, {\bar t}) & = \displaystyle{{{1} \over {2}}} a_{1} {\check \tau}_{{\rm s}}(\xi, {\bar t}) ([ {\check H}(\xi, {\bar t})]^{2}Z \cr & - \displaystyle{{1} \over {3}} \{[ {\check H}(\xi, {\bar t})]^{3}- [ {\check H} (\xi, {\bar t})-Z]^3\}) +O(\xi^4).} $$

Hence, by (19), for the power law $ {\check U}(\xi ,\,Z,\, {\bar t}) = {\check U}_{{\rm b}}(\xi ,\, {\bar t})$ neglecting terms of order ξ ⁿ⁺¹, which is independent of Z. This was deduced by Nye (2015) for the flow very close to the terminus, but the theoretical prediction shown in his Figure 4, reinforced by his comment in the text, implies a significant Z dependence further away, though still within the small ξ range corresponding to Glen (1961)'s measurements. That is, the actual flow solution is different from Nye (2015)'s approximate solution, independent of the measured data. For the polynomial law there is a Z dependence arising from $ {\check g}_{1}( \xi ,\,Z,\, {\bar t})$, but this is small, of order ξ ². Thus, the longitudinal velocity for both laws, neglecting order ξ ³, is expressed by

(32)$$ \eqalign{&{\check U}(\xi,Z, {\bar t}) = {\check U}_{{\rm b}}(\xi, {\bar t}) + K(\xi, {\bar t}), \cr & K(\xi, {\bar t}) = {\displaystyle{{1} \over {2}}} a_{1} {\check \tau}_{{\rm s}}(\xi, {\bar t}) \left\{[ {\check H}(\xi, {\bar t})]^2 - [ {\check H}(\xi, {\bar t})- Z]^2\right\},} $$

where a ₁ is zero for the power law, and is given by (29) for the polynomial law.

The later illustration will use Glen (1961)'s data which includes the down-surface velocity (the variation of the velocity at the surface, parallel to the surface), v _s(x _G) at time $ {\bar t} = 0$, where x _G is the distance along the curved surface. This interpretation, instead of the horizontal distance stated in a Figure caption, was advised by both Nye and Glen (personal communications). In dimensionless form

(33)$$ v_{\rm s}(x_{\rm G})/q_{0} = {\check V}_{\rm s}(\xi) = {\check U}[\xi, {\check H}(\xi,0),0], $$

where $ {\check {V}}_{\rm s}(\xi )$ can be approximated by an expansion

(34)$$ {\check V}_{\rm s}(\xi) = V_{0} +V_{1}\xi + V_{2}\xi^{2}, $$

with the coefficients V ₀, V ₁ and V ₂ determined by correlation with surface velocity data. The analysis applies for any such surface velocity data $ {\check {V}}_{\rm s}(\xi )$.

From (32), neglecting terms of order ε ³,

(35)$$ {\check V}_{\rm s}(\xi) = {\check U}_{\rm b}(\xi,0) + {{\displaystyle{1} \over {2}}} a_{1} {\check \tau}_{\rm s}(\xi,0) [ {\check H}(\xi,0)]^2 , $$

with a ₁ set to zero for the power law. Then by (34) and (35) the basal velocity at time $ {\bar t} = 0$ is

(36)$$ {\check U}_{\rm b}(\xi,0) = V_{0} +V_{1}\xi + V_{2}\xi^{2} - {{\displaystyle{1} \over {2}}} a_{1} {\check \tau}_{\rm s}(\xi,0)[ {\check H}(\xi,0)]^2, $$

and in turn the longitudinal velocity at time $ {\bar t} = 0$ becomes

(37)$$\eqalign{ {\check U}(\xi,Z,0) = V_{0} +V_{1}\xi + V_{2}\xi^{2} - {{\displaystyle{1} \over {2}}} a_{1} {\check \tau}_{\rm s}(\xi,0)[ {\check H}(\xi,0) - Z]^2,} $$

both determined completely by the profile and down-surface velocity data, independent of any basal sliding law. The corresponding longitudinal strain-rate is

(38)$$ \eqalign{{\check D}_{x x}(\xi,Z,0) & = -V_{1} - 2V_{2}\xi + a_{1} {\check \tau}_{\rm s}(\xi,0)[ {\check H}(\xi,0) - Z] {\check{\Gamma}}(\xi) \cr & + {{\displaystyle{1} \over {2}}} a_{1} {\check \tau}_{\rm s}^{\prime}(\xi,0)[ {\check H}(\xi,0) - Z]^2,} $$

with the longitudinal deviatoric stress $ {\check \sigma }_{x x}(\xi ,\,Z,\,0)$ given by the unique inverse of the monotonic relation (6)₁. Note that the longitudinal strain-rate is independent of Z for the power law, for which a ₁ = 0.

Now consider a net melting flux at $ {\bar t} = 0$ which varies with distance from the terminus, with expansion for small ξ

(39)$$ {\check B}(X,0) = B_{0} + B_{1}\xi + B_{2}\xi^{2}, $$

where B ₀ > 0 implies melting at the terminus, and apply the expansions to the profile Eqn (23), equating the coefficients of ξ ⁰, ξ, ξ ² to obtain the three identities

(40)$$ H_{1}(0)[X_{M}^{\prime}(0) - V_{0} ] = - B_{0}, $$

(41)$$\eqalign{ H_{1}^{\prime}(0) + 2H_{2}(0)X_{M}^{\prime}(0) - 2[V_{0}H_{2}(0) + V_{1}H_{1}(0)] = - B_{1},} $$

(42)$$ \eqalign{& H_{2}^{\prime}(0) + 3H_{3}(0)X_{M}^{\prime}(0) - 3[V_{0}H_{3}(0) + V_{1}H_{2}(0) + V_{2}H_{1}(0) \cr & \quad - {{\displaystyle{1} \over {2}}} a_{1} \kappa\,{cos}\,\chi H_{1}^{3}(0)] = -B_{2}.} $$

Immediately, by (40),

(43)$$ X_{M}^{\prime}(0) = V_0 - B_{0}/H_{1}(0), $$

which relates the terminus velocity $ X_{M}^{\prime }(0)$ linearly to the melt flux B ₀ at the terminus, and in particular the terminus is retreating if B ₀ > V ₀H ₁(0), which is a positive limit. The terminus can therefore be advancing for some range of positive melt flux. Then (41) and (42) in turn determine $ H_{1}^{\prime }(0)$ in terms of B ₁ and $H_{2}^{\prime }(0) $ in terms of B ₂.

In summary, a known profile and down-surface velocity in the terminus vicinity determine the current pressure, the longitudinal velocity, different from Nye (2015)'s approximation and hence the longitudinal deviatoric stress. Given the net melt (or accumulation) in the terminus vicinity, the current terminus velocity and instantaneous rate of change of the profile are determined. In total, a snapshot of the current flow and its rate of change in the terminus vicinity. Note that a solution independent of upstream conditions is possible only because the ice is assumed to be a viscous fluid relating only instantaneous variables, with no memory of its past.

5. AUSTERDALSBREEN ILLUSTRATIONS

The required data for correlation with (26)₁ and (34) to determine H ₁, H ₂, H ₃, V ₀, V ₁, V ₂ is given in Glen (1975)'s Table 1. Column 1 labels the stakes where measurements were made, column 2 the distance s _G from stake 1 along the ice surface, with s = s _G + 5 m an estimate of the distance from the terminus, column 3 is the height (vertical) h _G of the surface above stake 1, which, with Nye (2015)'s terminus slope 0.282, implies a height above the horizontal from the terminus of h _v = h _G + 1.4, and column 4 is the down-surface velocity v _s. Neglecting O(ε ³), $s = (x_M - x)(1 + (1/2)\, {\epsilon }^{2} \gamma _{M}^{2})$ where − γ _M is the terminus slope. Adopting Nye (2015)'s approximate surface slope 0.182 relative to the bed gives the approximation s = 1.016(x _M − x), so the distance along the x-axis from the terminus is x _M − x = 0.984 s. Finally the z surface co-ordinate is h = h _v sec χ(x _M − x) tan χ = 1.0045 h _v − 0.0953 (x _M − x). These relations determine the dimensionless co-ordinates $\bar {s}$, $ {\bar {x}}_M - {\bar {x}}$ and $ {\bar {h}}_{v}$, and the dimensionless down-surface velocity $ {\bar {v}}_{\rm s}$, then the scaled variables

(44)$$\eqalign{H &= {\bar h}, \quad \xi = {\epsilon}( {\bar x}_M - {\bar x}), \quad V = {\epsilon}\, {\bar v}_{\rm s}, \quad {\rm with} \cr & {\epsilon} = 0.04.} $$

A common scale factor ε, normalising ξ and $ {\bar v}_{\rm s}$ for both viscous laws, has been chosen. Table 1 below shows the data expressed in terms of these new variables.

Table 1. Surface profile and velocity data

Least squares correlations with the 12 surface data points for linear H(ξ), the constant surface slope approximation adopted by Nye (2015), and for a more accurate cubic H(ξ) adopted in the present analysis, give respectively

(45)$$ \eqalign{& H_1 = 4.6101,\quad H_2 = 0,\quad H_3 = 0; \cr & H_1 = 3.2752,\quad H_2 = 24.4453,\quad H_3 = -103.3843;} $$

with corresponding mean errors 0.0189 and 0.0140. The data points and correlations are shown in Figure 2 in unscaled (dimensionless) variables, with slopes exaggerated by the axes scaling. The linear correlation giving a constant slope ε H ₁ = 0.184 is nearly identical to Nye (2015)'s 0.182 approximation. While both correlations have similar error magnitudes, the corresponding surface slopes shown in Figure 3 reveal that the non-constant surface slope for the cubic correlation varies significantly from the constant slope of the linear correlation since the cubic correlation is able to correlate with data points both under and above the constant slope correlation.

Fig. 2. Profile (circles) data at $ {\bar t} = 0$ with linear (solid line) and cubic (dash-dot line) correlations. Horizontal inward from terminus shown as dashed line. Slopes exaggerated by axes scaling.

Fig. 3. Linear (solid line) and cubic (dash-dot line) correlation slopes.

Fig. 4. Down-surface velocity data (circles) at $ {\bar t} = 0$ with linear (solid line) and quadratic (dash-dot line) correlations.

However, these correlations are for data assumed accurate, but realistically there are significant data errors. In particular, Nye (2015) notes that the bed here assumed of constant slope, increases in the upper half of the measurement region, and he has suggested that depth measurements become less accurate moving away from the terminus. The presented illustration assumes accurate data points.

The linear and quadratic correlations for the 11 down-surface velocity points give respectively

(46)$$ \eqalign {& V_0 =0.3242,\quad V_1 = 6.0427, \quad V_2 = 0; \cr & V_0 = 0.3712,\quad V_1 = 3.9398,\quad V_2 = 14.2934;} $$

with corresponding mean errors 0.0330 and 0.0205. Figure 4 shows the down-surface velocity $ {\bar v}_{\rm s}$ data points and correlations in unscaled dimensionless variables, and while both correlations have similar error magnitudes, the quadratic correlation matches closely the points both under and over the linear correlation. The velocity at the terminus is approximately 8.1 m a⁻¹. Figure 5 shows the corresponding longitudinal strain-rates $-\partial \,v_{\rm s}/\partial \, {\bar {x}}$. The constant strain-rate −∂v _s/∂x = 6.043/60 = 1.067 a⁻¹ with the linear correlation corresponds to Nye (2015)'s constant r ₀ approximation, though he estimates this to be 0.064 a⁻¹. In fact his estimate is very close to the terminus value of the quadratic approximation 3.94/60 = 0.066 a⁻¹. Again, the more accurate quadratic correlation shows a significant variation from the constant strain-rate correlation. From (40),

(47)$$ \eqalign{& X_{M}^{\prime}(0) = 0.3242 - 0.2169B_0, \cr & {\rm or} \; X_{M}^{\prime}(0) = 0.3712 - 0.3053B_0,} $$

for the constant slope and quadratic slope correlations respectively, giving retreat/advance flux limits B ₀ = 1.4622 and B ₀ = 1.2158 m a⁻¹.

Fig. 5. Linear( solid line) and quadratic (dash-dot line) correlation longitudinal velocity gradients $\partial \,\bar {v}_{\rm s}/\partial \,\bar {x}$.

6. CONCLUSIONS

An asymptotic analysis of the unsteady flow in a glacier terminus vicinity based on the reduced model for small surface and bed slopes has shown that surface profile and down-surface velocity data at an observation time determine the longitudinal velocity independent of basal sliding and melt or accumulation boundary conditions. The resulting longitudinal strain-rate and deviatoric stress determined at that time impose a restriction on any basal sliding law, independent of the melt condition. In addition, the terminus velocity and the current rate of change of the surface profile are determined in terms of the melt/accumulation flux. However, since different profile and velocity correlations with errors of similar magnitude show significant differences in slope and velocity gradient, it must follow that errors in data points, which would change the correlations, could also yield significant differences. The present analysis determines accurate solutions for the given data.

One conclusion, independent of the data, is that the longitudinal velocity given by (37) shows negligible depth variation for the viscous power law, in contrast to Nye (2015)'s prediction of significant depth variation, and only small depth variation for the viscous polynomial law.

Illustrations are presented with the linear and cubic correlations of surface profile data, and linear and quadratic correlations of the velocity data, for Austerdalsbreen (Glen, 1961).