Bed topography and marine ice-sheet stability

Abstract

This paper examines the effect of basal topography and strength on the grounding-line position, flux and stability of rapidly-sliding ice streams. It does so by supposing that the buoyancy of the ice stream is small, and of the same order as the longitudinal stress gradient. Making this scaling assumption makes the role of the basal gradient and accumulation rate explicit in the lowest order expression for the ice flux at the grounding line and also provides the transcendental equation for the grounding-line position. It also introduces into the stability condition terms in the basal curvature and accumulation-rate gradient. These expressions revert to well-established expressions in circumstances in which the thickness gradient is large at the grounding line, a result which is shown to be the consequence of the non-linearity of the flow. The behaviour of the grounding-line flux is illustrated for a range of bed topographies and strengths. We show that, when bed topography at a horizontal scale of several tens of ice thicknesses is present, the grounding-line flux and stability have more complex dependencies on bed gradient than that associated with the ‘marine ice-sheet instability hypothesis’, and that unstable grounding-line positions can occur on prograde beds as well as stable positions on retrograde beds.

1. Introduction

The discovery that sections of the West Antarctic Ice Sheet are losing mass to the oceans (Wingham and others, 1998; Rignot, 1998; Shepherd and others, 2018) with significant implications for global sea levels (e.g. Shepherd and others, 2018) has given practical importance to understanding the evolution of the grounding line – the boundary between the grounded and floating parts – of marine ice sheets. The general problem of the waxing and waning of a marine ice sheet is a time-variant problem involving varying external forcings, the ice sheet's adjustment to these and the interactions between them (e.g. Fyke and others, 2018). Nonetheless, one can gain some insight into the dynamics by considering the much simpler question of whether, in conditions of constant external forcing, a steady-state marine ice sheet is possible at all, and if so, whether it is stable subjected to small perturbations. Since the original framing of Weertman (1974), the matter has received continuous theoretical attention (e.g. Chugunov and Wilchinsky, 1996; Hindmarsh and Le Meur, 2001; Schoof, 2007a, 2007b, 2012; Nowicki and Wingham, 2008; Robel and others, 2014; Tsai and others, 2015; Robel and others, 2016; Kowal and others, 2016; Schoof and others, 2017; Haseloff and Sergienko, 2018; Pegler, 2018), principally arising from the contention of Weertman (1974) that a marine ice sheet cannot exist in steady state on a retrograde bed (i.e. one that deepens towards the ice-sheet interior), a hypothesis that has become known as the ‘marine ice-sheet instability hypothesis’.

The most complete account has been given by Schoof (2007a, 2011), who, by restricting attention to rapidly-sliding ice stream flow, reduced the problem to an ordinary differential equation. In a later paper, Schoof (2012) showed that resulting grounding line locations require prograde bed slopes to be stable. While recent work (Kowal and others, 2016; Schoof and others, 2017; Haseloff and Sergienko, 2018; Pegler, 2018; Sergienko and Wingham, 2019) has questioned the instability hypothesis, it remains a widely-accepted description of the higher driving stress ice streams such as those in streams that characterise, for example, the Amundsen Sea sector. Since the original proposition of Mercer (1968), melt-triggered instability has been frequently invoked to describe the future evolution of sections of the West Antarctic ice sheet, as, for example, do recent descriptions of the future of the Thwaites Glacier basin (e.g. Joughin and others, 2014; Feldmann and Levermann, 2015; Cornford and others, 2015; Seroussi and others, 2017).

The instability hypothesis arises from there being a monotonically increasing relationship between the thickness at the grounding line and the discharge through it. This, in turn, depends on the assumption that the thickness gradient in the momentum equation is large in comparison with the basal gradient. While this is a reasonable description of ice streams flowing over smooth beds offering considerable resistance to the flow, ice streams can experience a variety of frictional resistance at their beds and flow over bed topography that contains a wide range of spatial scales as Figure 1 illustrates. Since surface expressions of shorter scale topography (e.g. several tens of ice thicknesses wavelength) is generally greatly subdued in the surface gradient (e.g. Gudmundsson, 2003), it follows that the basal and thickness gradients must generally be of the same order, and in such a case there is no longer good reason to suppose a monotonic relationship between grounding line thickness and flux, nor a correspondingly simple relationship between basal gradient and stability.

Fig. 1. Characteristics of the Pine Island Glacier and Thwaites Ice Stream region. (a) Magnitude of the surface gradient $\vert \vec \nabla S\vert$, $S = b + h$; (b) magnitude of the driving stress (kPa) computed as $\tau _d = \rho g h\vert \vec \nabla S\vert$ derived from the BedMachine Antarctica data set (Morlighem and others, 2020); (c) surface $S$ (blue line) and bed $b$ (red line) elevation along a flowline on Thwaites Glacier (red line in panel (b)).

In this paper, we examine the role of the basal topography in the existence and stability of steady grounding line locations. To do so, we take a scaling approach that assumes that the buoyancy force is small and similar in magnitude to the longitudinal-stress gradient. By doing so, the effect of the basal gradient in the momentum equation then naturally occurs in the lowest order relation between ice flux and thickness at the grounding line, which also supplies the transcendental equation for the steady grounding-line locations, assuming these exist. We also find that, in cases where the thickness gradient at the grounding line is large, our result reduces algebraically to the monotonic dependence of flux on the thickness of Schoof (2007a). We also consider the stability of the steady-state grounding line locations, and find that this now depends on the basal curvature in addition to the basal gradient and the simple relationship between basal gradient no longer holds. The stability condition that informs the ‘marine ice-sheet instability hypothesis’ no longer generally applies.

