Effects of calving and submarine melting on steady states and stability of buttressed marine ice sheets

Abstract

Mass loss from ice shelves is a strong control on grounding-line dynamics. Here we investigate how calving and submarine melt parameterizations affect steady-state grounding-line positions and their stability. Our results indicate that different calving laws with the same melt parameterization result in more diverse steady-state ice-sheet configurations than different melt parameterizations with the same calving law. We show that the backstress at the grounding line depends on the integrated ice-shelf mass flux. Consequently, ice shelves are most sensitive to high melt rates in the vicinity of their grounding lines. For the same shelf-averaged melt rates, different melt parameterizations can lead to very different ice-shelf configurations and grounding-line positions. If the melt rate depends on the slope of the ice-shelf draft, then the positive feedback between increased melting and steepening of the slope can lead to singular melt rates at the ice-shelf front, producing an apparent lower limit of the shelf front thickness as the ice thickness vanishes over a small boundary layer. Our results illustrate that the evolution of marine ice sheets is highly dependent on ice-shelf mass loss mechanisms, and that existing parameterizations can lead to a wide range of modelled grounding-line behaviours.

1. Introduction

Increased mass loss from marine-terminating outlet glaciers in Greenland and Antarctica is usually attributed to warming ocean waters (Holland and others, 2008; Joughin and others, 2010; Pritchard and others, 2012; Straneo and Heimbach, 2013). Key for linking ocean warming and glacier retreat is the buttressing effect of ice shelves, which transmit changes in ocean forcing to the grounded ice. Modelling studies investigating the dynamics of buttressed marine ice sheets have shown that laterally confined ice shelves can suppress a possible instability of marine ice sheets on retrograde slopes (Hughes, 1973; Mercer, 1978; Dupont and Alley, 2005; Schoof, 2007a; Goldberg and others, 2009; Gudmundsson and others, 2012; Gudmundsson, 2013; Pegler, 2018b). Consequently, a reduction of the ice-shelf buttressing caused, for example, by submarine melting might trigger instability (Mercer, 1978), and the grounding-line retreat simulated in large-scale ice-sheet models under projected ocean warming conditions (e.g. Joughin and others, 2014; Robel and others, 2019; Rosier and others, 2021) is often interpreted as a manifestation of this instability.

The marine ice-sheet instability, first proposed by Weertman (1974), is the consequence of a monotonic dependence of the ice flux at the grounding line on ice thickness, which is in turn linked to the depth of the sea floor through the floatation condition (Weertman, 1974; Thomas and Bentley, 1978; Schoof, 2007b). If such a monotonically increasing relationship between flux and ice thickness exists, then retreat of a grounding line onto a deeper sea floor leads to an increase in the ice flux, ice-sheet thinning and continued retreat of the grounding line until a position is reached where the ice flux decreases with further retreat or the ice sheet becomes grounded. Several studies have challenged the simplicity of this picture by demonstrating that a wide range of effects, including ice-shelf buttressing, alters this instability behaviour (e.g. Gomez and others, 2010; Pegler and Worster, 2012; Robel and others, 2014, 2016; Pegler, 2016; Brondex and others, 2017; Sergienko and Wingham, 2019, 2022). However, a systematic investigation of the effects of ice-shelf mass loss processes on the stability of buttressed marine ice sheets is still outstanding.

1.1. Ice-shelf mass loss

Ice shelves lose their mass by two processes: calving of icebergs at the ice-shelf front and submarine melting. Calving is a manifestation of the brittle nature of ice, and modelling individual calving events is beyond the scope of present-day ice-sheet models (e.g. Åström, 2013; Bassis and Jacobs, 2013). Thus, practical approaches to model calving in large-scale ice-sheet models include: (1) fixing the calving front position (e.g., Gudmundsson and others, 2012; Arthern and Williams, 2017), (2) keeping the length of the ice shelf fixed (Gagliardini and others, 2010) and (3) the application of strain-rate-based criteria (e.g. Vieli and others, 2000, 2001; Benn and others, 2007a; Nick and others, 2010; Levermann and others, 2012; Morlighem and others, 2016). In flowline models that do not resolve the flow transverse to the main flow direction, strain rates at the calving front can be linked to the ice thickness there (Vieli and others, 2000, 2001; Benn and others, 2007a; Nick and others, 2010); these calving criteria have some observational backing (Van der Veen, 1996). However, there is no clear understanding of how the choice of the calving law affects the grounding-line response to external forcing. Here, we aim to address this question by comparing the steady-state grounding-line positions obtained with different calving criteria.

Submarine melting results from the interaction of the ocean circulation with sub-shelf cavities. While fully coupled ice–ocean models are being developed (Goldberg and others, 2012, 2018; De Rydt and Gudmundsson, 2016; Seroussi, 2017; Jordan and others, 2018), the coupling of ice and ocean dynamics remains a major challenge. Such calculations require significant computational resources and detailed knowledge of the ocean conditions (e.g. Goldberg and others, 209) and atmospheric forcing, while at the same time providing complex outputs which make identification of the fundamental feedbacks difficult.

For use in large-scale ice-sheet models, a wide range of submarine melt parameterizations have been proposed. These include an explicit spatial dependence (i.e. a dependence on fixed coordinates (x,  y), e.g. Gagliardini and others, 2010; Aschwanden and others, 2019), a dependence on the ice-shelf draft (e.g. Joughin and others, 2010, 2014; Pollard and DeConto, 2012; Favier and others, 2014), or reduced representations of plume dynamics.

Subglacial plume models aim to capture the leading order effects of ice-shelf–ocean interactions caused by buoyant subglacial discharge at the grounding line or buoyant melt water at the ice–ocean interface (e.g. Jenkins, 1991, 2011; Sergienko and others, 2013; Slater and others, 2017; Hewitt, 2020). As this water rises along the ice-shelf face, it entrains sea water. This sea water provides heat to the plume, leading to further melting and increased buoyancy. Crucially, entrainment rates are often assumed to depend on the slope of the ice-shelf base, with larger slopes leading to higher melt rates. The existence of a positive feedback between melting and ice-shelf shape can theoretically explain the formation of undercut ice shelves, and this kind of feedback is explicitly incorporated in models that simulate the evolution of a buoyant plume at the ice-shelf base. Parameterizations of these dynamics typically include a dependence on the slope of the ice-shelf base (e.g. Lazeroms and others, 2018, 2019).

1.2. Study outline

With a diverse range of parameterizations of ice-shelf mass loss it becomes important to develop a basic understanding of how they affect the flow of the grounded ice sheet individually and in different combinations. To build a qualitative understanding, we use simplified representations of these parameterizations to investigate the interaction between submarine melting, calving, ice shelf shape and grounding-line position. From the ice–ocean interaction perspective, our study can be considered complementary to Slater and others (2017), who consider a full plume model and a drastically simplified ice model (by prescribing a fixed ice velocity). In comparison, we use a flowline ice model accounting for both mass and momentum balance, but we use a range of simplified melt parameterizations.

From an ice-dynamics perspective, our study builds on existing work by Schoof and others (2017), Haseloff and Sergienko (2018) and Pegler (2018a, 2018b) who consider different theoretical aspects of buttressed marine ice-sheet dynamics, but generally neglect melting. To do this, we follow the same approach as Schoof (2007b, 2011, 2012), Haseloff and Sergienko (2018) and Sergienko and Wingham (2019, 2022): we derive an expression for the flux at the grounding line, and use it to perform a linear stability analysis. Our results show that the stability criterion derived by Schoof (2012) only holds for sufficiently smooth beds and certain calving laws. Under these conditions, the relationship between the flux gradient and the accumulation rate determines the stability of steady states.

Like these other models we sacrifice model complexity (and thus applicability) for model tractability by using a flowline model with parameterized side drag. Flowline models cannot account for across-flow variations in geometric properties, flow or melt (e.g. Gudmundsson, 2003; Sergienko, 2012, 2013; Reese and others, 2018a; Zhang and others, 2020). Nevertheless, flowline models have been helpful in developing an understanding of buttressed grounding-line dynamics (Dupont and Alley, 2005; Nick and others, 2010; Hindmarsh, 2012; Jamieson and others, 2012; Pegler and others, 2013; Robel and others, 2014, 2016; Pegler, 2016, 2018a, 2018b; Schoof and others, 2017; Haseloff and Sergienko, 2018, Sergienko, 2022).

