‘Stable’ and ‘unstable’ are not useful descriptions of marine ice sheets in the Earth's climate system

Abstract

Investigations of the time-dependent behavior of marine ice sheets and their sensitivity to basal conditions require numerical models because existing theoretical analyses focus only on steady-state configurations primarily with a power-law basal shear stress. Numerical results indicate that the choice of the sliding law strongly affects ice-sheet dynamic behavior. Although observed or simulated grounding-line retreat is typically interpreted as an indication of marine ice sheet instability introduced by Weertman (1974), this (in)stability is a characteristic of the ice sheet's steady states – not time-variant behavior. To bridge the gap between theoretical and numerical results, we develop a framework to investigate grounding line dynamics with generalized basal and lateral stresses (i.e. the functional dependencies are not specified). Motivated by observations of internal variability of the Southern Ocean conditions we explore the grounding-line response to stochastic variability. We find that adding stochastic variability to submarine melt rates that produced stable steady-state configurations leads to intermittently advancing and retreating grounding lines. They can also retreat in an unstoppable manner on time-scales significantly longer than the stochastic correlation time-scales. These results suggest that at any given time of their evolution, the transient behavior of marine ice sheets cannot be described in terms of ‘stable’ or ‘unstable’.

1. Introduction

A heuristic study by Mercer (1968) suggesting that the West Antarctic Ice Sheet (WAIS) collapsed at least once in its history raised the question whether WAIS is currently experiencing another collapse (e.g., Hughes, 1975; Joughin and others, 2014) or whether it could collapse in the (near) future under a warming climate (e.g., Mercer, 1978; DeConto and Pollard, 2016). Resting on a bed below sea level, the fate of the WAIS is determined by the evolution of the grounding line, ‘the junction between ice sheet and ice shelf’ (Weertman, 1974). Since the first analytical analysis of marine ice sheets by Weertman (1974), grounding-line retreat either observed on the present-day Antarctic and Greenland ice sheets (e.g., Shepherd and others, 2018; Khan and others, 2020) or simulated with ice-sheet models (e.g., Cornford and others, 2015; Seroussi and others, 2020), has frequently been interpreted as an indication of ice-sheet ‘instability’.

Realistic ice sheets are subject to the atmospheric and oceanic forcing that constantly varies in time on a wide range of temporal scales. As a result, the ice sheets do not attain steady states, stable or not. To add to the confusion, the terms ‘stable’ and ‘unstable’ are frequently used to describe the migration of grounding lines, i.e. they are used to describe unsteady, time-variable behavior (e.g., Shepherd and others, 2018a). This mischaracterization might stem from the results of Weertman's(1974) original analysis, which related the existence of a steady-state ice-sheet configuration and its stability to the bed slope at the location of the grounding line. He concluded that ‘[a] stable ice sheet can occur if the bed slopes away from the center of the ice sheet.’ (Weertman, 1974). Because the beds in many locations under the WAIS, parts of the East Antarctic Ice Sheet, and parts of the Greenland Ice Sheet slope towards the interior of these ice sheets, the grounding line positions are described as ‘unstable’ as a corollary of Weertman's result, regardless of the temporal variability of either the ice sheets or their environmental conditions.

Existing theoretical analyzes of grounding line behavior are focused on steady states of idealized configurations (e.g., Schoof, 2007a; Wilchinsky, 2009; Tsai and others, 2015; Pegler and Worster, 2013; Pegler, 2018; Haseloff and Sergienko, 2018) and questions of their stability (e.g. Schoof, 2012; Sergienko and Wingham, 2019; Haseloff and Sergienko, 2022; Sergienko and Wingham, 2022; Sergienko, 2022b). These studies assume specific functional forms of basal and lateral shear parameterizations, most often a power-law sliding law (e.g., Schoof, 2007a, 2007b; Schoof, 2012; Sergienko and Wingham, 2019, 2022) and, in the case of laterally confined (‘buttressed’) marine ice sheets, a power-law lateral shear stress (e.g., Schoof and others, 2017; Pegler, 2018; Haseloff and Sergienko, 2018, 2022; Sergienko, 2022a). Consequently, results of these studies apply only to the steady-state behavior of these particular forms of basal and lateral shear stress formulations.

To date, investigations of the temporal evolution of ice sheets, especially under changing climatic forcing and for realistic regional or continental configurations, require the use of numerical simulations (e.g., Cornford and others, 2015; Seroussi and others, 2017, 2020). The results of these numerical studies are complex and are typically interpreted in terms of the conceptual results of the theoretical analyzes, even though those are applicable to steady-state conditions only. In particular, rapid grounding line retreat is often identified as a manifestation of the instability of grounding lines on retrograde beds. Key conclusions drawn from numerical studies include: (i) rapid grounding-line retreat is triggered by loss of buttressing ice shelves (e.g., Favier and others, 2016; Gudmundsson and others, 2019), (ii) the rate of grounding-line migration strongly depends on the choice of the sliding law (e.g., Brondex and others, 2017; Sun and others, 2020; Åkesson and others, 2021), (iii) stochastic variability to climate forcing can trigger rapid grounding line retreat (e.g., Christian and others, 2022).

Reconciliation of these modeling results with our conceptual understanding of grounding line dynamics requires an extension of the existing analytical framework to the time-evolving behavior of marine ice sheets with variable basal and lateral shear stress parameterizations. To develop this framework, we derive expressions for the ice flux at the grounding line, the rate of the grounding-line migration, and stability conditions for steady-state configurations. As an illustration, we apply it to compare a laterally confined marine ice sheet subject to two different sliding laws: i) a composite viscous and Coulomb-style sliding law proposed by Zoet and Iverson (2020) on the basis of the results of laboratory experiments (hereafter referred to as Zoet-Iverson sliding law) and ii) a power-law sliding law (Weertman, 1957; Fowler, 1981). The lateral shear is described by a power law (e.g., Raymond, 1996).

In this framework, all terms of the steady-state expressions can be written in closed form. For given bed topography, basal and lateral shear stresses and environmental characteristics (i.e. surface accumulation and submarine melting), steady-state configurations can be determined (if they exist) and their stability can be established. Capturing the temporal evolution of ice sheets analytically is not possible because there is no closed form expression for a time-variant buttressing parameter (Haseloff and Sergienko, 2018). Thus, investigations of the response of buttressed marine ice sheets to temporal variations in forcing need to rely on numerical simulations to produce quantitative results, but the expressions derived here provide insights into the controls on grounding line behavior and can guide our conceptual understanding.

In this study we choose to focus on the effects of stochastic variability in submarine melting. The Earth's climate has a wide range of modes of internal variability on a variety of temporal scales. For instance, the recent retreat of the grounding line of Pine Island Glacier has been attributed to decadal oceanic variability (Jenkins and others, 2018); retreats and advances of the grounding lines in the Pacific sector of West Antarctica exhibit strong variability on this timescale as well (Christie and others, 2023). Global circulation models indicate that conditions in the Southern Ocean exhibit variability on centennial timescales (e.g., Latif and others, 2013). Numerical studies investigating the response of grounding lines to variability in climate forcing using realistic (e.g., Hoffman and others, 2019; Robel and others, 2019) and idealized configurations (e.g., Christian and others, 2022; Felikson and others, 2022) indicate a substantially different grounding-line behavior compared to that without any variability in the climate forcing.

Our results show that a classification of the grounding-line behavior into ‘stable’ or ‘unstable’ loses its meaning when considering its time-dependent behavior. We find that starting from a stable steady-state configuration, the grounding line of a buttressed marine ice sheet can exhibit irregular oscillations as well as an unstoppable retreat depending on the specifics of the stochastic variability. The results of our linear stability analysis indicate that stability conditions of steady-state marine ice sheets depend on a large number of parameters combined in a complex way, and in contrast to the results of Weertman (1974) and Schoof (2012), stability or instability of a steady state cannot be easily inferred from a simple external characteristic, such as the bed slope at the grounding line. Consequently, terms ‘stable’ and ‘unstable’ are not useful for description of the grounding line behavior of realistic ice sheets.