2. Model description

The model is the same as that proposed by Schoof (2007a). Here we provide a brief description. The flow of an unconfined ice stream into an unconfined ice shelf (Fig. 2) can be described by a vertically integrated momentum balance under assumptions of a negligible vertical shear appropriate for ice-stream and ice-shelf flow (MacAyeal, 1989)

(1a)$$2\left(A^{-( {1/n}) }h\left\vert u_x\right\vert ^{( {1/n}) -1}u_x\right)_{x} - \tau_{b} -\rho g h\left(h + b\right)_x = 0 \quad \rm{\,for}\ {\it x}_{\rm d}\leq {\it x}\leq {\it x}_{\rm g}$$

and

(1b)$$ 2\left(A^{-( {1/n}) }h\left\vert u_x\right\vert ^{( 1/n) -1}u_x\right)_x \;-\;\rho g^\prime hh_x = 0 \quad \rm{\,for}\ {\it x}_{\rm g}\leq {\it x}\leq {\it x}_{\rm c}.$$

Here $u( x)$ is the depth-averaged ice velocity, $h( x)$ ice thickness, $b( x)$ is bed elevation (negative below sea level and positive above sea level), $A^{-1/n}$ is the ice stiffness parameter (assumed to be constant), $n$ is the exponent of Glen's flow law, $g$ is the acceleration due to gravity, $\tau _{\rm b}$ is basal shear. $g^\prime$ is the reduced gravity defined as

(2)$$g^\prime = \delta g$$

where

(3)$$\delta = {\rho_{\rm w}-\rho\over \rho_{\rm w}}$$

is the buoyancy parameter. $\rho$ and $\rho _{\rm {w}}$ are the densities of ice and water, respectively, $x_{\rm {d}}$ is the location of the ice divide (assumed to be $x_{\rm {d}} = 0$), $x_{\rm c}$ is the location of the calving front and $x_{\rm {g}}$ is the location of the grounding line. Basal sliding is described by a power-law function

(4)$$\tau_{\rm b} = C\left\vert u\right\vert ^{m-1}u,\; $$

where $C$ and $m$ are constant parameters.

Fig. 2. Model geometry: $b$ – bed elevation ($b\lt$0), $h$ – ice thickness, $x_{\rm {d}}$ – the ice divide location, $x_{\rm {g}}$ – the grounding line location; $x_{\rm c}$ – the calving front location.

The mass balance is