The paper is organized as follows: we present the model in Section 2. To build an understanding of the qualitative differences between different melt rate parameterizations, we briefly discuss the response of unconfined ice shelves to these in Section 3, before turning to buttressed marine ice sheets in Section 4. In Section 4.1 we show that different calving laws can render the same steady-state grounding-line position stable or unstable. In Section 4.2 we demonstrate that the buttressed grounding-line flux depends at leading order on the integrated ice flux and illustrate that this results in melting focused closer to the grounding line reducing buttressing more than the same amount of melting applied further away from it. Finally, in Section 4.3 we compare the steady-state grounding-line positions obtained with different melt rate and calving parameterizations, and show that the variability in grounding-line positions due to different calving laws is greater than the variability due to different melt rate parameterizations.

2. The model

We consider a laterally confined marine ice sheet of constant width W flowing in the positive x-direction from an ice divide at x = 0 to the calving front at x = x _c (Fig. 1). At the grounding line x = x _g the ice-shelf forms.

Fig. 1. Cross section of the model.

2.1. Full model equations

We model the velocity u along the centre line of the ice stream/shelf and the ice thickness h with a widely used laterally averaged version of a plug flow or ‘shallow stream’ model (e.g. MacAyeal, 1989; Dupont and Alley, 2005; Nick and others, 2009; Hindmarsh, 2012; Pegler, 2016; Haseloff and Sergienko, 2018)

(1a)$$\eqalign{& 2 { {\rm d}\over {\rm d} {x}}\left(A^{-1/n}h \left\vert {{\rm d} {u}\over {\rm d} {x}} \right\vert ^{1/n-1} {{\rm d} {u}\over {\rm d} {x}}\right)- {C_{\rm w} A^{-1/n} \over W^{1/n + 1}} h \vert u\vert ^{1/n-1}u \cr& - C\vert u\vert ^{m-1}u -\rho_{\rm i} g h{{\rm d} {( h + b) }\over {\rm d} {x}} = 0,\; \quad{\rm if }\, 0< x\leq x_{\rm g}}$$

(1b)$$\eqalign{& 2 {{\rm d} {}\over {\rm d} {x}}\left(A^{-1/n}h \left\vert {{\rm d} {u}\over {\rm d} {x}} \right\vert ^{1/n-1} {{\rm d} {u}\over {\rm d} {x}}\right)- {C_{\rm w} A^{-1/n} \over W^{1/n + 1}} h \vert u\vert ^{1/n-1}u \cr& - \rho_{\rm i} g \left(1-{\rho_{\rm i}\over \rho_{\rm w}} \right)h{{\rm d} {h}\over {\rm d} {x}} = 0,\; \quad {\rm if }\, x_{\rm g}< x\leq x_{\rm c}}$$

with the first term in (1a) the divergence of the longitudinal shear stress, the second term the lateral shear stress, the third term the basal shear stress and the last term the driving stress. The parameter A is the rate factor of the ice viscosity, n is a rheology parameter, C _w is a lateral shear stress parameter, ρ _i is the density of ice, ρ _w is the density of water, g is the gravitational constant, m and C are sliding coefficients, and

(2)$$b( x) = 100 - 2184.8\times\left({x\over 155\, {\rm km}}\right)^{2} + 1031.72\times\left({x\over 155\, {\rm km}}\right)^{4} - 151.72\times\left({x\over 155\, {\rm km}}\right)^{6}$$

is the elevation of the bed, for which we use a scaled version of the MISMIP and the MISMIP+ bed geometries (Cornford and others, 2020). Parameter values used in this study are listed in Table 1. The equation for mass conservation is