The manuscript is organized as follows: After a model description in section 2, we introduce the analytical framework for derivations of the expressions of the ice flux at the grounding line, the rate of the grounding line migration and the steady-state stability conditions for arbitrary forms of the basal and lateral shears (section 3). In section 4 we apply the developed framework to a marine ice sheet with Zoet-Iverson sliding law and compare its steady-states with a marine ice sheet with Weertman sliding law. These configurations are exposed to stochastic variability in submarine melting in section 5. Readers less interested in the mathematical aspects of the analysis can proceed to sections 6–7 where we summarize the results and provide their physical interpretation.

2. Model description

Following many previous studies (e.g., Schoof and others, 2017; Pegler, 2018; Haseloff and Sergienko, 2018, 2022), we consider a laterally confined ice stream flowing into a laterally confined ice shelf (Fig. 1).

Figure 1. Model geometry: b - bed elevation (b < 0), h - ice thickness, x _d - the ice divide location, x _g - the grounding line location; x _c - the calving front location.

It is described by a vertically integrated and laterally averaged 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^{w} -\tau^b -\rho gh\left(h + b\right)_x = 0,\; \quad x_d\leq x\leq x_g$$

(1b)$$2\;\left(A^{-{1}/{n}}h\left\vert u_x\right\vert ^{{1}/{n}-1}u_x\right)_x - \tau^{w} -\rho g^\prime hh_x = 0,\; \quad x_g\leq x\leq x_c.$$

Subscripts indicate partial derivatives, e.g., u _x = ∂u/∂x. u(x,  t) is the depth- and width-averaged ice velocity; h(x,  t) is ice thickness; b(x,  t) is bed elevation (negative below sea level and positive above sea level); A 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. τ ^b and τ ^w are the laterally averaged basal and the lateral shear, respectively, which can have different functional forms and dependencies on the ice stream and ice shelf. g ^′ is the reduced gravity defined as

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

where

(3)$$\delta = {\rho_w-\rho\over \rho_w}$$

is the buoyancy parameter, ρ and ρ _w are the densities of ice and water, respectively. x _d = 0 is the location of the ice divide, x _c(t) is the location of the calving front, and x _g(t) is the location of the grounding line – the location where the grounded ice starts to float and form the ice shelf. The first term on the left hand side of (1a) and (1b) is the divergence of the longitudinal stress τ _x = (2A ^−1/nh|u _x|^1/n−1u _x)_x. The basal shear τ ^b, and lateral shear τ ^w of the ice stream and ice shelf are treated as arbitrary functions, i.e.

(4a)$$\tau^b = \tau^b( \xi^b) ,\; $$

(4b)$$\tau^{w} = \tau^{w}( \xi^w) ,\; $$

where ξ ^b and ξ ^w represent a set of all variables that τ ^b and τ ^w may depend on – ice velocity, effective pressure (in the case of τ ^b) ice thickness, the ice-stream or ice-shelf width (in the case of τ ^w), etc. For brevity, we combine all these variables into one set and denote it ξ = {ξ ^b,  ξ ^w}, with understanding that the respective function (either τ ^b or τ ^w) depends on its respective set of variables (either ξ ^b or ξ ^w).

The mass balance is

(5a)$$h_t + \left(uh\right)_x \;= \dot a \qquad 0\leq x< x_g$$

(5b)$$h_t + \left(uh\right)_x \;= \dot a-\dot m \qquad x_g\leq x\leq x_c,\; $$

where $\dot a( x,\; \, t)$ is the net accumulation/ablation at the ice top surface and $\dot m( x,\; \, t)$ is the net melting/refreezing rate at the ice-shelf base (we assume that melting/refreezing at the ice-stream base is negligible).

Boundary conditions at the divide x _d and the calving front x _c are

(6a)$$( h + b) _x = 0,\; \quad u = 0,\; \ x = x_d$$

