2. Model description

Following many previous studies (e.g., Schoof and others, Reference Schoof, Devis and Popa2017; Pegler, Reference Pegler2018; Haseloff and Sergienko, Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022), 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, Reference MacAyeal1989)

(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 (Reference Haseloff and Sergienko2018) 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, Reference Haseloff and Sergienko2018). 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, Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022). 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, Reference Holt2006), and laterally confined glaciers like Greenland marine outlet glaciers, Pine Island Glacier in the Amundsen Sea Embayment (e.g., Wingham and others, Reference Wingham, Wallis and Shepherd2009), ice streams flowing into the Filchner-Ronne and Amery ice shelves (e.g., King and others, Reference King, Woodward and Smith2007), and many others. Such a stress regime has been considered by Weertman (Reference Weertman1974); Wilchinsky (Reference Wilchinsky2009); Schoof (Reference Schoof2007a, Reference Schoof2007b); Schoof (Reference Schoof2012); Tsai and others (Reference Tsai, Stewart and Thompson2015); Pegler (Reference Pegler2016, Reference Pegler2018); Schoof and others (Reference Schoof, Devis and Popa2017); Haseloff and Sergienko (Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022); Sergienko and Wingham (Reference Sergienko and Wingham2022); Sergienko (Reference Sergienko2022b). 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, Reference Rose1979; Bindschadler and Vornberger, Reference Bindschadler and Vornberger1998); it has been considered by Sergienko and Wingham (Reference Sergienko and Wingham2019). 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 (Reference Haseloff and Sergienko2018) 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 (Reference Schoof2007a, Reference Schoof2007b); and Sergienko and Wingham (Reference Sergienko and Wingham2022) 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 (Reference Haseloff and Sergienko2022) and Sergienko and Wingham (Reference Sergienko and Wingham2022), and in contrast to studies based on the boundary-layer theory (Schoof, Reference Schoof2007a, Reference Schoof2012; Tsai and others, Reference Tsai, Stewart and Thompson2015) 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 (Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022) for Weertman sliding law in the case of buttressed marine ice sheets, by Sergienko and Wingham (Reference Sergienko and Wingham2022) in the case of unconfined marine ice sheets, and in case of negligible bed slopes and accumulation rate by Schoof (Reference Schoof2007a, Reference Schoof2007b).

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 (Reference Schoof2012), 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, Reference Sergienko and Wingham2019). 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, Reference Bassis and Ultee2019), the existing formulations for τ ^w depend on ice velocity (e.g., Raymond, Reference Raymond1996). 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, Reference Thomas1977), 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, Reference Fyke, Sergienko, Löfveström, Price and Lenaerts2018) 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, Reference Sergienko2022b).

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, Reference Raymond1996; Schoof and others, Reference Schoof, Devis and Popa2017; Haseloff and Sergienko, Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022; Sergienko, Reference Sergienko2022a)

(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 (Reference Raymond1996).) 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 (Reference Zoet and Iverson2020) 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 (Reference Zoet and Iverson2020), we obtain α ≈ 1.1 × 10⁻¹¹ Pa⁻¹ m s⁻¹.

As Haseloff and Sergienko (Reference Haseloff and Sergienko2022) 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 (Reference Zoet and Iverson2020) 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 (Reference Joughin, Smith and Schoof2019), 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, Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022). 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, Reference Asay-Davis2016), 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, Reference Powell2009). 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 (Reference Haseloff and Sergienko2022) 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 (Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022) 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 (Reference Haseloff and Sergienko2022)(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 (Reference Sergienko and Wingham2019, Reference Sergienko and Wingham2022) and Sergienko (Reference Sergienko2022b) 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 (Reference Haseloff and Sergienko2018) 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 (Reference Tsai, Stewart and Thompson2015) 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ć, Reference Milanković1941), the Earth's climate has numerous modes of internal variability. Steig and others (Reference Steig, Ding, Battisti and Jenkins2012) 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 (Reference Jenkins2018); Christie and others (Reference Christie, Steig, Gourmelen, Tett and Bingham2023) have attributed the grounding line migration of these ice shelves to decadal variability. The results of numerous observational (e.g., Shepherd and others, Reference Shepherd, Wingham and Rignot2004) and modeling (e.g., Cornford and others, Reference Cornford2015; Seroussi and others, Reference Seroussi2017) 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, Reference Christian, Robel and Catania2022; Felikson and others, Reference Felikson, Nowicki, Nias, Morlighem and Seroussi2022) and realistic (Hoffman and others, Reference Hoffman, Asay-Davis, Price, Fyke and Perego2019; Robel and others, Reference Robel, Seroussi and Roe2019) 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, Reference Jenkins2018; Christie and others, Reference Christie, Steig, Gourmelen, Tett and Bingham2023) and modeling (Latif and others, Reference Latif, Martin and Park2013).

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, Reference Donea, Huerta, Ponthot and Rodríguez-Ferran2017). 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.