(5)$$\matrix{ h_{\rm t} + ( uh)_x = \left\{{\matrix{ \dot a \hfill\qquad \rm{for}\ 0\leq {\it x}\leq {\it x}_{\rm g},\hfill \cr \dot m \hfill \qquad \rm{for}\ {\it x}_{\rm g}< {\it x}\leq {\it x}_{\rm c}, \cr } } \right.} $$

where $\dot a$ and $\dot m$ are net accumulation/ablation (positive for accumulation) rates of the ice stream and ice shelf, respectively.

Boundary conditions at the divide $x_{\rm {d}}$ and the calving front $x_{\rm c}$ are

(6a)$$( h + b) _x = 0,\; \quad u = 0 \quad \rm{at}\ \quad {\it x} = {\it x}_{\rm d},\; $$

and

(6b)$$ 2 A^{-( {1/n}) }h\left\vert u_x\right\vert ^{( 1/n) -1}u_x = {1\over 2}\rho g^\prime h^2 \quad \rm{at}\ \quad {\it x} = {\it x}_{\rm{c}}.$$

At the grounding line $x_{\rm {g}}$ the continuity conditions

(7a)$$u_{\rm stream}( x_{\rm g}) = u_{\rm shelf}( x_{\rm g}) ,\; $$

(7b)$$h_{\rm stream}( x_{\rm g}) = h_{\rm shelf}( x_{\rm g}) $$

and

(7c)$$ \tau_{\rm stream}( x_{\rm g}) = \tau_{\rm shelf}( x_{\rm g}) ,\; $$

where $\tau = 2 A^{-( {1/n}) }h\left \vert u_x\right \vert ^{( 1/n) -1}u_x$ is the longitudinal stress, and the flotation condition

(8)$$h( x_{\rm g}) = -{\rho_{\rm w}\over \rho}b( x_{\rm g}) $$

are satisfied. The fact that the ice is grounded upstream of the grounding line and is floating downstream of it is reflected by two inequalities

(9a)$$h( x) \ge -{\rho_{\rm w}\over \rho}b( x) \quad\rm{\,for}\ {\it x}_{\rm d}< {\it x}< {\it x}_{\rm g},\; $$

and

(9b)$$ h( x) < -{\rho_{\rm w}\over \rho}b( x) \quad\rm{\,for}\ {\it x}_{\rm g}< {\it x}< {\it x}_{\rm c}.$$

In circumstances where ice shelves are unconfined, the momentum balance of the ice shelf (1b) can be integrated analytically with the boundary condition at the calving front, $x_{\rm c}$, (6b) and the continuity conditions (7), and the problem can be reduced to the ice-stream part only with the boundary conditions at the grounding line provided by the flotation condition (8) and the stress condition

(10)$$2 A^{-( {1/n}) }h\left\vert u_x\right\vert ^{( 1/n) -1}u_x = {1\over 2}\rho g^\prime h^2 \quad \rm{at}\ {\it x} = {\it x}_{\rm g}.$$

3. Grounding-line position, flux and stability

In this paper, we wish to study the effect of the basal topography and strength on the grounding line flux and stability. To gain some analytic insight requires us to approximate the model of section 2, which we achieve in this section by supposing the longitudinal stress gradient and the buoyancy to be small parameters. We introduce the characteristic scales for the problem: the bathymetry of a continental shelf $[ b]$, the length of an ice steam $[ x]$ and its velocity $[ u]$. With these scales, the scales of other parameters of the problem are

(11)$$[ h] = [ b] ,\; \quad [ \dot a] = {[ h] [ u] \over [ x] } \quad{\rm and}\ [ t] = {[ x] \over [ u] },\; $$

and three non-dimensional parameters then emerge. These are $\delta$, given by (3),

(12a)$$\varepsilon = {2A^{-( {1/n}) }[ u] ^{1/n}\over \rho g[ h] [ x] ^{1/n}}$$

and

(12b)$$ \gamma = {[ C] [ u] ^m[ x] \over \rho g[ h] ^2},\; $$

where $\varepsilon$ is characteristic of the ice-stream deformation, and $\gamma$ is characteristic of the ice-stream sliding.

The non-dimensional problem is

(13a)$$\varepsilon\left(H\left\vert U_X\right\vert ^{( 1/n) -1}U_X\right)_X - \gamma\left\vert U\right\vert ^{m-1}U-H( H + B) _X = 0 \quad{\rm for}\ 0\leq X\leq X_{\rm g},\; $$

(13b)$$H_T + \left(UH\right)_x = \dot a,\; $$

(13c)$$( H + B) _X = 0,\; \quad U = 0 \quad {\rm at}\ X = X_{\rm d},\; $$

(13d)$$\varepsilon H\left\vert U_X\right\vert ^{( 1/n) -1}U_X = {1\over 2}\delta H^2 \quad {\rm at}\ X = X_{\rm g}$$

and

(13e)$$ H = -{B\over 1-\delta} \quad {\rm at}\ X = X_{\rm g}.$$

The behaviour of the system (13) depends on the relative magnitudes of the three non-dimensional parameters $\delta$, $\gamma$ and $\varepsilon$. In this paper, we will concern ourselves with ice streams for which $\gamma$ is O(1), describing a flow in which the driving stress and the basal friction are of the same order. Indeed, we will set $\gamma$ equal to one. Equation (12b) can then be regarded as a constraint on the independence of the three scales; non-dimensionally it insists that the driving stress and basal friction are of the same order. In this and the next section, we will also limit ourselves to situations in which $\varepsilon$ is small, that is, in which the longitudinal stress gradient is small in comparison with the two other terms of the momentum equation.

There is then finally the choice of the parameter $\delta$. Numerically, its value is 0.1, and in his treatment of (13), Schoof (2007a, 2011) took $\delta$ to be large in comparison with $\varepsilon$. This has the consequence that at zeroth order in $\varepsilon$, the stress condition at the grounding line (13d) becomes $h$ = 0 at the grounding line, and the velocity $u$ becomes unbounded, leading Schoof (2007a) to introduce a boundary layer solution near the grounding line. As a consequence of the flotation condition, it then follows that the bed depth is also small at the grounding line. We discuss his approximate solution in greater detail in Appendix A. Here, however, we note that this solution is sensitive to the bed topography only through its value at the grounding line, and not, for example, to its gradient and curvature. The alternative, which we explore here, is to take $\delta$ to be a small parameter of the same order as $\varepsilon$, the approach taken by Sergienko and Wingham (2019) in analysing the corresponding problem when the bed is very weak.

When $\delta \sim \varepsilon$ the terms on either side of the stress condition (13d) are of the same order, and, as noted by Schoof (2007a), no boundary layer is needed to obtain a solution at lowest order.

3.1. Ice flux at the grounding line

It is convenient to use the flux $Q = UH$ instead of the ice velocity. Taking $\varepsilon \sim \delta \ll$1, (13a)–(13d) then become to lowest order

(14a)$$-\left\vert Q\right\vert ^{m-1}Q-H^{m + 1}( H + B) _X = 0,\; $$

(14b)$$H_T + Q_X = \dot a,\; $$

(14c)$$( H + B) _X = 0,\; \quad U = 0 \quad {\rm at}\ X = X_{\rm d}$$

and

(14d)$$ \eqalign{Q_{X}H^{m + 2} + Q^{m + 1} + QH^{m + 1}B_X = \left({\delta\over 2\varepsilon}\right)^nH^{n + m + 3} \quad {\rm at}\ X = X_{\rm g}.}$$

To obtain the stress condition (14d) we have used

(15)$$U_X = {Q_X\over H}-Q{H_X\over H^2}$$

in (13d). After rearranging terms, and using the momentum balance (14a) to express $H_X$

(16)$$H_{X} = -{1\over H}\left({Q\over H}\right)^m-B_X,\; $$

(13d) becomes (14d), which is the general form of the relationship between the ice flux and the ice thickness at the grounding line. For ice-sheet configurations evolving with time, according to (14b), $Q_X = \dot a - H_T$, and (14d) becomes

(17)$$( \dot a- H_T) H^{m + 2} + Q^{m + 1} + QH^{m + 1}B_X = \left({\delta\over 2\varepsilon}\right)^nH^{n + m + 3}\quad {\rm at}\ X = X_{\rm g}.$$

3.2. Steady-state grounding-line ice flux, position and stability

In a steady state, $H_T = 0$, and the relationship between ice flux and ice thickness at the grounding line (17) becomes

(18)$$Q^{m + 1} + QH^{m + 1}B_X = \left({\delta\over 2\varepsilon}\right)^nH^{n + m + 3}-\dot aH^{m + 2}.$$

This expression (18) is not new. The same equation can be found in Schoof (2007a, 2007b) in slightly different guises.

If, at the grounding line, the ice thickness gradient $H_X$ is large in comparison with the bed gradient $B_X$ in (16), and is also large enough to dominate the velocity gradient in (15), then $B_X$ and $\dot a$ may be dropped from (18) to obtain

(19)$$Q_{{\rm S}}^{m + 1} = \left({\delta\over 2\varepsilon}\right)^{n}H^{m + n + 3},\; $$

where $Q_{\rm S}$ is used to distinguish it from (18). This is the approximate solution obtained by Schoof (2007a, 2011) to his boundary layer problem when $\varepsilon \ll \delta \ll$ 1. (18) provides (19) as a special case when $\dot a$ and $B_X$ are small at the grounding line. That the two approaches provide the same algebraic solution is a consequence of Schoof's approximation, which is equivalent to supposing the longitudinal stress is negligible even within the boundary layer. The fact that this approximation is successful even when $\delta$ is O(1) is a consequence of the non-linearity of the ice flow and the chosen form of the sliding law. This is described in more detail in Appendix A; for our purposes here it is sufficient to note that (18) covers both cases with small and non-negligible basal slopes in the vicinity of the grounding line.

The location of the grounding line $X_{\rm {g}}$ is determined through the addition of a flotation condition to (14). Values of $X_{\rm {g}}$ are then determined as the roots of the transcendental equation (18) that satisfy the condition (of which there may be none, one or several). When $\delta \sim$ o(1), the flotation condition (13e), becomes $H = -B$ to lowest order. We have found, however, that in contrast to (18), this order of approximation does not provide accurate values of $X_{\rm {g}}$. Formally, this can be overcome by developing higher order approximations to (13). We have found that while adding considerable algebraic complexity to (18), the additional terms are very small numerically. It is only in determining $X_{\rm {g}}$ that the increased accuracy in $\delta$ is needed and we therefore use (18) together with (13e). Physically, this is not surprising. With gentle beds, the location of the steady-state grounding line is very sensitive to the buoyancy.

We now turn to examining the stability of the steady-state grounding-line locations. We use the same approach taken by Schoof (2012). This is conceptually straightforward, in that we subject the steady states to small, time-dependent perturbations about the steady-state solutions, which we now distinguish with the $\hat{}$ decoration, so that, for example, $X_{\rm {g}}( T) = \hat X_{{\rm g}} + \tilde X_\rm {g}( T)$. We substitute such expressions into (14), and by assuming that the time-variant perturbations are small, formulate a linearised perturbation problem. This results in solutions of the form

(20)$$\tilde X_{\rm g}( T) = \tilde X_{\rm g}( 0) {\rm e}^{\Lambda T}$$

for the grounding-line evolution. While the resulting problem is an eigenvalue problem for $\Lambda$, and we are not able to provide a closed form expression for $\Lambda$, we can as Schoof (2012) provide conditions which, if satisfied, ensure all eigenvalues are negative, and as a consequence the grounding-line position is a stable one. The derivation is algebraically complex, and we relegate it to Appendix B. It concludes that provided

(21a)$$\hat H_{X}< -{1\over 1-\delta}B_X\ {\rm at}\ X = {\hat X}_{{\rm g}}$$

and

(21b)$$ \hat H_{X}< -{m\over m + 1}B_X\ \ {\rm at}\ X = {\hat X}_{{\rm g}}$$

the grounding line is stable if, and only if,

(22)$$\eqalignno{& -{B_X\over 1-\delta}\Big\{\left({\delta\over 2\varepsilon}\right)^n( m + n + 3) \hat H^{n + m + 2}-\dot a( m + 2) \hat H^{m + 1}\cr & \qquad- ( m + 1) \hat Q\hat H^mB_X\Big\}\cr & > \Big\{\hat H^{m + 2}\dot a_{X} + \dot a\left[( m + 1) \hat Q^{m} + \hat H^{m + 1}B_X\right] + \hat Q\hat H^{m + 1}B_{XX}\Big\}\quad\cr & \qquad{\rm at}\ X = {\hat X}_{{\rm g}}.&}$$

If the conditions (21) are not met, it is not possible to determine the stability without solving the eigenvalue problem itself. Condition (21a) is met if the stream remains grounded upstream of the grounding line. It is satisfied for any steady-state solution of interest. With $m > 0$ and $\delta < 1$, if $B_X( \hat X_{{\rm g}}) \gt$ 0, and the bed is retrograde, then (21b), is satisfied if (21a) is met. On the other hand, if $B_X( \hat X_{{\rm g}}) \lt$ 0 and the bed slope is pro-grade, (21b) must be met independently of (21a) for (21) to be satisfied.

Under the same conditions of small $B_X$ and negligible $\dot a$, under which (18) reduces to (19), the only remaining terms in (22) are

$${{B_X\over 1-\delta}\left({\delta\over 2\varepsilon}\right)^n( m + n + 3) H^{n + m + 2}}$$

on the left-hand side, and $\displaystyle {\dot a( m + 1) Q^{m}}$ on the right-hand side. With (19), (22) reduces to

(23)$${m + n + 3\over m + 1}\left({\delta\over 2\varepsilon}\right)^{\textstyle{n\over m + 1}}\hat H^{\textstyle{m + n + 3\over m + 1}-1}\left(-{B_X\over 1-\delta}\right)> \dot a.$$

The left-hand side can be recognised as the derivative of $Q_{{\rm S}}$ of (19), with the flotation condition applied, and so (23) may be written

(24)$$Q^\prime_{\rm S}> \dot a$$

which is the stability condition determined by Schoof (2012). Schoof's result is also subject to the condition (21a), but (21b) is a new condition that arises from the additional terms in (18) over those in (19).

4. The effect of bed shape and strength on steady-state configurations

Through supposing $\varepsilon$ and $\delta$ to be of the same order, section 3 provides expressions for the grounding line flux and stability that acknowledge the role of the basal topography. These expressions are considerably more complex than those that result if the thickness gradient near the grounding line is sufficiently large that the contributions of the basal topography and the accumulation rate to the local strain rate $U_X$ are small. We might expect that the thickness gradient will increase as the strength of the bed, described here by the parameter $C$, increases, and therefore that the significance of the bed topography will become more marked as the bed weakens and the thickness gradient reduces.

In this section, we examine the steady-state grounding line flux and the grounding-line's stability as the bed strength and bed topography are varied, leaving investigations of time-evolving ice sheets to future studies. We examine first the case in which the bed varies only at the scale of the stream length, and then examine the case in which short-scale bed topography is also present.

In this section, we work with dimensional coordinates and variables. In particular, we use the dimensional form of (18), namely

(25)$$\eqalign{{C\over \rho g}q^{m + 1} + qh^{m + 1}b_x = \left({A^{1\over n}\over 4}\rho g^\prime\right)^nh^{n + m + 3}-\dot ah^{m + 2}\quad {\rm at}\ x = x_{{\rm g}},\;}$$

of (19)

(26)$$q_{{\rm S}} = \left({A\left(\rho g\right)^{n + 1}\delta^n\over 4^nC}\right)^{( {1}/( {m + 1}) ) }h^{{m + n + 3\over m + 1}} \quad {\rm at}\ x = x_{{\rm g}},\; $$

of (22)

(27)$$\eqalign{&-{b_x\over 1-\delta}\left[\left({A^{( 1/n) }\over 4}\rho g^\prime\right)^n( m + n + 3) h^{n + m + 2}-\dot a( m + 2) h^{m + 1}\right.\cr &\left.\qquad- ( m + 1) qh^mb_x\right]\cr &> \left[h^{m + 2}\dot a_{x} + \dot a\left(( m + 1) {C\over \rho g}q^{m} + h^{m + 1}b_x\right) + qh^{m + 1}b_{xx}\right]\quad {\rm at}\ x = x_{g},\;}$$

and of (24)

(28)$$q^\prime_S> \dot a.$$

To validate these approximate solutions, we compare them with numerical solutions of the problem (1), (5), (6), (8) and (10) obtained with the finite-element solver Comsol™ (COMSOL, 2021) that we term ‘exact’, acknowledging their numerical origin. These equations are solved in a steady-state form, that is, $h_t = 0$ in (5). The steady-state solutions are obtained by solving an optimization problem using a minimization procedure based on the Bound Optimization by Quadratic Approximation optimization algorithm (Powell, 2009). In all simulations, the grid resolution is spatially variable: it is 200 m through 95% of the length of the domain, and 1 m in the 5% closest to the grounding line position $x_{\rm {g}}$.

4.1. Long-wavelength bed topography

Figure 3 shows grounding line positions and the corresponding ice fluxes for the smooth, sinusoidal bed shown in Figure 3a for which the maximum magnitude of the basal slope is 1.6$\times$10$^{-3}$. Each grounding line location corresponds to a particular value of $\dot a$, which varies from 2 cm a$^{-1}$ to 14.1 m a$^{-1}$: as (25) implies, the larger thicknesses required at the deeper grounding line locations demand higher fluxes at the grounding line. Figure 3b shows the values of the driving stress and the thickness gradient at the grounding line for a set of grounding line positions and for a range of bed strengths; for $C = 7.6\times 10^6$ Pa m$^{-1/3}$s$^{1/3}$ (the value used by Schoof (2007b)) the locations are those shown in Figure 3a. Figure 3b illustrates that the thickness gradient at the grounding line is primarily determined by the strength of the bed, and that with the range of $C$ we have presented, the thickness gradient ranges from a value of the same order of magnitude as the basal topography (with a maximum slope 1.6$\times$10$^{-3}$) to one more than an order of magnitude higher. At the upper end of the range of ice-thickness gradient, we might expect Eqn (26), in which the grounding line flux is a monotonic function of thickness, to be a sufficient approximation. This is confirmed by Figure 3c, which shows flux as a function of thickness, and Figure 3d, which shows the ratio of the numerically computed flux to that predicted by (26). At the higher value of $C$, the flux has the monotonic character of (26), but, as the bed weakens, the nature of (25) and its dependence on the bed slope $b_x$ become apparent, and the grounding line thickness no longer uniquely determines the flux. Nonetheless, the disparity between (25) and (26) is not large: as Figure 3d shows, it reaches only 40% for the weakest bed with (26) consistently overestimating the flux values. Figure 3 also includes a comparison between (25) and ‘exact’ (numerically computed) solutions. In Figure 3a, the ‘exact’ grounding locations are plotted as circles over the locations determined as the roots of (25) plotted as diamonds. The difference between the grounding line positions computed numerically (circles) and as roots of expression (25) (diamonds) ranges from 20 to 920 m as the grounding line location varies from 100 to 800 km. The grounding line positions determined from (25) are significantly more accurate than those obtained from (26) for which the corresponding discrepancies are 3 and 40 km. Similarly, Figure 3c provides the same comparison of the relationship between flux and thickness at the grounding line. It is apparent that (25), and the scale choice that leads to it, provides an accurate approximation of the ‘exact’ solutions.

Fig. 3. Grounding line locations and flux with a smooth bed topography. (a) The grounding line positions for sliding parameter $C = 7.6\times 10^6$ Pa m$^{-1/3}$s$^{1/3}$ and various values of the accumulation; (b) the magnitude of driving stress, $\vert \tau _{\rm d}\vert$ (kPa) as a function of the magnitude of the ice thickness gradient, $\vert h_x\vert$, at the grounding line. Note the logarithmic scale; $\tau _{\rm d}$ and $h_x$ are both negative. (c) The relationship between ice flux and ice thickness computed for various values of the accumulation rate $\dot a$ and sliding parameter $C$. The values of $\dot a$ are discretely chosen, ranging from 2 cm a$^{-1}$ to 14.1 m a$^{-1}$. The size of symbols in panels (a) and (c) correspond to the values of $\dot a$. In panels (a) and (c), circles are ‘exact’ solutions, diamonds are the roots of expression (25). (d) Ratio of the numerically computed flux at the grounding line to the one computed with expression (26) derived by Schoof (2007a).The simulation parameters are the following: bed elevation ${b( x) = b_0 + b_{\rm a} \cos ( ( {\pi x}) /{L}) }$, with $b_0 = -500\, {\rm m}$, $b_{\rm a} = 250$ m and $L = 500$ km; ice stiffness parameter $A = 1.35\times$10$^{-25}$Pa$^{-3}$s$^{-1}$(which corresponds to ice temperature $\approx -20^\circ$C). In all simulations the ice flow is from left to right in (a).

Turning now to the stability of the grounding line positions, Figure 4a shows the grounding line positions and their corresponding stability as determined by (27). In general, for this smooth bed, (27) is in agreement with (28), the negative bed gradients are sufficient to ensure stability, while the positive bed gradients are sufficient to ensure an unstable grounding line position. However, there is one case, indicated by an arrow in Figure 4a, in which (27) and (28) disagree as to the stability. This is the consequence, in (27), of the terms in $\dot a$, which peak for the deepest grounding line, coinciding with a grounding close to the largest bed curvature. Here, although on a prograde slope, (27) predicts that the grounding line position is unstable. The accuracy of this prediction is confirmed by an ‘exact’ calculation illustrated in Figure 4b, in which a slight displacement of the grounding line upstream from its steady-state location results in an ongoing retreat.

Fig. 4. Stability of grounding line locations with a smooth bed topography. (a) Grounding line positions for various values of $C$ and $\dot a$ shown in Figure 3c. Open symbols denote stable positions ($\Lambda \lt$0), crossed symbols denote unstable positions ($\Lambda \gt$0). A red arrow indicates an unstable position that nonetheless satisfies the $q^\prime ( h( X_{\rm {g}}) ) > \dot a( X_{\rm {g}})$ stability criterion (Schoof, 2012), the unstable grounding line position indicated by the red arrow in panel (a). (b) Time evolution of an unstable grounding line position. The steady-state position marked by the black dashed line is for a sliding parameter $C = 3.8\times$10$^5$ Pa m$^{-1/3}$s$^{1/3}$ and an accumulation rate $\dot a = 9.3$ m a$^{-1}$.

4.2. Short-wavelength bed topography

The presence of the basal gradient in (25), and of the basal curvature, in (27), suggests that as the length scale of the basal topography shortens, its effects on grounding-line flux and stability will be correspondingly larger. It is the purpose of this section to investigate this hypothesis. To do so, we nonetheless need to limit any shorter scale topography to undulations whose length-scale remains large in comparison with the ice thickness, in keeping with the scale assumptions that underlie (1). Secondly, when short-scale topography is present, the longitudinal stress gradient can be as important as the basal friction and driving stress, in addition, it may no longer be appropriate to drop $B_X$ from (25) to reach (28). Therefore, the investigation must primarily depend on the ‘exact’ solutions.

Figure 5a shows a bed topography in which we have superimposed periodic undulations with a wavelength of 42 km and an amplitude of 150 m on a smooth basal topography. Figure 5a also shows the steady grounding line locations that result from a sliding coefficient of 7.6$\times$10$^6$ Pa m$^{-1/3}$s$^{1/3}$ and an accumulation rate of 0.6 m a$^{-1}$. In contrast to the smooth topography of Figure 3, the undulating bed provides many solutions for the grounding line location even at a fixed accumulation rate. The bed undulations provide local gradients of up to 3.5$\times$10$^{-2}$, and, as Figure 5b shows, the bed strength remains the dominant determinant of the driving stress at the grounding line. However, the driving stress is no longer sensitive to the thickness gradient irrespective of the sliding coefficient $C$.

Fig. 5. Grounding line locations and flux with an undulating bed topography. (a) The grounding line positions for a sliding parameter $C = 7.6\times 10^6$ Pa m$^{-1/3}$s$^{1/3}$ and an accumulation rate $\dot a = 0.6$ m a$^{-1}$; (b) the magnitude of the driving stress, $\vert \tau _{\rm d}\vert$ (kPa) as a function of the magnitude of the ice thickness gradient, $\vert h_x\vert$, at the grounding line. Note the logarithmic scale; both are negative. (c) The relationship between ice flux and ice thickness computed for various values of the sliding parameter $C$ and the accumulation rate $\dot a = 0.6$ m a$^{-1}$. In panels (a) and (a), diamonds are ‘exact’ solutions and circles are the roots of expression (18). (d) Ratio of the numerically computed flux at the grounding line to the one computed with expression (26) derived by Schoof (2007a). The simulation parameters are the following: bed elevation ${b( x) = b_0 + b_{\rm a} \cos {\pi x\over L} + b_{{\rm a}1} {\rm cos} {12\pi x\over L}- b_{{\rm a}2} \cos {24\pi x\over L}},\;$ with $b_0 = -800$ m, $b_{\rm a} = 600$, $b_{{\rm a}1} = b_{{\rm a}2} = 150$ m and $L = 500$ km; ice stiffness parameter $A = 1.35\times 10^{-25}$Pa$^{-3}$s$^{-1}$ (which corresponds to $T_{{\rm ice}}\approx -20^\circ$C). In all simulations the ice flow is from left to right in panel (a).

In general, for this bed, the bed gradient is of the same order as thickness gradient, and one might expect the distinctions between the ice flux and and that predicted by (26) to become marked. This is confirmed in Figures 5c and 5d, which show the relationship between flux and thickness at the grounding line, and the ratio of the grounding line flux from the ‘exact’ solutions to that provided by (26). As apparent from Figure 5c, the dependence of the ice flux on the ice thickness at the grounding line is much less pronounced than is the case for a smooth basal topography. Instead, the ice flux depends on a combination of the ice thickness and the bed slope. The two branches in Figure 5c ($q_{\rm g}\lt$0.005 m$^2$s$^{-1}$ and $q_{\rm g}\gt$0.01 m$^2$s$^{-1}$) correspond to the grounding line positions on the up and down slope of the underlying large-scale topography respectively. On the lower branch, the ice flux increases with increasing ice thickness, while on the upper branch it decreases. As Figure 5d shows, the flux no longer bears the simple relation with the thickness that (26) suggests. At the upper end of the bed strength, the flux ranges over values that differ by an order of magnitude from (26); as the bed weakens, this variation rises to over three orders of magnitude. Such behaviour illustrates the interplay between the bed shape (its elevation and slope) and ice flow.

Figures 5a and 5c also show the grounding line positions and fluxes determined as roots of (25). As can be seen, (25) continues to provide accurate estimates of the grounding line flux and position, even though, in these examples, the longitudinal stress gradient can in places be of equal magnitude to the basal friction. The maximum difference between the ‘exact’ grounding line positions and those determined from the roots of (25) range from 420 to 1200 m as the sliding coefficient $C$ varies from 7.6$\times$10$^6$ to 7.6$\times$10$^5$ Pa m$^{-1/3}$s$^{1/3}$. These differences are substantially larger for the grounding line positions computed with (26), whose corresponding range is 2.4–20.8 km. That (25) remains accurate reflects that, in all of these examples, the property of the smooth bed solutions – that the longitudinal stress gradient is small near the grounding line – is retained when the basal gradient and curvature become larger. As we show in Appendix C, this again arises from the non-linearity of the flow and the basal friction, and provided $\varepsilon < \delta \lt$ 1, one can expect this to be the case.

The accuracy of (25) when short-scale topography is present suggests that the stability criteria (27) will also correctly predict the stability. In Figure 6, we illustrate the stability of the grounding line locations shown in Figure 5a. Those locations for which (27) and (28) agree are shown in blue, and this is frequently the case, with 10 of the 14 cases in agreement, both stable and unstable. These cases conform to the expectation that stable grounding line positions occur on downsloping beds ($b_x \lt$ 0), unstable ones on upsloping beds ($b_x \gt$ 0). However, where grounding line locations are close to topographic maxima, and the basal curvature is large, this is no longer the case. (27) and (28) no longer agree as to their stability, and these locations are shown as green in Figure 6, with two of them, zoomed, in Figures 7a and 7b. We have verified that (27) correctly predicts the stability by examining ‘exact’, time-variant solutions that are initialised with a small perturbation from the steady states. Figure 7c shows that the green grounding line location in Figure 7a is stable even though the basal gradient is positive. Figure 7d shows that the green location in Figure 7b is unstable even though the bed gradient is negative. These examples illustrate that when short-scale topography is present, bed-slope alone is not a determinant of grounding line stability.

Fig. 6. Stability of grounding line locations with an undulating bed topography. (a) The bed topography and grounding line locations are those shown in Figure 5a. Blue open symbols are stable positions and crossed symbols are unstable positions which satisfy (22) and (24), derived by Schoof (2012). Green symbols are positions, for which stability conditions (22) and (24) produce opposite results. Zooms of regions I and II are shown in Figures 7(a–b), the red and blue triangles indicate the green symbols shown in Figures 7(a–b).

Fig. 7. Steady-state locations (a–b) and the time evolution (c–d) of perturbed grounding line positions. (a), (c) Zooms of region I in Figure 6; (b), (d) zooms of region II in Figure 6. In panels (c–d) dashed lines indicate steady-state positions. Open symbols are stable positions, crossed symbols are unstable positions which satisfy both conditions (22) and (24), green symbols are positions, for which the stability conditions (22) and (24) produce opposite results.

5. Discussion and conclusions

We have revisited investigations of marine ice sheets considered by Schoof (2007a, 2012) for which the driving stress and basal friction are of the same order, and for which the longitudinal stress gradient is, on average, smaller. We have done so to understand what role the basal topography and strength may have in determining the grounding line flux and stability. To make the role of bed topography explicit, we make a distinct scale choice by taking the buoyancy parameter $\delta$ and ice deformation parameter $\varepsilon$ to be small and of the same order, which considerably simplifies the problem by doing away with the need for a boundary layer in the vicinity of the grounding line. This approach then leads to expressions (18) and (22) which express explicitly the role of basal gradient and curvature, and in their dimensional form the basal strength, in determining the grounding line flux and stability.

These equations are derived in the limit $\varepsilon \sim \delta$. There are features of the physical situation that argue for this choice in practical circumstances, and notably that the thickness at the grounding line can be of the same order as the interior without implying run-away retreat. Figure 1 shows an example. Nonetheless, the equations derived in this limit have the property that when the thickness gradient at the grounding line is large, they reduce to the simpler forms derived by Schoof (2007a, 2011, 2012) assuming $\varepsilon \ll \delta$. As Appendix A shows, this results naturally if $\delta$ is small, but more generally it is also a consequence of the non-linearity of the ice flow and the basal friction, and particularly an ice flow law exponent with a value around 3 which, it turns out, subdues the longitudinal stress gradient even when the scaling arguments would imply it to have an important role near the grounding line; a property which carries over at the grounding line even when the undulating bed topography causes longitudinal stress gradient to be locally significant over the length of the stream.

The original hypothesis of Weertman (1974) and further investigated by Schoof (2007a, 2012) was developed for negligible bed gradients. In consequence, the grounding line flux is a monotonic function of grounding line thickness, and this results in the conclusion that stable grounding lines do not occur on positive bed gradients ($b_x > 0$). However, the Amundsen sector streams, such as the Thwaites Glacier, have topography at a wide range of scales (as is apparent in Fig. 1 and also described by Hogan and others (2020)). By taking $\varepsilon \sim \delta$, we introduce into the ice flux at the grounding line terms in the basal gradient $b_x$ arising from the momentum equation and in the accumulation rate $\dot a$ from the mass equation, and in addition introduce into the stability criteria (27) terms in their derivatives the basal curvature $b_{xx}$ and accumulation-rate gradient $\dot a_x$. This allows us to extend considerations of grounding line flux and stability to a greater range of bed geometries than considered in earlier work. The resulting expressions even appear to hold when the bed has scales of variation that are smaller than the scaling assumptions would appear to allow.

For relatively strong ($C =$ 7.6$\times$10$^6$ Pa m$^{-1/3}$s$^{1/3}$) and smooth beds ($\vert b_x\vert \sim$10$^{-3}$) the more complex flux and stability criteria do not numerically differ greatly from those of Schoof (2007a, 2012), but as the bed weakens, and/or the accumulation rate, bed gradient and curvature increase, the ice flux loses its monotonic dependence on thickness and assumes its more general implicit form (14d), which is (18) for steady-state configurations. This, in particular circumstances, can differ from the simple monotonic form by orders of magnitude. For this reason, the simplified monotonic form is not suitable as a universal parameterisation of grounding line dynamics in ice-sheet models, although it has been used that way (e.g. Pollard and DeConto, 2009; Pattyn, 2017).

Similarly, the more general stability criteria is also dependent in a complex fashion on accumulation rate, bed geometry and strength, and it is possible for stable grounding line positions to lie on up-slopes and unstable positions on down-slopes. Bed slope is no longer the determinant of stability. The gradients of the bed profile shown in Figure 1 are two-to-three times larger than the largest bed slopes considered here, suggesting that Thwaites Glacier is most likely in the regime where the flux through the grounding line is modulated by its topographic gradients, and where it is in a steady state, its stability would be determined by its bed curvature and spatial gradients of the net accumulation/ablation rate in addition to the bed slope.

The results of this study as well as of earlier studies by Chugunov and Wilchinsky (1996), Schoof (2007a, 2007b, 2012) and Tsai and others (2015) are obtained in the scaling regime $\varepsilon \ll \gamma \sim$1, that is, for marine ice streams whose longitudinal stress gradient is significantly smaller than the driving stress and basal friction and which can be neglected in the momentum balance. While the results in this paper are more general than its predecessors, for extremely weak beds a different regime is appropriate, in which $\varepsilon \sim \gamma \ll$1. This regime has been investigated by Sergienko and Wingham (2019) with the outcome that the corresponding expressions for the grounding line flux and stability are yet more complex than (25) and (27). In another important respect too this study is limited. In one dimension, the flow must go over bed undulations, whereas in two dimensions, it also has the option of flowing around obstacles (e.g. Gudmundsson, 2003; Sergienko, 2012). Particularly at shorter scales, it is questionable how informative a one-dimensional study can be. In addition, once the dimensionality of the problem is acknowledged, the implications of lateral confinement of the grounded flow, and ‘buttressing’ of the ice shelf flow substantially complicate the grounding dynamics in comparison with unconfined (one-dimensional) flow (e.g. Kowal and others, 2016; Pegler, 2018). For laterally confined configurations, the ice flux at the grounding line is an implicit function of the ice thickness (Schoof and others, 2017; Haseloff and Sergienko, 2018), and as numerical computations show (Reese and others, 2018), lead to non-physical values of the ‘buttressing parameter’ introduced by Schoof (2007b).

These previous studies on lateral confinement, our previous work considering weak beds (Sergienko and Wingham, 2019) and the present results examining the role of bed topography, together indicate that the simple form of the marine ice-sheet instability – with grounding-line flux monotonically increasing with bed depth – is the exception rather than the rule of the marine ice-sheet behaviour.