(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 x = x_c.$$

The conditions (6a) state that the driving stress at the divide is zero and there is no flow that enters or leaves from the left. The condition (6b) requires the longitudinal stress in the ice to balance the pressure deficit at the calving front.

At the grounding line x _g the continuity conditions

(7a)$$u_{{\rm stream}}( x_g) = u_{{\rm shelf}}( x_g) $$

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

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

(the last condition is on the longitudinal stress τ = 2A ^−1/nh|u _x|^1/n−1u _x), and the flotation condition

(8)$$h( x_g) = h_f = -{\rho_w\over \rho}b( x_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_w\over \rho}b( x) ,\; \quad x_d< x\leq x_g,\; $$

(9b)$$h( x) < -{\rho_w\over \rho}b( x) ,\; \quad x_g< x< x_c.$$

Haseloff and Sergienko (2018) have demonstrated that the momentum balance of the ice shelf (1b) can be integrated with the boundary condition at the calving front, x _c, (6b) and the continuity conditions (7). As a result, the problem can be split into two – one for the ice stream and one for the ice shelf. The boundary conditions at the grounding line are the flotation condition (8), velocity and thickness continuity conditions (7a)–(7b), and the stress continuity condition (7c) written in a form

(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\theta,\; \quad x = x_g,\; $$

where θ ≤ 1 is the buttressing parameter – the ratio between the backstress at the grounding line and the backstress in the absence of the ice-shelf lateral confinement, which is $\displaystyle {{1}/{2}\rho g^\prime h^2}$. This parameter represents the effects of the laterally confined ice shelf on the stress-regime of the grounding line. It encapsulates the integrated effects of the ice-shelf dynamics, the ice-shelf mass balance, other ice-shelf parameters and calving conditions on the grounding line. Hence, θ depends on the length of the ice shelf, lateral shear τ ^w and its parameters ξ, the functional forms of the net accumulation/ablation and melting/refreezing rates, $\dot a$ and $\dot m$, respectively, calving conditions (e.g., a prescribed calving-front position, a prescribed value of ice thickness at the calving front) and their parameters (Haseloff and Sergienko, 2018). It also depends on the ice flux

(11)$$q = uh$$

in the ice shelf as well as at the grounding line, q(x _g). In some cases (e.g., for particular forms of sub-ice-shelf melting $\dot m$ and calving conditions), it is possible to write a closed-form expressions of the buttressing parameter (Haseloff and Sergienko, 2018, 2022). However, there is no general expression for θ. Consequently, we treat it as an arbitrary function of a set of parameters $\varpi$ that includes all parameters mentioned above and write $\theta ( \varpi )$.

The grounded parts of marine ice sheets can experience distinct dynamic regimes. In the first one, the divergence of the longitudinal stress (the first term on the left-hand side of (1a)) has much lower magnitude than all other terms. The driving stress (the last term on the left-hand side of (1a)) is typically large, on the order of hundreds of kPa. This dynamic regime is characteristic of very wide glaciers like Thwaites Glacier in the Amundsen Sea Embayment (West Antarctica) (e.g., Holt and others, 2006), and laterally confined glaciers like Greenland marine outlet glaciers, Pine Island Glacier in the Amundsen Sea Embayment (e.g., Wingham and others, 2009), ice streams flowing into the Filchner-Ronne and Amery ice shelves (e.g., King and others, 2007), and many others. Such a stress regime has been considered by Weertman (1974); Wilchinsky (2009); Schoof (2007a, 2007b); Schoof (2012); Tsai and others (2015); Pegler (2016, 2018); Schoof and others (2017); Haseloff and Sergienko (2018, 2022); Sergienko and Wingham (2022); Sergienko (2022b). In the second regime, the divergence of the longitudinal stress is of the same order of magnitude as the other components of the momentum balance (1a). The driving stress in this regime is typically low, on the order of several kPa. This dynamic regime is characteristic of Siple Coast ice streams flowing into the Ross Ice Shelf (West Antarctica) (e.g., Rose, 1979; Bindschadler and Vornberger, 1998); it has been considered by Sergienko and Wingham (2019). Here, we focus on the first dynamic regime and provide a brief description of the treatment of the second dynamic regime in section 3.4.

3. Framework for arbitrary forms of sliding and lateral shear

As Haseloff and Sergienko (2018) showed, the behavior of a buttressed marine ice sheet is described by the momentum and mass balances of the grounded part (1a) and (5a), respectively, and the boundary conditions at the ice divide (6a) and at the grounding line: the flotation condition (8) and the stress condition (10).

A scaling analysis of Schoof (2007a, 2007b); and Sergienko and Wingham (2022) can be directly applied here, resulting in the leading order problem of the grounded part (1a), (5)–(10)

(12a)$$-\tau^{w}-\tau^b-\rho g h( h + b) _x = 0,\; \quad x_d\leq x< x_g,\; $$

(12b)$$h_t + ( uh) _x = \dot a,\; \quad x_d\leq x< x_g,\; $$

(12c)$$( h + b) _x = 0,\; \quad u = 0,\; \quad x = x_d$$

(12d)$$2 A^{-{1}/{n}}h\left\vert u_x\right\vert ^{{1}/{n}-1}u_x = {1\over 2}\rho g^\prime h^2\theta,\; \quad x = x_g$$

(12e)$$h = -{b\over 1-\delta} \quad x = x_g.$$

This is an approximate model that describes the behavior of the grounded ice stream, which is buttressed by a laterally confined ice shelf. Below we derive a relationship between the ice flux at the grounding line and dynamic (stress), geometric (bed elevation and slope) and environmental conditions at the grounding line; an expression for the grounding line migration and stability conditions for steady-state configurations in terms of the general forms of the basal and lateral shears. The only condition imposed on these general forms is that they are continuously differentiable with respect to variables they depend on. We follow derivations of Haseloff and Sergienko (2022) and Sergienko and Wingham (2022), and in contrast to studies based on the boundary-layer theory (Schoof, 2007a, 2012; Tsai and others, 2015) that consider configurations with shallow beds, negligible bed slopes and small accumulation rates at the grounding line, we do not make any additional assumptions apart from those under which the vertically integrated momentum balance (1) is derived.

3.1 Ice flux at the grounding line

Observing that

(13)$$u_x = {q_x\over h}-{q\over h^2}h_x,\; $$

where q is the ice flux (11), rearranging terms in (12a), the problem (12) can be written as

(14a)$$h_x = -{\tau^{w}\over \rho gh}-{\tau^b\over \rho gh}-b_x,\; \quad x_d\leq x< x_g,\; $$

(14b)$$h_t + q_x = \dot a,\; \quad x_d\leq x< x_g,\; $$

(14c)$$( h + b) _x = 0,\; \quad q = 0,\; \quad x = x_d,\; $$

(14d)$$q_x\rho gh^2 + q( \tau^{w} + \tau^b + \rho ghb_x) = \rho gh^{n + 3}\left({1\over 4}\rho g^\prime A^{{1}/{n}}\theta\right)^n,\; \quad x = x_g,\; $$

(14e)$$h = -{b\over 1-\delta} \quad x = x_g,\; $$

where the stress-condition (12d) written in terms of q was obtained by using (13) and (14a). We assume that θ > 0. Although it is possible to write (12d) in a form permitting a negative buttressing parameter, it is unclear how realistic such a situation could be. Physically, θ < 0 implies that the ice at the grounding line is under compression and u _x < 0. It is unclear which processes can determine the limit of compression, i.e. how negative u _x can be. Consequently, we consider only the case of θ > 0.

For a dynamic ice sheet that evolves with time, (14d) is

(15)$$\left(\dot a-h_t\right)\rho gh^2 + q( \tau^{w} + \tau^b + \rho ghb_x) = \rho gh^{n + 3}\left({1\over 4}\rho g^\prime A^{{1}/{n}}\theta \right)^n,\; \quad x = x_g.$$

As apparent from these expressions, if specific forms of the lateral shear τ ^w and the sliding law τ ^b depend on the ice velocity $\displaystyle {u = {q}/{h}}$, the ice flux at the grounding line will be an implicit function of ice thickness h and thinning rate h _t. Only if these forms are independent of the ice velocity, and the ice shelf is unbuttressed, θ=1, then the ice flux is can be written as an explicit function of the ice thickness

(16)$$q = {\displaystyle{\rho gh^{n + 3}\left(( {1}/{4}) \rho g^\prime A^{{1}/{n}}\right)^n-\left(\dot a-h_t\right)\rho gh^2}\over \displaystyle{\tau^{w} + \tau^b + \rho ghb_x}},\; \quad x = x_g.$$

If on the other hand $\displaystyle {\tau ^{w} + \tau ^b + \rho ghb_x} = 0$, then Eqn. (15) determines the grounding line position x _g(t), and the flux through the grounding line is given by $\displaystyle {q_g( t) = q( x_g( t) ,\ t) = {\int } _0^{\hskip1pt x_g( t) }( \dot a-h_t) {\rm d}x}$.

3.2 The rate of grounding line migration

The rate of grounding line migration $\dot x_g$ is determined by taking the total time derivative of the flotation condition (14e).

(17)$$\dot x_g = -{h_t + ( {b_t}/( {1-\delta}) ) \over h_x + ( {b_x}/( {1-\delta}) ) },\; $$

where the term b _t is the bed elevation changes due to, for instance, erosion/deposition processes, glaciostatic adjustment, changes in the local or global sea level. Combining this expression with (15) and rearranging terms leads to

(18)$$\eqalign{& \dot x_g = \cr& {\rho gh^{n + 3}\left({1 \over 4}\rho g^\prime A^{{1}/{n}} \theta\right)^n-\displaystyle{q( \tau^{w} + \tau^b + \rho ghb_x) -\rho g h^2\left(\dot a + ( {b_t}/( {1-\delta}) ) \right)}\over \displaystyle{h\left(\tau^{w} + \tau^b-( {\delta}/( {1-\delta}) ) \rho ghb_x\right)}}.}$$

With δ ≈ 0.1 (3), the last term in the denominator will typically be small compared to the other terms and can be neglected. This expression can also be written in terms of the ice velocity u

(19)$$\dot x_g \approx {\displaystyle{\rho gh^{n + 2}\left({1}/{4}\rho g^\prime A^{{1}/{n}}\theta\right)^n-u( \tau^w + \tau^b + \rho ghb_x) -\rho g h\left(\dot a + b_t\right)}\over \displaystyle{\tau^{w} + \tau^b}},\; $$

where τ ^w, τ ^b and θ are expressed in terms of u.

Expressions (15) and (18) reduce to previously derived expressions by Haseloff and Sergienko (2018, 2022) for Weertman sliding law in the case of buttressed marine ice sheets, by Sergienko and Wingham (2022) in the case of unconfined marine ice sheets, and in case of negligible bed slopes and accumulation rate by Schoof (2007a, 2007b).

3.3 Steady states and their stability

In a steady state, h _t = 0, and the ice-flux expression (15) is

(20)$$\rho gh^2\dot a + q( \tau^{w} + \tau^b + \rho ghb_x) = \rho gh^{n + 3}\left({1\over 4}\rho g^\prime A^{{1}/{n}}\theta\right)^n,\; \quad x = x_g.$$

The dependence of the sliding and lateral shear laws on ice velocity has consequences for the stability of steady-state configurations, and particularly whether only the sign or both the sign and the magnitude of a stability parameter (denoted here λ) can be determined. As derivations of Appendix A show, if the sliding law depends on the ice velocity, it is only possible to make inferences about the sign of λ, based on properties of a steady-state configuration at the grounding line. If λ is negative, the steady-state configuration is stable, and if λ is positive, it is unstable. These conditions imply that provided that

(21a)$$h_{x}< -\displaystyle{{b_x\over 1-\delta}}$$

(21b)$$\left(G_1( \xi) + G_2( \xi) q-F_1( \varpi) \theta_q( \varpi) \right)> 0$$

the grounding line is stable, if

(22)$$\eqalign{& {b_x\over 1-\delta}\Bigg[F_2( \varpi) -F_1( \varpi) \theta_h( \varpi) -2\dot a\rho gh-q\rho gb_x\Bigg]\cr & \quad + \dot aG_1( \xi) + qG_3( \xi) + \dot a_x\rho g h^2-F_1( \varpi) \theta_{\varpi( \varpi\neq h) }\varpi_x < 0,\; }$$

where

(23a)$$G_1( \xi) = \tau^{w} + \tau^b + \rho ghb_x$$

(23b)$$G_2( \xi) = \left(\tau^{w} + \tau^b\right)_q$$

(23c)$$G_3( \xi) = \tau^{w}_\xi\xi_x + \tau^b_\xi\xi_x + \rho ghb_{xx}$$

(23d)$$F_1( \varpi) = n\rho gh^{n + 3}\left({1\over 4}\rho g^\prime A^{{1}/{n}} \theta( \varpi) \right)^{n-1}$$

(23e)$$F_2( \varpi) = ( n + 3) \rho gh^{n + 2}\left({1\over 4}\rho g^\prime A^{{1}/{n}} \theta( \varpi) \right)^n,\; $$

and where all partial derivatives with respect to ξ or $\varpi$ indicate a sum of partial derivatives with respect to all variables that τ ^b, τ ^w or θ depend on, and ξ _x and $\varpi _x$ are the spatial derivatives of these variables; a subscript $\varpi ( \varpi \neq h)$ indicates partial derivatives with respect to all variables except h. If the conditions (21) are not satisfied, it is not possible to determine whether the grounding line is stable or unstable without solving a corresponding eigenvalue problem described in Appendix A. Condition (21a) requires that the ice stream remains grounded upstream of the grounding line, and, consequently is satisfied for any steady-state solution of interest. The second condition (21b) relates specifics of the sliding law dependence of the ice velocity and the gradient of the ice thickness. If the inequality in (22) is reversed, the grounding line is unstable.

It is possible to write these stability conditions in the same form as Schoof (2012), i.e.

(24)$${ {\rm d} q_g\over {\rm d}x_g}> \dot a,\; $$

provided

(25)$$h_{x}< -\displaystyle{{b_x\over 1-\delta}},\; \quad {\partial F\over \partial q}> 0,\; \qquad {\partial F\over \partial q_x}> 0,\; $$

where (20) is written in terms of an implicit function F

(26)$$\eqalign{F( q,\; \dot a,\; h,\; x,\; \cdot) \, \colon = \dot a \rho gh^2 + q( \tau^{w} + \tau^b + \rho ghb_x) \cr- \rho gh^{n + 3}\left({1\over 4}\rho g^\prime A^{{1}/{n}}\theta\right)^n = 0}$$

with ( ⋅ ) being a placeholder indicates dependencies on other variables than (q,  q _x,  h,  x). A derivation of this form is described in Appendix A.

In some special cases, it is possible to obtain an analytic expression for the stability parameter λ (see Appendix A). The marine ice sheet is stable if λ < 0 and unstable if λ > 0, and the magnitude of λ allows to compute the e-folding time t _e = |λ|⁻¹ (Sergienko and Wingham, 2019). We provide a closed-form expression for λ in (A.15) for the case where both τ ^b and τ ^w are independent of u simultaneously. While there are formulations of the basal shear τ ^b independent of u that are based on the arguments that subglacial till can behave plasticaly, (e.g., Bassis and Ultee, 2019), the existing formulations for τ ^w depend on ice velocity (e.g., Raymond, 1996). Potentially, an argument that the lateral shear could be described by the yield stress in a similar way as the basal shear could be made (e.g., Thomas, 1977), however, such formulations have not being tested in numerical models.

The stability conditions (22) are derived under the assumption that there are no feedbacks between the marine ice sheet and its environmental conditions, for example, the surface net accumulation $\dot a$ and submarine melting $\dot m$ do not depend on characteristics of the ice sheet (e.g., its surface elevation or thickness), but can be functions of the ice-sheet horizontal extent and other parameters. If the feedbacks between the marine ice sheets and their environmental conditions are taken into account (e.g., Fyke and others, 2018) then it is not possible to derive stability conditions that are based on the characteristics of the steady-state configurations and the environmental conditions (Sergienko, 2022b).

3.4 Scope of the framework

The derived expressions and conditions can be applied to unconfined and buttressed marine ice sheets as long as the momentum balance of the grounded part can be approximated by (12a) with sufficient accuracy. For an unconfined or very wide ice shelf, θ → 1 and the behavior of the marine ice sheet is determined by the properties of the grounded part only. In the case of a laterally unconfined marine ice sheet, in addition to θ = 1, lateral shear stresses are negligible, i.e. τ ^w = 0 in all above expressions.

The approximate expression of the momentum balance (12a) is valid if the driving stress is balanced by the basal and lateral shears and the divergence of the longitudinal stress is smaller than these components. If, for instance, the ice-stream bed is very weak, and τ ^b is significantly smaller than the lateral shear and driving stress, then τ ^b can be neglected. However, conditions where both the basal and lateral shear stresses are of the same order of magnitude as the divergence of the longitudinal stress cannot be treated with this developed framework. For ice streams flowing into unconfined ice shelves in such a stress regime, the ice flux at the grounding line and stability conditions can be obtained from those derived by Sergienko and Wingham (2019), their eqn (29) and eqn (36), respectively, where in the first-order solution h ₁(x) eqn (30), the last term $\displaystyle {C( {q}/{\vert b\vert }) ^m}$ is replaced with τ ^w + τ ^b in which the dependence on the ice velocity is expressed in terms of the ice flux, i.e. as $\displaystyle {u = {q}/{h_0}}$, where h ₀ = −b.

The only case that is not covered either by the framework presented here or the analysis of Sergienko and Wingham (2019), is the case in which the leading order balance is between the divergence of the longitudinal stress and the driving stress, and the lateral and basal shear stresses are negligible (either through the length of the ice stream or in the vicinity of the grounding line). This case corresponds to a laterally unconfined marine ice sheet with no basal or lateral resistance in the vicinity of the grounding line; it requires a separate treatment. A potential approach to develop such a treatment could be based on using the approximate form of the momentum balance (12a) far away from the grounding line, the approximate form of the momentum balance considered by Sergienko and Wingham (2019)(Eqn. (19)) with Γ → 0 (where Γ is a non-dimensional parameter related to the basal shear), and matching these formulations by means of the matched asymptotic analysis.

4. Application to Zoet-Iverson and Weertman sliding laws

As an example of the use of the framework described above, we apply it to specific forms of the lateral and basal shear stresses. We choose the same lateral-shear form as used in previous studies of laterally confined marine ice sheets and outlet glaciers (e.g., Raymond, 1996; Schoof and others, 2017; Haseloff and Sergienko, 2018, 2022; Sergienko, 2022a)

(27)$$\tau^w = {C_w A^{-1/n} \over W^{1/n + 1}} h \vert u\vert ^{1/n-1}u,\; $$

where C _w = 2(n + 1)^1/n is a lateral shear stress parameter and W is the width of the marine ice sheet and ice shelf. (The derivation of this form can be found in Raymond (1996).) Here, we use the same formulation of the lateral shear for the grounded part and the ice shelf. All parameters are listed in Table 1.

Table 1. Model parameters

.

Based on laboratory experiments, Zoet and Iverson (2020) have proposed a sliding law of the form

(28)$$\tau^b = N\tan\phi\left({1\over u + u^t}\right)^{1/p}\vert u\vert ^{1/p-1}u,\; $$

where N is the effective pressure, ϕ is the till friction angle, u ^t is a transition speed at which basal shear transitions from the viscous-like deformation of underlying till to Coulomb style, and p ~ 5 is a slip exponent. The transition speed u ^t is

(29)$$u^t = \alpha N$$

with $\alpha = {N_f}/( 2 + N_fk) ( 1/( {\eta ( Ra) ^2k_0^3}) + {4C_1}/( ( Ra) ^2k_0) )$ determined by the geological and rheological conditions of the underlying till (N _f is the till-bearing capacity factor; k is a pressure-shadow factor; η is the effective dynamic viscosity of ice; R is the clast radius, a is the fraction of the clast radius that extends above the ice-till interface; C ₁ = C_pK/L is the regelation parameter, where C _p is the depression of the melting temperature of ice with pressure, K in the mean thermal conductivity of ice and rock, and L is the volumetric latent heat of ice, and k ₀ = 2π/(4R), with 4R chosen to approximate the wavelength of bumps in Nye's model). Using the values listed for these parameters in the Supplemental Materials of Zoet and Iverson (2020), we obtain α ≈ 1.1 × 10⁻¹¹ Pa⁻¹ m s⁻¹.

As Haseloff and Sergienko (2022) point out, in limited cases it is possible to construct analytic expressions for the buttressing parameter θ. For instance, for submarine melting that depends on the horizontal position x

(30)$$\eqalign{\theta = 1 - \left[{( 4^n C_w) ^{1/( 1 + n) } [ q( x_c) ] ^{1/n} + \displaystyle{( ( {n + 1}) /{n}) {C_w}/{W} {\int}_{x_g}^{\,x_c}[ q( x') ] ^{1/n} {\rm d}x'}\over \rho g^\prime A^{1/n}W^{1/n}h^{1/n + 1}} \right]^{2n/( n + 1)}}$$

where

(31)$$q( x) = q_g + { \int }_{x_g}^ {\hskip1pt x} \left(\dot a-\dot m\right){\rm d}x^\prime .$$

The effective pressure N is approximately the difference between the ice overburden and the subglacial water pressure

(32)$$N\approx\rho g h-p_w.$$

At the grounding line x = x _g, the subglacial bed is in contact with the ocean, and we assume

(33)$$N( x_g) \approx 0.$$

Zoet and Iverson (2020) performed their experiments for constant values of N = 136 kPa and N = 158 kPa. To investigate the effects of reducing N to zero at the grounding line, we use an approach similar to that used by Joughin and others (2019), and linearly reduce N from a constant value (136 kPa) to zero in a zone (10 km and 50  km) upstream from the grounding line:

(34)$$N( x) = \left\{\matrix{N_0 \hfill & {\rm if } \, {x} < {x}_{g}- {L}_{{N}_{0}} \hfill \cr N_0 \displaystyle{\left(1 - { {x}_g-x\over L_{N_0}}\right)} \hfill & {\rm if } \, {x}_{g}- {L}_{N_0}\leq {x} < {x}_{g}. \hfill }\right.$$

We refer to $L_{N_0}$ as the reduction length scale.

We compare the steady-state and dynamic behavior of a marine ice sheet with the above described forms of lateral and Zoet-Iverson basal shear stress to the behavior of a marine ice sheet with the same lateral shear and Weertman sliding law

(35)$$\tau^b = C\left\vert u\right\vert ^{m-1}u$$

that has been investigated in earlier studies (Haseloff and Sergienko, 2018, 2022). We do so by using the constructed analytic expressions and numerical simulations of the steady-state formulation of the full model (1)–(8), i.e. we set h _t = 0. To be consistent with the analysis of the marine ice-sheet behavior with Zoet-Iverson sliding law, we consider a linear reduction of the Weertman sliding parameter C from its constant value to zero in the 10 km and 50 km zones upstream of the grounding line (as in Eqn. (34) but with N replaced by C). We use the bed shape

(36)$$b( x) = b_0 + b_2\tilde x^2 + b_4\tilde x^4 + b_6\tilde x^6,\; $$

where $\tilde x = x/L_s$, which is based on the shape used in MISMIP+ experiments (Asay-Davis and others, 2016), with minor differences in the coefficients b ₀ and L. For the submarine melt rate we use

(37)$$\dot m( x) = \dot m_0\left[1-\left({x-x_g\over L_c-x_g}\right)^{1/3}\right] + \dot a,\; \quad x_g< x \leq L_c.$$

Melting is largest at the grounding line, and at the calving front L _c the net ablation/accumulation $\dot a-\dot m$ is always zero; this form is chosen purely for its simplicity.

Numerical simulations that solve the full model (1)–(8) are performed with the finite-element solver Comsol^TM(COMSOL, 2023). 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 _g. The high resolution in the vicinity of the grounding line is used to ensure its accurate position where both the flotation and the stress conditions are satisfied simultaneously.

4.1 Steady-state configurations

For τ ^w and τ ^b given by (27) and (28), the expression for the ice flux at the grounding line (20) is

(38a)$$\dot a\rho gh^2 + q\left[{C_w A^{-1/n} \over W^{1/n + 1}} {q^{1/n}\over h^{1/n-1}} + N\tan\phi\left({q\over q + hu^t}\right)^{1/p} + \;\rho ghb_x\right] = \rho gh^{n + 3}\left({1\over 4}\rho g^\prime A^{{1}/{n}}\theta\right)^n.$$

If N(x _g) = 0 then this expression simplifies and becomes

(38b)$$\dot a\rho gh^2 + q\left[{C_w A^{-1/n} \over W^{1/n + 1}} {q^{1/n}\over h^{1/n-1}} + \rho ghb_x\right] = \rho gh^{n + 3}\left({1\over 4}\rho g^\prime A^{{1}/{n}}\theta\right)^n.$$

These expressions are the steady-state stress conditions that are satisfied together with the flotation condition (8) at the junction between the grounded and floating parts of the marine ice sheet. They are implicit relationships between the steady-state ice flux and ice thickness. According to the mass balance (5a), in a steady state

(39)$$\displaystyle{q( x_g) } = { \int }_{x_d}^{{\hskip1pt x}_g}\dot a {\rm \, d}x.$$

With this condition (39), (38a) and (38b) are transcendental equations whose roots determine the grounding line positions, of which there can be none, one or several.

The results of the analysis by Haseloff and Sergienko (2022) show that the choice of the conditions at the calving front has a strong effect on the steady-state configurations of marine ice sheets. For the chosen bed shape (36) and the values of parameters described in Table 1, there are unique steady-state positions on the up-sloping part of the bed with a fixed position of the calving front (Figs. 2(a,b)). Out of three calving-front conditions considered by Haseloff and Sergienko (2018, 2022) this is the only one that yields steady-state configurations for this bed shape and the parameter set. With both sliding laws the ice sheets attain similar steady-state configurations. The grounding-line position for the case with zero basal shear is upstream of the grounding-line position with non-zero basal shear. The upstream displacement of the grounding line position for the 10 and 50 km reduction zones are fairly similar – 5.7 km for the 10 km and 6.7 km for the 50 km reduction zones for both sliding laws (6.5 km for Weertman's). The differences between the ice thicknesses of the grounded ice streams computed with the approximate and the full model do not exceed 3 m in all simulations. The differences between the ice-shelf thicknesses are larger: within 20 m directly downstream of the grounding line, and they progressively reduce to under 2 m towards the calving front.

Figure 2. Steady-state configurations of ice sheets. (a) profiles with the Zoet-Iverson sliding law with N reducing to 0 in the 50 km upstream of the grounding line (green line), 10 km (red line) and constant through the length of the ice stream; (b) profiles with the Weertman sliding law with C reducing to 0 in the 50 km upstream of the grounding line (green line), 10 km (red line) and constant through the length of the ice stream. Solid lines are numerical simulations of the full model, dashed lines are results of the approximate model (14). (c) The ice flux at the grounding line as a function of the ice thickness. Lines are solutions of the approximate analytic expression (38a) with the corresponding sliding laws; τ _b(x _g) = 0 is a solution of (38b); symbols indicate steady-state values from the numerical simulations: filled symbols indicate Zoet-Iverson sliding law open symbols indicate Weertman sliding law, colors indicate the reduction of τ ^b in the zone upstream of the grounding line similar to those shown in panels (a) and b. (d) The ice flux at the grounding line as a function of the position x _g. The light magenta, light green and turquoise lines are solutions of (38a) and (38b), the black line is the integrated accumulation rate through the grounding line (39); symbols are the same as in panel (c).

One of the advantages of the analytic expressions (38) is the possibility to examine the steady-state ice flux for a range of possible grounding line positions, see Figures 2c,d. The solid lines are obtained by solving (38) with the flotation condition (8) for x _g between 350 km and 550 km. Due to the appearance of b _x and θ in (38), these curves depend on the ice shelf and bed geometry. For $x_g \lesssim$ 350 km all ice on the ice shelf is removed by melting and the ice shelf front is located where the ice thickness goes to zero. We exclude this case from our analysis.

In the presence of a sufficiently long ice shelf (≈100 − 250 km, $x_g\lesssim 475$  km) the form and value of the sliding law do not greatly affect the steady-state ice flux at the grounding line. However, the differences are significantly more pronounced as the ice shelf becomes shorter and buttressing is smaller (Fig. 2(d)).

In steady state, the ice flux at the grounding line matches the integrated surface accumulation (39), and the steady-state grounding line position is determined by the intersection of the flux determined from (38) and the mass balance flux given by (39) (black line in Fig. 2(d)); these intersections are marked by the symbols in Fig. 2(d). The values obtained with the approximate model (14) are in a good agreement with the values obtained in numerical solutions of the full model: the differences between positions obtained with the numerically solved full model (1)–(9) and the roots of (38) are between 650 m a 1200 m for both sliding laws, with smaller values for the configurations with the non-vanishing basal shear and larger values for the configurations for the 50 km reduction zone. The reduction of τ ^b to zero at the grounding line results in a small reduction (less then 2%) of the steady-state ice flux through the grounding line Figure 2(c)). This reduction is due to the fact the grounding line position is slightly upstream, and in a steady state the flux at the grounding line is the integral of the surface accumulation rate $\dot a$ from the divide to the grounding line; with the constant values of $\dot a$, this integral is smaller for the grounding line position located upstream than for the grounding line position located downstream.

As the basal and lateral shear stress depend on the ice velocity, the stability of the steady-state configurations can be determined from (21)–(23). The expressions for the stability coefficients for the Zoet-Iverson sliding law are listed in Appendix B; the expressions for the Weertman sliding law can be found in Haseloff and Sergienko (2022)(Eqn. (C20)). According to (22), all configurations shown in Figs. 2(a,b) are stable. This is confirmed with numerical simulations in the same way as in similar studies by Sergienko and Wingham (2019, 2022) and Sergienko (2022b) by displacing the grounding line positions 1 km upstream from their steady-state positions and solving the full time-dependent problem (1)–(8). In all simulations, the displaced grounding lines return to their steady-state positions, confirming that the steady-state configurations are stable.

Expressions (15)–(20) were derived under assumptions that on the grounded part the divergence of the longitudinal stress is smaller than other components of the momentum balance. It is instructive to analyze the components of the ice stream (1a) and ice shelf (1b) momentum balances in order to verify our assumption that for the ice stream the divergence of the longitudinal stress is smaller than other components of the momentum balance is valid. We use the terms from the numerical simulations with the full model (Fig. 3).

Figure 3. Terms of the steady-state momentum balance. (a) Zoet-Iverson sliding law; solid lines correspond to a constant value of N, dashed lines correspond to reduction of τ ^b to zero in the 10 km upstream of the grounding line, dotted-dashed lines correspond to reduction of τ ^b to zero in the 50 km upstream of the grounding line; (b) Weertman sliding law, the lines are the same as in panel (a); (c) zoom-in of τ _x terms in the vicinity of the grounding line for Zoet-Iverson sliding law; (d) the same as panel (c) for Weertman sliding law.

The driving stress τ ^d (last terms on the left-hand side of Eqns. (1a)–(1b)) is discontinuous at the grounding line. This discontinuity results in a discontinuity of the divergence of the longitudinal stress τ _x at the grounding line, even though the longitudinal stress itself is continuous in accordance with the continuity condition (7c). The lateral shear stress is continuous at the grounding line because of the continuity conditions (7a)–(7b). As Haseloff and Sergienko (2018) surmised, there is a boundary layer downstream of the grounding line, in which the longitudinal stress divergence rapidly reduces to zero, and outside this boundary layer the momentum balance of the ice shelf is between the driving stress and the lateral shear stress.

Upstream of the grounding line, the dominant balance is between the basal shear, lateral shear and driving stress. The imposed reduction of the basal shear is accompanied by the reduction of similar magnitude in the driving stress, rather than significant increase of the magnitude of the longitudinal-stress divergence. It is only in the immediate vicinity of the grounding line that τ _x reaches a magnitude of 1–2 kPa, with the value slightly larger for the 10 km reduction zone than for the 50 km reduction zone. This relative smallness of the longitudinal-stress divergence provides justification for the approximated momentum balance (14a). Both sliding laws result in a similar behavior of the stress terms. This behavior of the momentum-balance terms is qualitatively similar to that obtained by Tsai and others (2015) in their analysis of an unconfined marine ice sheet with Coulomb sliding law. The linear stability analysis, and hence a derivation of stability conditions of their specific formulation of the sliding law is problematic because it is written in terms of a min function, and its derivatives may not exist.

5. Response to the stochastic variability in submarine melting

While most previous theoretical studies investigate steady-state configurations that assume that all internal and external parameters remain constant in time, realistic ice sheets experience environmental conditions that change in time. This temporal variability happens on a wide range of timescales, because in addition to the long-term changes in the orbital forcing (Milanković, 1941), the Earth's climate has numerous modes of internal variability. Steig and others (2012) linked the variability in the tropical Pacific to variability in submarine melting of the ice shelves in the Pacific region of West Antarctica, and Jenkins and others (2018); Christie and others (2023) have attributed the grounding line migration of these ice shelves to decadal variability. The results of numerous observational (e.g., Shepherd and others, 2004) and modeling (e.g., Cornford and others, 2015; Seroussi and others, 2017) studies indicate that the grounding-line dynamics is sensitive to submarine melting. Consequently, we chose to explore how the marine ice sheets respond to stochastic variability in submarine melting. Previous modeling studies with idealized (Christian and others, 2022; Felikson and others, 2022) and realistic (Hoffman and others, 2019; Robel and others, 2019) geometries show that stochastic variability in the climate forcing have strong effects on the long-term behavior of the grounding lines.

We use the steady-state configurations obtained in the previous section as initial conditions; as discussed above they are all stable configurations. We modify the expression for the submarine melt rate (37) in the following way

(40)$$\dot m( x,\; t) = \Bigg(\dot m_0 + \Delta\dot m( t) \Bigg)\left[1-\left({x-x_g\over L_c-x_g}\right)^{1/3}\right] + \dot a,\; \quad x_g< x\le L_c,\; $$

where $\Delta \dot m( t)$ varies stochastically

(41)$$\Delta\dot m( t) = \dot m_{s} {\cal N}( {t}/{T_{s}}) + \dot m_{l} {\cal N}( {t}/{T_{l}}) .$$

Subscripts s and l refer to ‘short’ and ‘long’ time-scales, $\dot m_{s, l}$ are the amplitudes of the melt-rate variability; T _s,l are the correlation time-scales, and ${\cal N}( t)$ is a noise function with a uniform distribution and zero mean value. For the ‘short’ time-scale we choose the decadal and for the ‘long’ time scale we choose centennial time scales as suggested by in situ and remote sensing observations (Jenkins and others, 2018; Christie and others, 2023) and modeling (Latif and others, 2013).

For each sliding law we perform four sets of simulations (see table 2): with constant values of N or C or with a reduction in the 10 km zone upstream of the grounding line; and with T _l = 200 yrs or T _l = 500 yrs. The values of m _s,l are the same in all simulations. Each set consists of five simulations with the noise function ${\cal N}( t)$ initialized with a random seed to ensure independence of the simulation results.

Table 2. Performed time-dependent simulations

X indicates the occurrence of an irreversible ice sheet collapse, 0 indicates no ice sheet collapse within the 100 kyrs simulation period. For T _l = 200 yrs and T _l = 500 yrs, 5 simulations were performed for each configuration. Ice sheet collapses occur at different times during these simulations, depending on the seed in the random noise function ${\cal N}$

The problem (1)–(8) is solved numerically on domains with a moving boundary (the grounding line) using an arbitrary Lagrangian-Eulerian (ALE) method (Donea and others, 2017). The time-step is adaptive and does not exceed two years in order to resolve variability on T _s = 10 yrs. The grounding line moves with velocity (17), where all terms are computed by the model and b _t = 0. The analysis relies on numerical simulations because it is not possible to write a closed-form expression for θ ((30) is only valid for steady-states). Instead, θ is determined numerically from (10).

Simulations with the Zoet-Iverson sliding law, the 10 km reduction zone and T _l=200 yr indicate that stochastic variations in the submarine melting can cause migration of the grounding line on the order of ~40 km with excursions to almost 70 km (Fig. 4(a)). Although the time-averaged values of the rate of the grounding-line migration $\dot x_g$ over the 100 kyr simulation period are low (~−4.5 m yr⁻¹, Fig. 4(b)), there are infrequent instances when the rate of the grounding-line retreat exceeds 1 km yr⁻¹. Together in all five simulations that happened 27 times. The average periodicity of such events is ~18.5 kyr – significantly longer then the long correlation time-scale T _l = 200yrs.

Figure 4. Response of the grounding line to stochastic perturbation in the melt rate with T _l = 200 yrs (left column) and T _l = 500 yrs for Zoet-Iverson sliding law and with $L_{N_0} = 10$ km. (a) the grounding line position (km); (b) the rate of the grounding line migration (m yr⁻¹); (c) the buttressing parameter θ. Different colors correspond to simulations with different seeds in the noise functions. Dashed lines indicate steady-state values.

In steady state, melting reduces the ice-flux divergence ($\dot a-\dot m = q_x$), which can lead to both thinner and slower ice downstream of a melt perturbation. As the steady-state buttressing parameter depends on the integrated ice flux (see (30)), in either case melting leads to a reduction in buttressing if everything else is fixed. In time-dependent calculations as those performed here, melting can either reduce the ice-flux divergence or lead to thinning ($\dot a - \dot m = q_x + h_t$) and it is not obvious how the partition between these two processes affects the buttressing parameter and consequently the grounding line evolution. Moreover, as the location of the calving front is fixed in our simulations, a retreating grounding line lengthens the ice shelf, which could in turn increase the amount of buttressing experienced at the grounding line, with the counter-intuitive result that grounding line retreat and an increase in buttressing go hand-in-hand. This is indeed what we observe: the temporal evolution of the buttressing parameter θ is shown in Figure 4(c) and θ decreases (buttressing increases) as the grounding-line retreats, and conversely. Simulations with the same parameters and no reduction of τ ^b, as well as simulations with a Weertman sliding law with and without reduction of τ ^b produce quantitatively similar results.

Increase of the long correlation time-scale T _l to 500 yr produces markedly different results in simulations with and without reduction of τ ^b near the grounding line. If the basal shear stress is reduced at the grounding line ($L_{N_0} = 10$ km) then the marine ice sheet collapses in all simulations, i.e. the grounding line retreats to the ice divide at x _d = 0, see Figure 4(d). The timing of this collapse is significantly longer (~5 kyr to ~75 kyr) than the correlation time-scale T _l of stochastic variability of submarine melting (500 yr). The rates of grounding line migration prior to the ice-sheet collapse are lower than in simulations with T _l =200 yrs. During the collapse they exceed 1 km yr⁻¹ (Fig. 4(e)). Similar to the simulations with the shorter correlation time-scale T _l =200 yrs, prior to collapse, the buttressing parameter θ increases when the grounding line retreats, and decreases when the grounding line advances. It is only during the collapse that θ significantly increases indicating strong loss of buttressing (Fig. 4(f)). Simulations with the Weertman sliding law and $L_{N_0} = 10$ km exhibit quantitatively similar behavior – the ice sheet collapses in ~8 kyr – 120 kyr depending on the simulation.

If τ ^b does not vanish at the grounding line and T _l = 500 yrs, the marine ice sheet does not undergo collapse; instead, it experiences an intermittent advance and retreat of the grounding line. In these simulations, the rate of the grounding line migration is overall lower than in simulations with T _l = 200 yr. Only in one simulation (for each sliding law) does the rate of the grounding line migration exceeds 1 km yr⁻¹. However, if the correlation time scale T _l is increased to 1000 yrs, then a collapse-like behavior is observed in simulations with both Zoet-Iverson and Weertman sliding laws and constant values of N and C as well. The ice-sheet behaviors in these simulations are similar to those in simulations with τ ^b vanishing at the grounding line and T _l = 500 yrs, described above. Table 2 summarizes these results.

To get a better understanding how the values of the correlation time-scales affect the dynamics of the grounding line, we analyze the co-evolution of the melt rate, the grounding line displacement and the buttressing parameter in single simulations with different values of T _l. As Figs. 5(a,b) and 5(d,e) show, in general, the grounding line retreat coincides with a rapid increase of the melt rate and vice versa. Changes in θ overall mimic changes in the grounding-line position: when the grounding line retreats, θ decreases (buttressing increases), and when the grounding line advances, θ increases (buttressing decreases). As pointed out above, this behavior of θ is a consequence of the chosen conditions at the calving front and the fixed values of the surface ablation/accumulation rate $\dot a$. The difference between simulations with different values of T _l is that in simulations with T _l = 500 yrs the ice shelf can experience melt rates substantially larger than their mean values for longer time periods compared to simulations with T _l = 200 yrs. The longer periods of sustained high melt rates imply longer periods during which the grounding line is upstream of its time-averaged position, and consequently shorter horizontal extent of the grounded part of the ice sheet, through which it accumulates mass.

Figure 5. Response of the grounding line to stochastic perturbations in the melt rate in a single simulation for Zoet-Iverson sliding law with the 10 km reduction zone (a)–(c) for T _l = 200 yrs and (d)–(f) for T _l = 500 yrs. (a), (d) melt rate at the grounding line (m yr⁻¹); (b), (e) displacement of the grounding line from the time-mean position (km); (c), (f) the buttressing parameter θ. Note different horizontal scales in panels (a)–(c) and (d)–(f).

Ultimately, the unstoppable retreat of the grounding line is caused by the loss of the ice flux through the floating ice shelf exceeding the mass gain upstream of the grounding line for prolonged periods of time (Fig. 6). This leads to continuous thinning and further grounding line retreat, until the ice sheet has completely disappeared.

Figure 6. Evolution of the ice flux q(x) through the length of the ice sheet during unstoppable retreat for Zoet-Iverson sliding law with the 10 km reduction zone for T _l = 500 yrs. Dashed lines indicate the mass flux gained upstream of the grounded part.

6. Discussion

The existing theoretical analyzes have focused thus far only on steady states of marine ice sheets and their stability. In order to improve our conceptual understanding of the dynamic behavior of marine ice sheets and extend the theoretical analyzes to time-variant behavior that is characteristic of realistic ice sheets experiencing environmental conditions that vary in time, we have developed a framework to analyze the effects of different forms of basal and lateral shear stresses on marine ice sheet dynamics. This framework can be applied to analyze the role of lateral and basal shear stress parameterizations in steady state configurations. Written with general forms of τ ^w and τ ^b, the derived expressions also provide insights into how basal and lateral shear (or their absence) affect the temporal evolution of grounding lines. The stability conditions presented here are derived assuming no feedbacks between ice sheet and surface accumulation or submarine melting. Even in the absence of these additional feedbacks, the stability conditions depend on a large number of parameters in a non-trivial fashion. This, and the presence of feedbacks between realistic ice sheets and their environmental conditions, questions generalized statements about their (in)stability and the usefulness of such statements.

The rate of grounding line migration (18) (or (19)) is determined by the imbalance between the stresses on the grounded side of the ice stream (the basal, lateral and form drag associated with the shape of the bed), the ice-shelf buttressing, whose effect is encapsulated in the buttressing parameter θ, the ice velocity at the grounding line, which reflects the effects of lateral and basal shear through the length of the grounded part (19) and the net accumulation/ablation at the grounding line. In steady state these terms balance each other, $\dot x_g = 0$, and the grounding line does not move. The fact that $\dot x_g$ is determined by the imbalance of several terms reflects that it is controlled by the competing effects of several processes – both local, at the grounding line, as well as integral, through the lengths of the grounded part and the ice shelf. This indicates that the observed migration of present-day grounding lines cannot be attributed to a single cause (e.g., submarine melting), and all other processes and interactions between them have to be taken into account.

To illustrate the developed framework, we have applied it to a marine ice sheet whose lateral shear increases with velocity and ice thickness and whose basal shear is described by either a Zoet-Iversion sliding law or a Weertman sliding law. The only calving law that yields a steady-state configuration for the chosen bed topography is a fixed calving front position – a widely used condition in idealized and realistic studies. For both sliding laws, we find a good agreement between the steady-state configurations obtained with the full and the approximate models. Analysis of the momentum balances of these ice sheets shows that the divergence of the longitudinal stress at the grounding line is discontinuous (Fig. 3) because the driving stress at the grounding line is discontinuous (Eqn. (1)). Although the discontinuity of the driving stress is obvious from the momentum balance (1), the discontinuity of the divergence of the longitudinal stress is less so. Its existence should be kept in mind in analyzes of the large-scale ice-sheet simulations that are typically based on the vertically integrated momentum balances. Most likely, the requirements of high spatial resolution in the vicinity of the grounding line (Seroussi and others, 2014; Gladstone and others, 2017) are associated with the need to resolve this discontinuity of the divergence of the longitudinal stress with sufficient accuracy.

The steady-state configurations with the two sliding laws we compare do not differ qualitatively. There are also no substantial differences in the steady-state configurations of ice sheets whose basal shear reduces to zero or remains large at the grounding line, regardless of the sliding law.

However, if the grounding line is dynamic and driven by stochastic variability in submarine melting, then the ice sheet response depends on the basal shear stress at the grounding line and the magnitude of the correlation time-scales. Our results show that on the given bed topography, the marine ice sheet can exhibit intermittent advance and retreat of their grounding lines in response to up to 40% variations compared to steady-state melt rates if such variations are sustained on the order of 200 yrs. If such melt-rate variations are sustained for longer periods (~500yrs) simulated ice sheets with vanishing basal shear experience unstoppable retreat of their grounding lines, while simulated ice sheets with large basal shear at the grounding lines continue to exhibit intermittent advance and retreat of their grounding lines. They experience such an unstoppable retreat if melt-rate variations are sustained for longer periods (~1 kyr). We point out that the exact values of periods are specific to the chosen parameters - the bed shape, the magnitudes of the basal and lateral shears and the amplitudes of the melt-rate variations. For instance, for the amplitudes of variations much larger than the steady-state values, the correlation timescales may be shorter, e.g., of the order of several decades. However, our finding that for the same amplitudes of the melt-rate variability, the longer correlation timescales can eventually lead to unstoppable retreat of the grounding line is a robust result.

The unstoppable retreat observed in our simulations is a consequence of an excess mass loss at the grounding line which exceeds the integrated mass gain upstream of the grounding line for prolonged periods of time, leading to sustained thinning and retreat of the ice sheet. The retreat occurs at times significantly longer than the stochastic correlation time-scale, on the order of thousands to hundreds of thousands of years, after many periods of melt-rates higher than their time-averaged values. During intermittent grounding-line retreats, the ice-shelf buttressing increases and vice versa. This behavior is a result of holding the calving front position fixed.

Even though our assumption that all parameters except submarine melt rates remain constant in time is unrealistic, it simplifies the analysis and allows us to focus on the effects of time-varying melting. A drastically different response to temporal variability in submarine melting observed in our simulations with vanishing and non-vanishing basal shear stress at the grounding line suggests that the basal traction in the vicinity of the grounding line exhibits a strong control on the transient behavior of marine ice sheets. It also suggests that the temporal variability of basal conditions can have a critical role in its dynamics. Although inferences about the basal shear spatial distribution from surface observations by means of inverse modeling have been done routinely (e.g., Sergienko and others, 2008; Morlighem and others, 2013), there is very little understanding of its temporal variability. The presence of high basal shear structures in the vicinity of the grounding lines of many Antarctic ice streams and Greenland outlet glaciers has been attributed to the interactions between subglacial hydraulic systems, sediments and the ice flow (Sergienko and Hindmarsh, 2013; Sergienko and others, 2014). The highly dynamic nature of the subglacial hydraulic system means that assuming sliding law parameters do not change in time is very limiting. Establishing characteristics of their variability is equally important as establishing a functional form of the sliding law, and is necessary for improvements in understanding of grounding line dynamics.

The concepts of marine ice-sheet stability and instability introduced by Weertman (1974) and further developed by Schoof (2007a, 2007b, 2012) apply to steady-state configurations of ice sheets that require that all environmental conditions (in this study expressed via the net surface accumulation/ablation and submarine melting) do not change with time. Neither the present-day ice sheets nor paleo ice sheets have ever attained steady states, because the Earth's climate changes in time on a variety of temporal scales. Hence, the concepts of marine ice-sheet (in)stability cannot be used to describe their behavior. This does not imply, however, that analyzes of steady states are useless. The purpose of such analyzes are different from those focused on the ice-sheets’ response to time-variant environmental conditions (past or future). By ignoring the variability of external conditions, the steady-state analyzes focus on the internal dynamics of the ice sheets and provide insights into their dynamics and conceptual understanding of complex behavior controlled by poorly known subglacial processes, which effects are encapsulated in sliding laws. Additionally, these steady-state analyzes have been a driving force for the development of ice-sheet numerical models, providing benchmarks and establishing the requirements for numerical algorithms and the model resolutions.

Understanding of the behavior of the realistic ice sheets requires analyzes that take into account changing environmental conditions. This requires shifting the focus from steady-state considerations to the temporal evolution of the ice sheets and their characteristics such as the grounding line migration. If its rate is small or close to zero that means that the ice sheets are close or at equilibrium with their environmental conditions; nevertheless, they are not in steady states, because if these conditions change the rate of the grounding line migration changes in response to them (Eqn. (19)). Currently, there is neither quantitative nor conceptual understanding of such dynamic equilibria, their characteristic timescales, controls, etc. Developments of such approaches need to become the focus of future studies.

7. Conclusions

We have developed a generalized framework for the analysis of the steady-state behavior and dynamic evolution of marine ice sheets, and applied it to a configuration with a power-law lateral shear and Zoet-Iversion and Weertman sliding laws. When driven with stochastic variability in submarine melt rates, the grounding line evolution strongly depends on the specific modeling choices (e.g., the size of the modeled domain, the shape of the bed, the conditions at the calving front). Our results illustrate that classifications in stable or unstable grounding lines are not useful concepts in understanding the transient behavior of ice sheets. Although all our simulations started from stable steady-state configurations, their response to temporal variability, even of a stochastic nature with no trends, cannot be described in terms of ‘stable’ or ‘unstable’ at any given time. Understanding the dynamic behavior of marine ice sheets requires a new conceptual approach that takes into account the ice-sheet environmental conditions and their variability on a wide range of temporal scales.