(3)$${{\rm d} {( uh) }\over {\rm d} {x}} = \left\{\matrix{ \dot a \hfill & {\rm if }\, 0\leq x\leq x_{\rm g},\; \hfill \cr \dot m \hfill & {\rm if }\, x_{\rm g}< x< x_{\rm c} \hfill }\right.$$

with $\dot a$ the net accumulation rate on the ice sheet and $\dot m$ the net mass balance term of the ice shelf, combining both surface and basal mass balance contributions (positive: accumulation/freeze-on). We are strictly interested in the effect of changing the ice-shelf mass balance and therefore assume a constant accumulation rate $\dot a$ on the grounded ice sheet but different models for the submarine melt rate, given in Eqn (5). To simplify the derivation of asymptotic solutions, we neglect accumulation on the ice shelf; including a constant accumulation term does not change the results qualitatively.

Table 1. Model parameters with their value, where applicable

The boundary conditions are

(4a)$${{\rm d} {( h + b) }\over {\rm d} {x}} = u = 0 \quad {\rm at }\; x = 0$$

(4b)$$h = h_{\rm g} = -{\rho_{\rm i}\over \rho_{\rm w}}b( x_{\rm g}) \quad {\rm at }\; x = x_{\rm g}$$

(4c)$$2 A^{-1/n} h\left\vert {{\rm d} {u}\over {\rm d} {x}} \right\vert ^{{1}/{n}-1}{{\rm d} {u}\over {\rm d} {x}} = {1\over 2}\rho_{\rm i} g\left(1-{\rho_{\rm i}\over \rho_{\rm w}}\right)h^{2} \;{\rm at }\; x = x_{\rm g} + L_{\rm s}.$$

Additionally, we require the velocity u, the stress (and consequently the velocity gradient ${\rm d}u/{\rm d}x$), and the ice thickness h to be continuous at the grounding line x _g. To close this model, we still require a condition to describe the length of the ice shelf L _s and a model for the sub-shelf melt rate $\dot m$ which we describe next.

2.2. Ice-shelf mass loss

2.2.1 Melting

As direct coupling of ice-sheet models with ocean-circulation models is computationally expensive (e.g. Jordan and others, 2018), a multitude of submarine melt rate parameterizations exist which attempt to capture the relevant processes with either location-dependent parameterizations (e.g. Gagliardini and others, 2010), ice-thickness-dependent parameterizations (e.g. Joughin and others, 2010; Favier and others, 2014) or parameterizations which include a dependence on the ice-shelf slope (e.g. Little and others, 2012; Lazeroms and others, 2018, 2019; Favier and others, 2019).

The focus of this study is to gain insights into the qualitative differences in the grounding-line dynamics resulting from the use of different melt rate parameterizations. We therefore consider a limited subset of idealized parameterizations representative of typical choices:

(5)$$\dot m = \left\{\matrix{ & \hskip-10pt f( x) \hfill \cr & \hskip-12pt \gamma_2 h^{2} \hfill \cr & \hskip-12pt \gamma_3 ( h_{\rm g}-h) \left\vert {{\rm d} {h}\over {\rm d} {x}}\right\vert {1\over \sqrt{1 + \epsilon^{2} ( {{\rm d} {h}}/{{\rm d} {x}}) ^{2}}},\; \quad \epsilon\ll 1. }\right.$$

The first parameterization assumes a dependence on the downstream coordinate x only, that is, we assume melting to be independent of the ice-shelf thickness h or slope dh/dx. The second parameterization (5) ₂ depends on the ice-shelf draft (depth of the ice-shelf base below sea level) and is similar to two parameterizations that have been used in studies of Pine Island (Joughin and others, 2010; Favier and others, 2014). The important assumption underlying this parameterization is that melting increases with depth; consequently, melt rates are largest directly at the grounding line.

The slope-dependent parameterization (5) ₃ is a strong simplification of the melt rate derived by Lazeroms and others (2018), who derive a general parameterization of the melt rates obtained with a plume model (Jenkins, 1991). In their parameterization, the melt rate is proportional to sinα where α is the slope of the ice draft, and an 11th-order polynomial in the distance from the grounding line. We include the former dependence through an explicit slope-dependent term dh/dx, and the latter through the dependence on (h _g − h). The term $1/\sqrt {1 + \epsilon ^{2} ( {{\rm d} {h}}/{{\rm d} {x}}) ^{2}}$ with $0< {\epsilon}\ll 1$ is included as regularization of the melt rate in the limit of dh/dx → ∞.

2.2.2. Calving

Unless all ice-shelf mass is removed through melting, solution of the ice-flow model (1)–(4) additionally requires a condition that determines the position of the calving front x _c. This is generally done through a calving law, and similarly to melt rate parameterizations, a large range of options exist (e.g. Van der Veen, 1998; Dupont and Alley, 2005; Benn and others, 2007b; Goldberg and others, 2009; Nick and others, 2010; Gagliardini and others, 2010; Jamieson and others, 2012; Gudmundsson and others, 2012).

The most numerically convenient of these choices is the assumption of a fixed calving front position x _c = const. (e.g. Gudmundsson and others, 2012; Arthern and Williams, 2017) as this does not require explicitly tracking the movement of the calving front. An alternative class of calving laws is based on stress or ice thickness criteria, effectively prescribing a certain ice thickness at the calving front as extensional stress and ice thickness there are linked through (4c) (e.g. Vieli and others, 2000; Benn and others, 2007a; Nick and others, 2010; Choi and others, 2021). A third calving choice we will consider assumes a fixed ice-shelf length, i.e. assuming x _c = x _g + L ₀ with L ₀ = const. While this choice would be computationally difficult to implement in laterally extended numerical ice-sheet models, it is particularly suitable for our investigation of the roles of melting and calving on grounding-line dynamics through a simplified flowline model as it does not introduce dynamical feedbacks through a varying ice-shelf length.

To summarize, the position of the calving front is here set through:

(6)$$x_{\rm c} =\left\{ \matrix{x_g + L_0 \hfill& {\rm where}\ L_0={\rm const.} \hfill\cr {\rm const.}\hfill & \cr x^{{\dagger}} \hfill&{\rm where}\ h(x^{{\dagger}}) = h_{c} = {\rm const. }\hfill }\right.$$

with the first option describing a constant ice-shelf length L ₀, the second option describing a constant calving front position x _c, and the third option describing a constant ice thickness at the calving front h _c. Note however that in the case of strong melting it is possible for the ice thickness to go to zero before this location in which case the ice-shelf length is set by the location where h goes to zero. The length of the ice shelf is thus given by

(7)$$L_{\rm s} = \min( x_{\rm c}-x_{\rm g},\; x^{\ast }-x_{\rm g}) \quad {\rm with}\quad h( x^{\ast }) = 0.$$

2.3 Reduced model equations

In the grounded part of the marine ice-sheet longitudinal stress gradients can be neglected in the momentum balance and the above problem can be described by a reduced model (Schoof, 2007b; Haseloff and Sergienko, 2018; Sergienko and Wingham, 2022):

(8a)$$- {C_{\rm w} A^{-1/n} \over W^{1/n + 1}} h \vert u\vert ^{1/n-1}u - C\vert u\vert ^{m-1}u -\rho_{\rm i} g h{{\rm d} {( h + b) }\over {\rm d} {x}} = 0,\; \quad{\rm if }\, 0< x\leq x_{\rm g}$$

(8b)$${{\rm d} {q}\over {\rm d} {x}} = {{\rm d} {( uh) }\over {\rm d} {x}} = \dot a,\; \quad{\rm if }\, 0< x\leq x_{\rm g}.$$

This is essentially a ‘shallow ice’ model (Fowler and Larson, 1978; Morland and Johnson, 1980), and the boundary conditions of this system are:

(9a)$${{\rm d} {( h + b) }\over {\rm d} {x}} = u = 0 \quad\quad\quad{\rm at }\;\; x = 0$$

(9b)$$h = h_{\rm g} = -{\rho_{\rm i}\over \rho_{\rm w}}b( x_{\rm g}) \quad\quad\quad{\rm at }\;\; x = x_{\rm g}$$

(9c)$$uh = q_{\rm g}\quad\quad\quad\quad {\rm at }\;\; x = x_{\rm g}.$$

with (9c) effectively replacing the stress boundary condition (4c) (see, e.g. Schoof, 2007b, for details). It is straightforward to see that for the steady-state problems considered here, the flux at grounding line has to match the integrated mass balance (8b):

(10)$$q_{\rm g} = \int_0^{x_{\rm g}}\dot a \;{\rm d} x.$$

Determining an expression for the flux at the grounding line therefore closes the reduced system.

2.3.1 Grounding-line flux

If the longitudinal stress gradient is small at the grounding line, an expression for the ice flux can be derived from the continuity condition of the stress at the grounding line (see e.g. Hindmarsh, 2012; Sergienko and Wingham, 2022), i.e.

(11a)$$\tau_{{\rm stream}} = \tau_{{\rm shelf}} \quad {\rm at }\quad x = x_{\rm g},\; $$

where

(11b)$$\tau_{{\rm stream}} = {2\over A^{1/n}} h\left\vert {{\rm d} {u}\over {\rm d} {x}} \right\vert ^{{1}/{n}-1}{{\rm d} {u}\over {\rm d} {x}}$$

and

(11c)$$\tau_{{\rm shelf}} = \tau_0\times \Theta\quad {\rm with }\quad \tau_0 = {1\over 2} \, \rho_{\rm i} g \left(1-{\rho_{\rm i}\over \rho_{\rm w}}\right)h_{\rm g}^{2}.$$

Θ is the ratio between the backstress at the grounding line and the unbuttressed backstress; we discuss its calculation in detail in Section 2.3.2.

From (8a), and reasonably assuming u > 0 (i.e. ice flows towards the calving front) and du/dx > 0

(12)$${{\rm d} {h}\over {\rm d} {x}} = -\left[{C_{\rm w}\over A^{{1}/{n}}W^{1/n + 1}\rho_{\rm i} g} u^{{1}/{n}} + {C\over \rho_{\rm i} g} {u^{m}\over h} + {{\rm d} {b}\over {\rm d} {x}}\right].$$

Substituting this expression into (11a) and rearranging terms leads to an implicit expression for the ice flux at the grounding line

(13)$$\eqalign{& {{\rm d} {q}\over {\rm d} {x}}h^{1/n + m + 2} + {C_{\rm w}\over A^{{1}/{n}}W^{1/n + 1}\rho_{\rm i} g} q^{1/n + 1}h^{m + 1} \cr & \quad+ {C\over \rho_{\rm i} g}q^{m + 1}h^{1/n} + q h^{1/n + m + 1}{{\rm d} {b}\over {\rm d} {x}} \cr & \quad = \left[{A^{1/n}\over 4}\rho_{\rm i} g\left(1-{\rho_{\rm i}\over \rho_{\rm w}}\right)\right]^{n} h^{1/n + m + 3+n}\Theta^{n}}$$

where dq/dx is determined by (8b) at x = x _g; for steady-states ${{\rm d} {q}}/{{\rm d} {x}} = \dot a$. Note that we only used the momentum balance to determine (13), therefore it is valid for both steady-state and time-evolving problems.

The implicit flux expression (13) contains the flux expression proposed by Schoof (2007a, b)

(14)$$q_{\rm g} = q_0 \times \Theta^{{n}/{( m + 1) }} \quad {\rm with } \quad q_0 = \left({A ( \rho_{\rm i} g) ^{n + 1} ( 1-\rho_{\rm i}/\rho_{\rm w}) ^{n}\over 4^{n} C} \right)^{{1}/{( m + 1) }} h_{\rm g}^{{( 3 + m + n) }/{( m + 1) }}$$

as a limiting case. Writing (13) in similar form, one can obtain

(15)$$q_g= \frac{q_0}{\left(1+\Gamma\right)^{1/(m+1)}} \times \Theta^{{n}/{(m+1)}} \quad\text{ with }\quad \Gamma= q_x\frac{\rho_i g}{C} h_g^{m+2}q_g^{-m-1} + \frac{\rho_i g}{C}\frac{\text{d}b}{\text{d}x}h_g^{m+1}q_g^{-m} + \frac{C_w A^{-1/n}}{C W^{1/n+1}}h_g^{m+1-1/n}q_g^{1/n-m} \quad $$

with the non-unity parameters in the denominator accounting for the effects of flux gradients (first term in Γ), bed gradients (second term) and lateral shear (third term) on the grounded part, upstream of the grounding line.

We have specifically chosen our model configuration and parameters to satisfy the original assumptions underlying Schoof's asymptotic theory, so that the two expressions (14) and (13) yield almost indistinguishable results in most of the domain. We point out that at the far downstream side of the bed, past ≈ 300 km, the solutions of (13) cease to exist as the bed slope becomes too large. Mathematically, this is the result of the bracketed term in the denominator of (15) becoming negative in this regime. Consequently, we do not consider configurations with grounding-line positions past this point.

By adopting the grounding-line flux expression (13), we neglect a wide range of dynamic behaviour that is caused by varying basal boundary conditions on the grounded part of the ice sheet (e.g. Tsai and others, 2015; Brondex and others, 2017; Sergienko and Wingham, 2019). We also restrict ourselves to the simplest possible geometry of a laterally confined parallel-sided outlet glacier, as the model considered here is only suitable under such conditions (e.g. Reese and others, 2018b). While this is a clear limitation of our model, neglecting the influence of the basal ice-sheet properties in particular facilitates the interpretation of our results considerably.

2.3.2. Backstress at the grounding line

In either of the flux expressions above Θ must be determined from the ice-shelf equations. In some cases this is possible analytically, through construction of asymptotic solutions of the ice-shelf equations (1b) and (3) ₂ with (10) and (4c), see Appendix B. There we find that the steady-state grounding-line stress for a spatially variable melt rate $\dot m = f( x)$ is:

(16)$${\Theta = 1 {-} \left[{( 4^{n} C_{\rm w}) ^{1/( 1 + n) } [ q( x_{\rm c}) ] ^{1/n} + {n + 1 \over n}{C_{\rm w} \over W}\int_{x_{\rm g}}^{x_{\rm c}}[ q( x') ] ^{1/n}{\rm d} x'\over \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) A^{1/n}W^{1/n}h_{\rm g}^{1/n + 1}} \right]^{2n/( n + 1) }}$$

with q(x) the steady-state flux on the ice shelf $q( x) = q_{\rm g} + \int _{x_{\rm g}}^{x} f( x')\; {\rm d} x'.$ We discuss this expression in detail in Section 4.2.

As it is not always possible to construct analytical solutions, we also calculate Θ numerically with Matlab ODE solvers. In this case we use a Newton method to determine the backstress τ _shelf that simultaneously satisfies the ice-shelf momentum balance (1b) and mass balance condition (3), with boundary conditions (10) and (4c). We refer to this approach as semi-analytical.

Finally, verification of our solutions with full numerical solutions of (1)–(7) is obtained with Comsol and an adapted version of the time-dependent grounding-line code described in Schoof (2007a) and Robel and others (2014), which has been extended with an ice-shelf solver to account for buttressing; these solutions are referred to as numerical solutions. To summarize, we use three different ways to determine grounding-line position:

1. (1) analytical solution of (10) and (13) with suitably constructed Θ (e.g. (16), (18) and (25))

2. (2) semi-analytical solution of (10) and (13) with numerically determined Θ

3. (3) numerical solution of the full system of equations ((1)–(7)). In this case, the time-dependent mass balance equation

   $${\partial {h}\over \partial {t}} + {\partial {( uh) }\over \partial {x}} = \dot a \quad {\rm if }\, 0\leq x \leq x_{\rm g}$$
   is solved instead of (3) ₁, and the model is run into steady state.

We compare the different approaches in Appendix A and Section 4.2, where we show good agreement between them.

3. Unconfined ice-shelf profiles for different melt rate parameterizations

Before discussing the effect of melting on buttressed marine ice sheets below, we start by considering unconfined ice shelves to build intuition for the interpretation of the buttressed solutions. For an unconfined ice shelf, Eqns (1b) with (3) ₂ can be solved without recourse to the equations for the grounded ice sheet, if we prescribe the ice flux q and ice thickness h at the grounding line:

(17)$$q = q_{\rm g},\; \quad h = h_{\rm g}\; {\rm at }\; x = x_{\rm g}.$$

As outlined in the previous section, we additionally need to prescribe the location of the calving front x _c. For an unconfined ice shelf whose grounding-line position is insensitive to the ice-shelf processes, assuming a fixed calving front position or a constant ice-shelf length are identical choices, and we prescribe x _c to be fixed at x _c = x _g + L ₀. Note however that in the case of strong melting it is possible for the ice thickness to go to zero upstream of this location. In this case x _c is set by the location where h goes to zero. The ice-shelf length L _s is thus given by (7).

Figure 2 shows typical ice-shelf profiles and velocities obtained for the three different melt rate parameterizations given in (5). Each of the parameters γ₁, γ₂ and γ₃ represents a different physical quantity. In order to be able to compare different melt rates, we therefore use the average melt rate on the shelf as reference: profiles plotted in the same colours in panels a₁–c₃ are for the same average melt rates, which are indicated in panels a₄–c₄.

Fig. 2. Examples for unconfined ice-shelf profiles (row 1) and velocities (row 2) for different melt rate distributions (row 3). Profiles in the same colour have the same average melt rates $\bar {\dot m} = L_{\rm s}^{-1}\int _{x_{\rm g}}^{x_{\rm g} + L_{\rm s}}\dot m \;{\rm d} x$ as indicated in panels a₄–c₄. Note that the y-axis of panel b₄ is logarithmic, as the ice thickness remains non-zero for all γ₂ ≤ 0.

We might expect melt to reduce the ice thickness until it goes to zero at the ice-shelf front, at which point the ice shelf starts to shorten. The solutions for a constant melt rate (column a) indeed exhibit this behaviour, with the ice-shelf shortening once γ₁ < −q _g/L ₀. However, for the second melt rate ($\dot m = \gamma _2 h^{2}$, column b) the ice thickness at the calving front can never reach zero, and the ice shelf thins most close to the grounding line and subsequently maintains an almost constant ice thickness towards the calving front. For this parameterization the average melt rate cannot exceed − q _g/L ₀, and as this limit is approached for γ₂ → −∞, the ice thickness at the calving front asymptotically approaches zero (Fig. 2b₄).

The third melt rate parameterization exhibits yet another behaviour (Fig. 2c): the ice shelf appears to retreat before the ice thickness at the shelf front reaches zero and there appears to be a lower bound on the ice thickness at the ice-shelf front (Fig. 2c₄). Closer examination of these solutions (Appendix B.1) shows that this is only an apparent lower bound: for γ₃ less than a critical γ_c, the melt rate at the ice-shelf front goes to − 1/ε, so that the ice thickness goes to zero in a small boundary layer of length O(ε) (see 8d–f), with ε the regularization parameter in (5).

These examples illustrate that using different melt rate parameterizations leads to very different ice-shelf profiles. With a constant melt rate, the ice shelf becomes almost triangular in the limit of strong melting. Conversely, the ice-thickness dependent melt rate leads to a shape that quickly thins away from the grounding line, as melting is strongest there. When we account for the positive feedback between submarine melt and slope, the ice shelf only thins up to a critical value given in Eqn (B7), beyond which the ice shelf retreats with an seemingly non-zero ice thickness at the ice-shelf front.

Melting affects not only the ice thickness, but also the velocity, as for example illustrated by the velocity profiles in Figure 2a₂: as melting increases the velocity decreases towards the ice-shelf front compared to solutions with weaker melting. This is a consequence of melting reducing the ice flux, not only the ice thickness in steady state. We will see in Section 4.2 that this effect is even more pronounced for buttressed ice shelves experiencing melt.

4. Buttressed marine ice sheets

4.1. The role of the calving law

Before considering the effect of different melt rate parameterizations on grounding-line positions, we focus on the effects of the calving law. We consider steady-state configurations of a 40 km wide marine ice sheet on an overdeepened bed in the absence of submarine melting (Fig. 3). We use the three different calving criteria described in (6), with the calving front position (x _c ≈ 380 km), ice-shelf length (L ₀ ≈ 155 km), and ice thickness at the calving front (h _c ≈ 416 m) chosen in such a way that all three calving laws produce the same steady-state configuration on the retrograde part of the bed (panels a₂, b₁ and c₂, respectively). All other parameter values are listed in Table 1.

Fig. 3. Extent and stability of a 40 km wide outlet glacier on a bed given by (2) and constant accumulation on the grounded ice sheet. Column a: results for fixed ice-shelf length of 155 km; column b: results for fixed calving front position at x _c = 380 km, column c: results for fixed ice thickness h _c = 415 m at the calving front. Panels a₄–c₄ show the solutions for the buttressed ice flux q _g at the grounding line and the integrated accumulation (dashed line). At the intersection of the flux q _g with the integrated accumulation steady-state grounding-line positions are possible, which are shown in panels a₁–a₃, b₁, and c₁–c₂.

4.1.1. Steady-state configurations

In steady state, two conditions are simultaneously satisfied at the grounding line: the ice flux matches the accumulation integrated over the grounded part of the ice sheet (10), and the ice thickness is at floatation (9b). To determine grounding-line positions analytically, we use the asymptotic solutions of Θ for constant melt rates derived in Haseloff and Sergienko (2018):

(18)$$\eqalign{ \Theta &= 1-{1\over h_{\rm g}^{2}} \left({q_{\rm g}\over \left[\rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) \right]^{n} A W}\right)^{2/( n + 1) } \cr &\quad \times \left[\left(4^{n} C_{\rm w}\right)^{1/( n + 1) } + {1 + n\over n}C_{\rm w} {L_{\rm s}\over W}\right]^{2n/( n + 1) }.}$$

For a fixed ice thickness at the calving front (6) ₃, we additionally need a relation between the ice thickness at the calving front and the ice-shelf length, which we approximate with (Haseloff and Sergienko, 2018)

(19a)$$\eqalign{& h_{\rm c} = \left[ h_{{\rm c}, {\rm unconf}}^{{( n + 1) ^{2}}/{n}}\, {\rm erfc}\left({C_{\rm w}\over 2}{L_{\rm s}^{1 + 1/n}\over W^{1 + 1/n}}\right)\right. \cr & \quad\quad\left. +\; h_{{\rm c}, {\rm conf}}^{{( n + 1) ^{2}}/{n}}{\rm erf}\left({C_{\rm w}\over 2}{L_{\rm s}^{1 + 1/n}\over W^{1 + 1/n}}\right)\right] ^{{n}/{( n + 1) ^{2}}}.}$$

h _c,unconf is the ice thickness at the calving front of unconfined ice shelves

(19b)$$h_{{\rm c}, {\rm unconf}} = {2^{2n/( n + 1) }h_{\rm g} q_{\rm g}^{1/( n + 1) }\over \left[2^{2n}q_{\rm g} + ( n + 1) [ \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) ] ^{n} h_{\rm g}^{1 + n} A L_{\rm s}\right]^{1/( n + 1) }}$$

and h _c,conf is the ice thickness at the calving front of strongly buttressed ice shelves

(19c)$$h_{{\rm c}, {\rm conf}} = \left[{ 4^{n} C_{\rm w} q_{\rm g}^{1/n + 1}\over \left[\rho_{\rm i} g( 1-\rho_{\rm i}/\rho_{\rm w}) \right]^{1 + n}A^{1/n + 1} W^{1/n + 1}} \right]^{n/( n + 1) ^{2}}.$$

For illustrative purposes, we solve (8)–(9) numerically with Matlab ODE solvers to obtain the ice-sheet profiles, and (1b) and (3) ₂ with (10) and (4c) to obtain the ice-shelf profiles.

In our example, the grounding-line flux obtained with a fixed ice-shelf length mirrors the shape of the bed and the analytical and semi-analytical solutions predict three possible steady states (Fig. 3a₁–a₃), two on the downwards sloping sections of the bed, and one on the upwards sloping part of the bed. Conversely, for a fixed calving front position x _c = const., we find only one steady-state grounding-line position, located on the retrograde section of the bed (Fig. 3b₁), as the flux appears to be monotonically increasing over the entire domain (Fig. 3b₄). Finally, for a fixed ice thickness at the calving front, the ice flux at the grounding line appears independent of the bed shape over most of the domain, with the upstream branch of the flux closely following the shape of the unbuttressed flux (not shown), and the downstream branch of the ice flux being constant. This leads to the existence of two equilibrium grounding-line positions (Fig. 3c₁–c₂), one close to the unbuttressed grounding-line solution with a very short ice shelf (Fig. 3c₁) and one on the retrograde part of the bed (Fig. 3c₂). Note that the analytical flux for a constant thickness at the calving front peaks at x ≈ 100 km. This is due to the approximation of the calving-front ice thickness with (19a), which is an ad-hoc approximation of the ice thickness at the calving front and only formally correct in the two limits of no and strong buttressing (see Haseloff and Sergienko, 2018).

For all three calving laws, we do not compute the flux at the grounding line q _g past $x\,\raise2pt{{\mathop{>}\limits_{\raise1pt{\approx} }}}\, 300\,{\rm km}$. As mentioned above, the slope of the bed db/dx has a large magnitude there, resulting in the bracketed term in the denominator of (13) being negative; such steep bedrock slopes also violate the assumptions of the shallow shelf approximation, which our model is based on (MacAyeal, 1989).

To provide more than one or two points of comparison per calving law, we compare the semi-analytical and analytical results with numerical results for a range of different widths in the appendix (see Fig. 7); there we generally find good agreement between semi-analytical solutions and the numerical solutions.

4.1.2. Stability of steady states

To determine the stability of these steady states, we perform a linear stability analysis described in Appendix C. It shows that the stability of these steady states cannot be determined from comparison of the flux gradient with the accumulation rate alone. Instead, the stability condition depends on the chosen calving law, melt rate parameterization (if we additionally take melting into account) and strength of the bed slope.

In the absence of melting, and assuming that the calving law can be expressed through a calving front position x _c(x _g,  h _g,  q _g) which is a function of grounding-line position x _g, ice thickness at the grounding line h _g, and ice flux at the grounding line q _g, the general stability condition (C20) can be written in terms of the flux gradient dq _g/dx _g and the accumulation rate $\dot a$ (see (C26))

(20)$$\eqalign{& {\rm if }\, \left(A_1 + B {q_{\rm g}^{1/n-1}\over n} \left(( 4^{n} C_{\rm w}) ^{{1}/{( n + 1) }} W + C_{\rm w} ( x_{\rm c}-x_{\rm g}) \right) \right.\cr & \left.\quad+ B C_{\rm w} q_{\rm g}^{1/n} {\partial {x_{\rm c}}\over \partial {q_{\rm g}}}\right)> 0,\; \, \zeta> 0,\; \left(h_x-{{\rm d} {h_{\rm g}}\over {\rm d} {x_{\rm g}}} \right)< 0 \cr & {\rm then\; for}\quad \left\{\matrix{ \displaystyle\dot a < {{\rm d} {q_{\rm g}}\over {\rm d} {x_{\rm g}}} \quad {\rm stable} \cr \displaystyle\dot a > {{\rm d} {q_{\rm g}}\over {\rm d} {x_{\rm g}}} \quad {\rm unstable.} }\right.}$$

We point out that the flux gradient depends not only on the bed slope (through dh _g/dx _g), but also on the accumulation gradient (${{\rm d} {\dot a}}/{{\rm d} {x_{\rm g}}}$), bed curvature (${{\rm d} {{}^{2} b}}/{{\rm d} {x_{\rm g}^{2}}}$) and buttressing; it is (C23)

(21)$${{{\rm d} {q_{\rm g}}\over {\rm d} {x_{\rm g}}} \!=\! { ( B_1-A_2 - B C_{\rm w} q_{\rm g}^{1/n}{\partial {x_{\rm c}} \over \partial {h_{\rm g}}}) {{\rm d} {h_{\rm g}} \over {\rm d} {x_{\rm g}}} - A_3{{\rm d} {\dot a} \over {\rm d} {x_{\rm g}}} - A_4{{\rm d} {{}^{2} b} \over {\rm d} {x_{\rm g}^{2}}} + B C_{\rm w} q_{\rm g}^{1/n}\left(1 - {\partial {x_{\rm c}} \over \partial {x_{\rm g}}} \right)\over A_1 + q_{\rm g}^{1/n-1} B \left[{1 \over n}( 4^{n} C_{\rm w}) ^{1/( n + 1) }W + {1 \over n}C_{\rm w} ( x_{\rm c}-x_{\rm g}) + q_{\rm g} {\partial {x_{\rm c}} \over \partial {q_{\rm g}}}\right]}.}$$

The parameters A ₁–A ₄ are defined in (C12), the parameters B and B ₁ are defined in (C14) and ζ is given by (C11).

This condition is significantly more complex than the stability conditions derived in Schoof (2012) and Sergienko and Wingham (2022), as it takes into account not only accumulation and bed gradients, but also the lateral confinement of the grounded ice sheet and buttressing by the ice shelf. The requirement $A_1 + B {q_{\rm g}^{1/n-1}\over n} ( ( 4^{n} C_{\rm w}) ^{{1}/{( n + 1) }} W\; +\; C_{\rm w} ( x_{\rm c}-x_{\rm g}) )\; +\; B C_{\rm w} q_{\rm g}^{1/n}{\partial {x_{\rm c}}}/{\partial {q_{\rm g}}} > 0$ is a necessary condition to make statements about stability with a linear stability analysis, i.e. without using a full numerical model. ζ > 0 is required for steady states to exist; otherwise the denominator in (13) is complex. The condition h _x − dh _g/dx _g is the requirement that the ice upstream of the grounding line remains grounded (Schoof, 2012). Provided these conditions are satisfied, the difference between the accumulation rate and the flux gradient can indeed be used to determine the stability of a steady state, in the same manner as in the case of Schoof's flux formula (Schoof, 2007b, 2012).

The terms $B{q_{\rm g}^{1/n-1}\over n} ( ( 4^{n} C_{\rm w}) ^{{1}/{( n + 1) }} W + C_{\rm w} ( x_{\rm c}-x_{\rm g}) )$ and $BC_{\rm w} q_{\rm g}^{1/n}$ are positive by definition, and we have explicitly chosen our geometry to ensure A ₁ > 0. For a constant calving front position (x _c = constant) or a constant ice-shelf length (x _c = x _g + L _s, L _s = constant), the partial derivative ∂x _c/∂q _g vanishes, and we can use a visible comparison between the flux and integrated accumulation to determine the stability of steady states in Figures 3a, b. For the solutions with a constant ice-shelf length, we can identify the two steady states on the prograde bed (a₁ and a₃) as stable grounding-line positions, while the steady-state grounding-line position on the retrograde part of the bed (a₂) is unstable. However, if the same steady state is obtained with a fixed calving front position, then the steady-state position on the retrograde bed (b₁) becomes stable. The position and stability of these steady states is also confirmed by numerical calculations, see Appendix A.

The flux gradient can be determined from (21), but this expression is complex and does not allow for a straightforward identification of the leading-order processes. To understand the role of the calving law better, we therefore use the limit of strong buttressing (formally when L _s ≫ W) to obtain a simplified analytic expression for the grounding-line flux from setting Θ = 0 in (18) (see Haseloff and Sergienko, 2018, for details):

(22a)$$q_{{\rm g}, L_{\rm s}} = {[ \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) ] ^{n} A W^{n + 1}\over ( 4^{n} C_{\rm w}) ^{1/( n + 1) }W + {1 + n \over n}C_{\rm w} L_{\rm s}}h_{\rm g}^{n + 1} \quad {\rm if }\, L_{\rm s} = {\rm const.},\; \, W\ll L_{\rm s}.$$

Equation (22a) is plotted as a dashed yellow line in Figures 3a₄. As it qualitatively agrees with the other solutions, we can use (22a) to better understand the behaviour shown in Figure 3a₁–a₃: (22a) predicts that all spatial variability in the ice flux is due to changes in the ice thickness at the grounding line h _g = −ρ _w/ρ _i b(x _g). This is qualitative similar to the ice-flux expression for unconfined ice sheets by Schoof (2007b) (14), and therefore explains why this calving law reproduces the instability of grounding lines on upwards sloping beds: ${{\rm d} {q_{{\rm g}, L_{\rm s}}}}/{{\rm d} {x_{\rm g}}}$ is negative for all areas where the bed slopes upwards (db/dx _g > 0), so that the grounding line is always unstable on retrograde beds for a constant, positive accumulation rate $\dot a$.

For a constant calving front position, the general stability condition (20) predicts that the sole steady state on the retrograde section of the bed is stable. The stability of this grounding line is also confirmed by the numerical solutions in Appendix A. As for the constant ice-shelf length, an explicit flux expression can be obtained in the limit of strong buttressing:

(22b)$$q_{{\rm g}, x_{\rm c}} = {[ \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) ] ^{n} A W^{n + 1}\over ( 4^{n} C_{\rm w}) ^{1/( n + 1) }W + {1 + n \over n}C_{\rm w} ( x_{\rm c}-x_{\rm g}) }h_{\rm g}^{n + 1} \quad {\rm if }\, x_{\rm c} = {\rm const.},\; \, W\ll L_{\rm s}.$$

Now, even as the depth of the seafloor decreases with increasing x _g (db/dx _g > 0) on upwards sloping parts of the bed, it is possible for the flux to increase due to the second term in the denominator, which decreases with increasing x _g. Physically, this corresponds to the fact that the length of the ice shelf and hence the amount of buttressing decreases as the grounding line advances. Less buttressing leads to an increase in the ice flux, and this increase in ice flux can exceed the reduction in ice flux due to a smaller ice thickness at the grounding line, as shown in Figure 3b₄.

When the ice thickness at the calving front is prescribed, Figure 3 shows that we have two branches, the upstream branch where the ice thickness at the calving front is set by (19b) and one where it is set by (19c). For the upstream branch, we can use (19b) to determine ∂x _c/∂q _g as

$${\partial {x_{\rm c}}\over \partial {q_{\rm g}}} = {4^{n} \left[1- ( h_{\rm c}/h_{\rm g}) ^{n + 1}\right]\over ( n + 1) [ \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) ] h_{\rm g}^{n + 1}A}> 0 \quad {\rm if }\, h_{\rm c} = h_{{\rm c}, {\rm unconf.}}$$

ensuring $A_1 + B {q_{\rm g}^{1/n-1}\over n} ( ( 4^{n} C_{\rm w}) ^{{1}/{( n + 1) }} W + C_{\rm w} ( x_{\rm c} -x_{\rm g}) )+ B C_{\rm w} q_{\rm g}^{1/n} $ ${\partial {x_{\rm c}}\over \partial {q_{\rm g}}} > 0) $. Consequently, we can determine from equation (20) that these configurations are stable.

The stability condition (20) does not extend to the case where the ice thickness at the calving front is set by (19c), as this ice thickness directly sets the flux at the grounding line (which one can obtain from rearranging (19c)):

(22c)$$q_{{\rm g}, h_{\rm c}} = {[ \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) ] ^{n}\over ( 4^{n} C_{\rm w}) ^{n/( n + 1) }} A W h_{\rm c}^{n + 1}.$$

For a simplified flux condition of this form the Sturm–Liouville problem considered in Appendix C reduces to the regular Sturm–Liouville problem considered in Schoof (2012), and his stability conditions remain valid. As ${{\rm d} {q_{{\rm g}, h_{\rm c}}}}/{{\rm d} {x_{\rm g}}} = 0$, this branch is always unstable for positive accumulation rates.

These examples illustrate that different calving laws produce very different grounding-line positions for buttressed marine ice sheets and that the choice of calving law can change the stability of the grounding line from stable to unstable and vice versa. This suggests that the grounding-line response of buttressed marine ice sheets to different melt forcings will depend on the calving law as well.

4.2. The role of the melt location

We have seen in Section 3 that different melt distributions result in very distinct ice-shelf profiles. In this section and the next, we aim to understand the response of buttressed ice sheets to different melt parameterizations. Numerical studies (e.g., Gagliardini and others, 2010) suggest that melting closer to the grounding-line affects the grounding-line response more than melting further away from it. We can understand this behaviour by determining the steady-state grounding-line stress for a spatially variable melt rate $\dot m = f( x)$ (see Appendix B.2):

(23a)$$\!\!\!\!\!\!\eqalign{\tau_{{\rm shelf}}\over \tau_{0} & = \Theta \cr &= 1\!-\! \left[{( 4^{n} C_{\rm w}) ^{1/( 1 + n) } [ q( x_{\rm c}) ] ^{1/n} + {n + 1 \over n}{C_{\rm w} \over W}\int_{x_{\rm g}}^{x_{\rm c}}[ q( x') ] ^{1/n}{\rm d} x'\over \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) A^{1/n}W^{1/n}h_{\rm g}^{1/n + 1}} \right]^{2n/( n + 1) }.}$$

τ ₀ is the grounding-line stress of an unconfined ice shelf, given by (11c) ₃, and q(x) is the steady-state flux on the ice shelf:

(23b)$$q( x) = q_{\rm g} + \int_{x_{\rm g}}^{x} f( x')\; {\rm d} x'.$$

Equation (23a) shows that the steady-state grounding-line stress depends on the integrated ice flux and the calving flux, that is, the ice flux at the ice-shelf front. Both of these quantities decrease with increasing melt, albeit in different ways. The calving flux is insensitive to the location of melting, and only depends on the total amount of melting experienced throughout the ice shelf, $q( x_{\rm c}) = q_{\rm g} + \int _{x_{\rm g}}^{x_{\rm c}}f( x')\; {\rm d} x'$ with the second term the integrated melt rate on the ice shelf. Conversely, the integrated flux depends on the melt location: the integral $\int _{x_{\rm g}}^{x_{\rm c}}q^{1/n}\; {\rm d} x$ is smaller if q is reduced closer to x _g for the same amount of reduction. A consequence of (23a) is therefore that melting closer to the grounding line reduces the backstress more than the same amount of melting farther away from the grounding line.

To illustrate the effect of the melt location on the grounding-line position, we consider a δ-function melt perturbation:

(24)$$\dot m = -\mu q_{\rm g} \times \delta ( x_{\rm m}-x) ,\; \quad 0< \mu< 1,\; \quad x_{\rm g}< x_{\rm m}< x_{\rm g} + L_{\rm s}$$

i.e. we melt only at a single point on the ice shelf and remove a fraction μ of the grounding-line flux there (note that δ(x − x _m) has units m⁻¹). We start our calculations from the stable steady state shown in Figure 3a₁ and keep the ice shelf length constant (L ₀ = 155 km) to avoid the non-linear effects introduced through more complicated calving laws (Section 4.1).

Figures 4a, b show steady-state profiles for μ = 0.5 and different relative melt locations

$$x_{\rm r} = ( x_{\rm m}-x_{\rm g}) /L_0,\;$$

that is, for different locations of melting relative to the total ice-shelf length; at the grounding line x _r = 0. Clearly, in this example moving the location of melting closer to the grounding line leads to steady-state grounding-line positions further inland. Note also the pronounced drop in ice-shelf velocities right downstream of the melt perturbation.

Fig. 4. Influence of relative melt location x _r = (x _m − x _g)/L ₀ on grounding-line position. x _r = 0 corresponds to melting at the grounding line, x _r = 1 for melting at the calving front. Panels show ice-sheet/shelf profiles (a), velocities (b) and grounding-line positions (c) for melt applied at different relative melt positions x _r (24). Note that melting closer to the grounding line leads to more grounding-line retreat than melting further farther away from it. Analytic solution obtained with (25).

For this example, the stress (23a) with melt rate (24) can be written as a function of relative melt-rate position x _r:

(25)$${\eqalign{& \Theta = 1-{q_{\rm g}^{2/( n + 1) }\over h_{\rm g}^{2}} \cr & \times \left[{\left(4^{n} C_{\rm w}\right)^{1/( n + 1) }( 1\!\!-\!\!\mu) ^{1/n} + {1 + n \over n}C_{\rm w} {L_{\rm s} \over W}[ x_{\rm r} + ( 1\!-\!\mu) ^{1/n}( 1\!\!-\!\!x_{\rm r}) ] \over \rho_{\rm i} g ( 1-\rho_{\rm i}/\rho_{\rm w}) A^{1/n} W^{1/n}}\right]^{2n/( n + 1) }.}$$

with (1 − μ) the fraction of the grounding-line ice flux not removed by melt. For constant μ, Eqn (25) describes the reduction in buttressing as a monotonically increasing function of the relative melt location x _r, with the change in τ _shelf being largest for melting directly at the grounding line (x _r = 0) and reducing towards the ice-shelf front (x _r = 1). Equation (25) together with (14) provides an analytic (albeit implicit) expression of the grounding-line flux, and these solutions are shown as black lines in Figure 4, together with semi-analytical (red lines) and numerical solutions (dashed blue lines).

While this example of strong localized melting is instructive in understanding some of the factors controlling grounding-line position, it is hardly realistic. We will therefore turn to more realistic (but still highly idealized) melt rates next.

4.3. Combination of melt rate and calving parameterizations

In this section, we explore how steady-state grounding-line positions depend on different combinations of the calving laws considered in Section 4.1 and the melt parameterizations considered in Section 3. We use the same parameters that we used for the examples shown in Figure 3 (listed in Table 1) and the different melt parameterizations introduced in Eqn (5). To calculate grounding-line positions, we use the semi-analytical method described above, which we have shown in Section 4.2 and Appendix A to agree well with the numerical solutions.

For all combinations, we start with the same melt-free configurations shown in Figure 3 (γ₁ = γ₂ = γ₃ = 0). We then vary the melt rate parameters γ₁, γ₂ or γ₃, respectively, while keeping either the calving front position fixed (green lines in Fig. 5), the ice-shelf length fixed (purple lines) or the ice thickness at the calving front fixed (magenta lines). Columns a–c in Figure 5 correspond to the different melt rate parameterizations considered, and the first three panels show the corresponding steady-state ice-sheet/shelf profiles for each calving law. For completeness, we also include solutions for positive values of γ_i in panels a₄–c₄ (but not in the profiles), even though these correspond to freeze-on, and the models are not applicable in this parameter range.

Fig. 5. Steady-state solutions for three different calving laws (indicated by colour) and three different sub-shelf melt parameterizations (indicated by column titles). Ice-sheet/shelf profiles (panels a₁–c₃) are shown in 50 km intervals, but only for negative values of γ_i, as the melt rate parameterizations are not applicable for freeze-on. Values of the melt rate parameters γ_i at corresponding steady-state grounding-line positions shown in panels a₄–c₄. Note that γ_i has different units for different melt rate parameterizations, therefore numerical values are not directly comparable between different parameterizations. For all melt rates more negative values correspond to more melting.

The ice-shelf profiles at one fixed grounding-line position can differ substantially both between different melt rate parameterizations (compare for example panels a₁–c₁) and different calving laws (compare for example panels b₁–b₃). At least the former should not come as a surprise, as this mirrors the behaviour seen in Section 3. There are, however, also some common characteristics. For example the solutions with a fixed ice thickness at the calving front (panels a₃–c₃) are virtually indistinguishable for different melt rate parameterizations. We can also note that the solutions for a fixed ice-shelf length and for a fixed calving front with a constant melt rate are identical in the cases of strong melting when the length of the ice shelf is set by the position where all ice is removed through melting (panels a₁ and a₂). The same is true for the slope-dependent melt rate when the melt rate parameter γ₃ exceeds the critical value γ_c (panels c₁ and c₂).

γ₁, γ₂ and γ₃ represent different physical quantities, therefore they are not directly comparable. Instead, we compare the average ice-shelf melt rate, Figure 6, with the colour indicating the calving law and the line type indicating the melt rate parameterization. Perhaps the most striking property of Figure 6 is that solutions with the same calving law are more similar to each other than solutions with the same melt rate parameterization, and that their shape roughly mirrors the shape of the grounding-line flux in the absence of melt, compare Figures 3a₄–c₄. Nevertheless, there is a large spread in average melt rates associated with different steady-state grounding-line positions. For example, the average melt rate $\bar {\dot m}$ associated with a grounding line close to the unbuttressed grounding line at x ≈ 80 km ranges from $\bar {m}\approx -6$ m a⁻¹ to $\bar {m}\approx -1$ m a⁻¹ (ignoring positive values, as these are only included for completeness).

Fig. 6. Average melt rate $\bar {\dot m} = L_{\rm s}^{-1}\int _{x_{\rm g}}^{x_{\rm g} + L_{\rm s}}\dot m \;{\rm d} x$ vs grounding-line position for solutions shown in Figure 5.

Figure 6 illustrates that the grounding-line positions obtained with the thickness-based calving law (magenta lines) follow qualitatively different relationships between average melt rate and steady-state grounding-line position than the other two calving laws. The grounding-line positions are largely insensitive to the melt distribution and the average melt rates required to change the steady-state grounding-line position are substantially lower than the melt rates required for the other calving laws. This behaviour is linked to the grounding-line flux for this calving law, see Figure 3, which is essentially bi-modal: either mostly following the flux expression for the unbuttressed case, or constant for different grounding-line positions.

The other two calving laws overlap when subjected to the strong melt forcing with constant ($\dot m = \gamma _1$, solid green and purple lines) or slope-dependent melting ($\dot m = \gamma _3( h_{\rm g}-h) \vert {\rm d}h/{\rm d} x\vert$, dash-dotted green and purple lines). This is not surprising, as in these cases the ice-shelf length is set by the location where the ice thickness goes to zero, rather than the calving law (compare profiles in a₁–a₂ and c₁–c₂, respectively).

The average melt rate associated with a particular steady-state grounding-line position is usually the highest for the plume-like parameterization, as it predicts strong frontal melting for large values of γ₃. Our results in 4.2 showed that melting at the ice-shelf front has less impact on the backstress at the grounding line, even though it leads to high average melt rates. Conversely, the average melt rates obtained with the thickness-based parameterization $\dot m = \gamma _2 h^{2}$ are lower than for the other melt rates, as in this case melting is strongest at the grounding line, where it has maximum impact because of the integral dependence of the buttressing parameter on melting.

For a constant melt rate $\dot m = \gamma _1$, we can use the results of the linear stability analysis in Appendix C to determine the stability of the steady states shown here, which generally reproduces the same stability behaviour that we have seen in Section 4.1. However, these results cannot be extended to the other two parameterizations, as the stability condition (C20) is specific to the flux parameterization (13) with (23a), the solution for Θ for $\dot m = f( x)$. If the expressions for Θ were known for the other two melt rate parameterizations, it would be possible to extend the stability analysis using the same methodology as presented in Appendix C.

5. Discussion

In this study, we have considered how different ice-shelf mass-loss processes affect the steady-state configurations and stability of confined marine ice sheets. To this end, we have derived an expression for the ice flux at the grounding line, which extends the results of earlier studies by Schoof (2007b), Haseloff and Sergienko (2018) and Sergienko and Wingham (2022), and accounts for the lateral confinement of the grounded ice sheet and ice-shelf buttressing. Using this expression, we have compared steady-state grounding-line positions obtained with different melt rate and calving parameterizations. By means of a linear stability analysis, we have derived stability conditions for a specific form of the buttressing parameter Θ.

5.1. The role of calving

Our results suggest that the variability in grounding-line positions due to different calving laws is greater than the variability due to different melt rate parameterizations. We find that using different calving parameterizations can change the response of grounding lines to forcing qualitatively – from a stable to an unstable configuration. In contrast, using different melt parameterizations only changed this response quantitatively by changing the average melt rate associated with a particular grounding-line position.

This difference is due to the different ways these parameterizations are formulated: the calving laws we introduced change the length of the ice shelf as a function of grounding-line position (due to their dependencies on x _g, h _g and q _g), compare for example the simplified flux expressions (22a)–(22c). In contrast, the melt parameterizations differ in their spatial distribution of melt rates, but do not alter the length of the ice shelf for most of the parameters we considered. Only when melting becomes strong enough to determine the ice-shelf length, do the differences between different calving laws vanish (e.g. see the grounding-line positions for x _c = const. and L _s = const. in the limit of strong melting, Figs 5a₄, c₄).

The grounding-line flux depends on the specifics of the calving law and melt parameterization, therefore a generalized stability criterion cannot be derived. Instead, stability conditions have to be evaluated on a case-by-case basis for specific formulations. Our linear stability analysis provides a template for this process and gives a stability criterion for melt rates of the form $\dot m = f( x)$. Provided the bed is sufficiently smooth and the steady-state calving front position does not retreat with increasing grounding-line flux, the stability condition can be formulated in terms of the gradient of the grounding-line flux dq _g/dx _g and the accumulation rate, in line with Schoof (2012). Application of these results illustrates that stable grounding-line positions on upwards sloping beds obtained with a constant calving front position are indeed due to the increase of ice-shelf length and hence buttressing with grounding-line retreat (as suggested by Schoof and others, 2017).

However, linear stability analyses that can provide stability conditions are only possible if the flux expression is explicitly known. Even the simple melt rate parameterizations considered here do not always provide a closed-form expression of the backstress at the grounding line. Therefore generalized inferences about the (in)stability of one particular glacier are impossible without detailed numerical studies.

5.2. The role of melting

We find that the buttressed grounding-line flux depends at leading order on the integrated ice-shelf flux or the double integral of the melt-rate distribution (Eqn (23)). Numerical studies suggest that melting closer to the grounding-line affects grounding-line dynamics more than melting away from it (e.g. Gagliardini and others, 2010), and our results provide an analytical confirmation and explanation of these numerical results: the closer to the grounding-line melting occurs, the stronger is the reduction in the integrated ice-shelf flux.

We explore the effects of using different melt parameterizations, which can lead to a large spread in ice-shelf profiles and grounding-line positions for the same average melt rate. In particular, one of the melt parameterizations we considered includes a positive feedback between melting and ice-shelf slope (see (5) ₃), whereas increased melting leads to steeper slopes of the ice-shelf base, which in turn increases melting there. In our model this positive feedback leads to a threshold γ_c of the melt rate parameter γ₃, which regulates the strength of melting. If γ₃ < γ_c the above feedback leads to singular melt rates at the ice-shelf front, akin to strong frontal melting. This gives the appearance of the ice only thinning up to an apparent minimum ice-shelf thickness at the ice-shelf front; past this minimum thickness, further melting appears to lead to ice-shelf retreat, rather than further thinning. We emphasize that this is only an apparent minimum ice-shelf thickness, as the ice thickness in fact goes to zero over a short boundary layer whose extent depends on a regularization parameter.

The existence of this positive feedback has been recognized by other authors (e.g. Jenkins, 1991; Little and others, 2012; Sergienko and others, 2013; Slater and others, 2017). The main difference between our study and these authors is the direct coupling with an ice-sheet/shelf flow model. The coupling with the ice-shelf flow model is essential for the feedback to operate fully, as both ice thickness and flow change in response to melting (see also Sergienko and others, 2013). The latter point is best illustrated by our results for strong localized melting (Section 4.2), which illustrates that the response to melting is not only thinning. In some of our examples the velocity response to melting is much more pronounced than the ice thickness response. This is because for the steady states considered here, melting reduces the ice flux, which is the product of ice thickness and velocity, and it is not immediately obvious which of these will show a larger response. In other words, without modelling ice-shelf flow and its mass balance simultaneously, one might expect the ice shelf to simply thin in response to melting, but clearly this is not always the case. This also reiterates that one might expect very different long-term and short-term responses to melting, and that further studies of this time-dependent behaviour are necessary to understand these dynamics.

Another difference to these earlier studies is that they use the full plume model, which means the singularity in (5) ₃, which we regularized through inclusion of an artificial regularization parameter ε, is removed and the positive feedback is automatically bounded. Nevertheless, the existence of this positive feedback poses questions about the importance melting plays in settings which allow for the formation of strong submarine plumes. If strong melting leads to severe undercutting, then this could trigger calving of the undercut ice-shelf parts. This could appear as ice-thickness or stress-threshold calving similar to observed calving behaviour of Columbia Glacier, Alaska, USA (Van der Veen, 1996). That said, observations and modelling studies also suggest that the positive feedback might be limited by other factors, for instance, the inclusion of lateral variations in melt (e.g. Gladish and others, 2012; Sergienko, 2013; Langley and others, 2014).

5.3. Model limitations

Systematic comparison of different melt parameterizations and calving laws in time-dependent numerical models is challenging, not least because the melt rate parameters used in different parameterizations are not directly comparable and the numerical effort to accurately calculate grounding-line positions remains significant. Our approach allows us to compare nine different combinations of melt rates and calving laws, but the price for this computational efficiency is a simplified model that neglects several important aspects.

We consider steady states only, and therefore cannot make any statements about the transient evolution of buttressed marine ice sheets. There are several timescales to consider for the ice shelf: in a time-dependent model the immediate effect of changing the melt rate is a change in the ice thickness h (as e.g. considered in Reese and others, 2018a; Zhang and others, 2020), rather than the ice flux q as the depth-integrated mass balance is ${\partial {h}}/{\partial {t}} + \nabla \cdot {\bf q} = \dot m.$ It is thus conceivable that transient changes in the grounding-line stress differ substantially from the steady-state response. As such, our steady-state results can only give an indication of the direction of change we anticipate. Moreover, in time-dependent models, calving rates are often applied, i.e. instead of prescribing the ice-shelf length directly, the rate of change in calving front position $\dot {x_{\rm c}}$ is prescribed, which introduces an additional timescale. Depending on this timescale, the relative importance of calving and melting processes might change if the timescale associated with calving is long in comparison with the timescale over which the ice-shelf geometry adjusts to melting.

We have strictly focused on geometric settings that allow us to apply the flux parameterizations developed in Schoof (2007a), Haseloff and Sergienko (2018) and Sergienko and Wingham (2022). This requires a constant width. We have also neglected the effects of varying basal boundary conditions. Relaxing any of these assumptions is known to alter the dynamics of buttressed marine ice sheets (Gomez and others, 2010; Robel and others, 2014, 2016; Brondex and others, 2017; Åkesson and others, 2018; Reese and others, 2018b; Sergienko and Wingham, 2019, 2022).

6. Conclusions

Calving and melt parameterizations strongly influence the grounding-line positions obtained for marine ice sheets with laterally confined ice shelves. With everything else held constant, different calving laws result in qualitatively different ice-sheet configurations and stability properties. Conversely, for a given calving parameterization, ice-sheet steady-state configurations and their stability are less sensitive to different melt parameterizations. If the melt rate depends on the ice-shelf slope, melt rates can become singular, leading to an apparent lower limit of the ice thickness at the ice-shelf front. Consequently, model results of marine ice-sheet evolution strongly depend on the way ice-shelf mass loss is parameterized.