1. Introduction
Increased mass loss from marine-terminating outlet glaciers in Greenland and Antarctica is usually attributed to warming ocean waters (Holland and others, Reference Holland, Thomas, De Young, Ribergaard and Lyberth2008; Joughin and others, Reference Joughin, Smith and Holland2010; Pritchard and others, Reference Pritchard2012; Straneo and Heimbach, Reference Straneo and Heimbach2013). 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, Reference Hughes1973; Mercer, Reference Mercer1978; Dupont and Alley, Reference Dupont and Alley2005; Schoof, Reference Schoof2007a; Goldberg and others, Reference Goldberg, Holland and Schoof2009; Gudmundsson and others, Reference Gudmundsson, Krug, Durand, Favier and Gagliardini2012; Gudmundsson, Reference Gudmundsson2013; Pegler, Reference Pegler2018b). Consequently, a reduction of the ice-shelf buttressing caused, for example, by submarine melting might trigger instability (Mercer, Reference Mercer1978), and the grounding-line retreat simulated in large-scale ice-sheet models under projected ocean warming conditions (e.g. Joughin and others, Reference Joughin, Smith and Medley2014; Robel and others, Reference Robel, Seroussi and Roe2019; Rosier and others, Reference Rosier2021) is often interpreted as a manifestation of this instability.
The marine ice-sheet instability, first proposed by Weertman (Reference Weertman1974), 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, Reference Weertman1974; Thomas and Bentley, Reference Thomas and Bentley1978; Schoof, Reference Schoof2007b). 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, Reference Gomez, Mitrovica, Huybers and Clark2010; Pegler and Worster, Reference Pegler and Worster2012; Robel and others, Reference Robel, Schoof and Tziperman2014, Reference Robel, Schoof and Tziperman2016; Pegler, Reference Pegler2016; Brondex and others, Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017; Sergienko and Wingham, Reference Sergienko and Wingham2019, Reference Sergienko and Wingham2022). 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, Reference Åström2013; Bassis and Jacobs, Reference Bassis and Jacobs2013). Thus, practical approaches to model calving in large-scale ice-sheet models include: (1) fixing the calving front position (e.g., Gudmundsson and others, Reference Gudmundsson, Krug, Durand, Favier and Gagliardini2012; Arthern and Williams, Reference Arthern and Williams2017), (2) keeping the length of the ice shelf fixed (Gagliardini and others, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010) and (3) the application of strain-rate-based criteria (e.g. Vieli and others, Reference Vieli, Funk and Blatter2000, Reference Vieli, Funk and Blatter2001; Benn and others, Reference Benn, Hulton and Mottram2007a; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Levermann and others, Reference Levermann2012; Morlighem and others, Reference Morlighem2016). 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, Reference Vieli, Funk and Blatter2000, Reference Vieli, Funk and Blatter2001; Benn and others, Reference Benn, Hulton and Mottram2007a; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010); these calving criteria have some observational backing (Van der Veen, Reference Van der Veen1996). 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, Reference Goldberg2012, Reference Goldberg2018; De Rydt and Gudmundsson, Reference De Rydt and Gudmundsson2016; Seroussi, Reference Seroussi2017; Jordan and others, Reference Jordan2018), 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, Reference Goldberg, Gourmelen, Kimura, Millan and Snow209) 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, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010; Aschwanden and others, Reference Aschwanden2019), a dependence on the ice-shelf draft (e.g. Joughin and others, Reference Joughin, Smith and Holland2010, Reference Joughin, Smith and Medley2014; Pollard and DeConto, Reference Pollard and DeConto2012; Favier and others, Reference Favier2014), 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, Reference Jenkins1991, Reference Jenkins2011; Sergienko and others, Reference Sergienko, Goldberg and Little2013; Slater and others, Reference Slater, Nienow, Goldberg, Cowton and Sole2017; Hewitt, Reference Hewitt2020). 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, Reference Lazeroms, Jenkins, Gudmundsson and van de Wal2018, Reference Lazeroms, Jenkins, Rienstra and van de Wal2019).
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 (Reference Slater, Nienow, Goldberg, Cowton and Sole2017), 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 (Reference Schoof, Davis and Popa2017), Haseloff and Sergienko (Reference Haseloff and Sergienko2018) and Pegler (Reference Pegler2018a, Reference Pegler2018b) 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 (Reference Schoof2007b, Reference Schoof2011, Reference Schoof2012), Haseloff and Sergienko (Reference Haseloff and Sergienko2018) and Sergienko and Wingham (Reference Sergienko and Wingham2019, Reference Sergienko and Wingham2022): 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 (Reference Schoof2012) 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, Reference Gudmundsson2003; Sergienko, Reference Sergienko2012, Reference Sergienko2013; Reese and others, Reference Reese, Gudmundsson, Levermann and Winkelmann2018a; Zhang and others, Reference Zhang, Price, Hoffman, Perego and Asay-Davis2020). Nevertheless, flowline models have been helpful in developing an understanding of buttressed grounding-line dynamics (Dupont and Alley, Reference Dupont and Alley2005; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Hindmarsh, Reference Hindmarsh2012; Jamieson and others, Reference Jamieson2012; Pegler and others, Reference Pegler, Kowal, Hasenclever and Worster2013; Robel and others, Reference Robel, Schoof and Tziperman2014, Reference Robel, Schoof and Tziperman2016; Pegler, Reference Pegler2016, Reference Pegler2018a, Reference Pegler2018b; Schoof and others, Reference Schoof, Davis and Popa2017; Haseloff and Sergienko, Reference Haseloff and Sergienko2018, Sergienko, Reference Sergienko2022).
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.
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, Reference MacAyeal1989; Dupont and Alley, Reference Dupont and Alley2005; Nick and others, Reference Nick, Vieli, Howat and Joughin2009; Hindmarsh, Reference Hindmarsh2012; Pegler, Reference Pegler2016; Haseloff and Sergienko, Reference Haseloff and Sergienko2018)
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
is the elevation of the bed, for which we use a scaled version of the MISMIP and the MISMIP+ bed geometries (Cornford and others, Reference Cornford2020). Parameter values used in this study are listed in Table 1. The equation for mass conservation is
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.
Model parameters with their value, where applicable
The boundary conditions are
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, Reference Jordan2018), 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, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010), ice-thickness-dependent parameterizations (e.g. Joughin and others, Reference Joughin, Smith and Holland2010; Favier and others, Reference Favier2014) or parameterizations which include a dependence on the ice-shelf slope (e.g. Little and others, Reference Little, Goldberg, Gnanadesikan and Oppenheimer2012; Lazeroms and others, Reference Lazeroms, Jenkins, Gudmundsson and van de Wal2018, Reference Lazeroms, Jenkins, Rienstra and van de Wal2019; Favier and others, Reference Favier2019).
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:
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) 2 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, Reference Joughin, Smith and Holland2010; Favier and others, Reference Favier2014). 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) 3 is a strong simplification of the melt rate derived by Lazeroms and others (Reference Lazeroms, Jenkins, Gudmundsson and van de Wal2018), who derive a general parameterization of the melt rates obtained with a plume model (Jenkins, Reference Jenkins1991). 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, Reference Van der Veen1998; Dupont and Alley, Reference Dupont and Alley2005; Benn and others, Reference Benn, Warren and Mottram2007b; Goldberg and others, Reference Goldberg, Holland and Schoof2009; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Gagliardini and others, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010; Jamieson and others, Reference Jamieson2012; Gudmundsson and others, Reference Gudmundsson, Krug, Durand, Favier and Gagliardini2012).
The most numerically convenient of these choices is the assumption of a fixed calving front position x c = const. (e.g. Gudmundsson and others, Reference Gudmundsson, Krug, Durand, Favier and Gagliardini2012; Arthern and Williams, Reference Arthern and Williams2017) 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, Reference Vieli, Funk and Blatter2000; Benn and others, Reference Benn, Hulton and Mottram2007a; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2010; Choi and others, Reference Choi, Morlighem, Rignot and Wood2021). A third calving choice we will consider assumes a fixed ice-shelf length, i.e. assuming x c = x g + L 0 with L 0 = 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:
with the first option describing a constant ice-shelf length L 0, 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
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, Reference Schoof2007b; Haseloff and Sergienko, Reference Haseloff and Sergienko2018; Sergienko and Wingham, Reference Sergienko and Wingham2022):
This is essentially a ‘shallow ice’ model (Fowler and Larson, Reference Fowler and Larson1978; Morland and Johnson, Reference Morland and Johnson1980), and the boundary conditions of this system are:
with (9c) effectively replacing the stress boundary condition (4c) (see, e.g. Schoof, Reference Schoof2007b, 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):
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, Reference Hindmarsh2012; Sergienko and Wingham, Reference Sergienko and Wingham2022), i.e.
where
and
Θ 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
Substituting this expression into (11a) and rearranging terms leads to an implicit expression for the ice flux at the grounding line
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 (Reference Schoof2007a, Reference Schoofb)
as a limiting case. Writing (13) in similar form, one can obtain
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, Reference Tsai, Stewart and Thompson2015; Brondex and others, Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017; Sergienko and Wingham, Reference Sergienko and Wingham2019). 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, Reference Reese, Winkelmann and Gudmundsson2018b). 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) 2 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:
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 (Reference Schoof2007a) and Robel and others (Reference Robel, Schoof and Tziperman2014), 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) analytical solution of (10) and (13) with suitably constructed Θ (e.g. (16), (18) and (25))
(2) semi-analytical solution of (10) and (13) with numerically determined Θ
(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) 1, 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) 2 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:
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 0. 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 γ1, γ2 and γ3 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 a1–c3 are for the same average melt rates, which are indicated in panels a4–c4.
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 a4–c4. Note that the y-axis of panel b4 is logarithmic, as the ice thickness remains non-zero for all γ2 ≤ 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 γ1 < −q g/L 0. 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 0, and as this limit is approached for γ2 → −∞, the ice thickness at the calving front asymptotically approaches zero (Fig. 2b4).
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. 2c4). Closer examination of these solutions (Appendix B.1) shows that this is only an apparent lower bound: for γ3 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 2a2: 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 0 ≈ 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 a2, b1 and c2, respectively). All other parameter values are listed in Table 1.
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 a4–c4 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 a1–a3, b1, and c1–c2.
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 (Reference Haseloff and Sergienko2018):
For a fixed ice thickness at the calving front (6) 3, 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, Reference Haseloff and Sergienko2018)
h c,unconf is the ice thickness at the calving front of unconfined ice shelves
and h c,conf is the ice thickness at the calving front of strongly buttressed ice shelves
For illustrative purposes, we solve (8)–(9) numerically with Matlab ODE solvers to obtain the ice-sheet profiles, and (1b) and (3) 2 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. 3a1–a3), 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. 3b1), as the flux appears to be monotonically increasing over the entire domain (Fig. 3b4). 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. 3c1–c2), one close to the unbuttressed grounding-line solution with a very short ice shelf (Fig. 3c1) and one on the retrograde part of the bed (Fig. 3c2). 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, Reference Haseloff and Sergienko2018).
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, Reference MacAyeal1989).
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))
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)
The parameters A 1–A 4 are defined in (C12), the parameters B and B 1 are defined in (C14) and ζ is given by (C11).
This condition is significantly more complex than the stability conditions derived in Schoof (Reference Schoof2012) and Sergienko and Wingham (Reference Sergienko and Wingham2022), 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, Reference Schoof2012). 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, Reference Schoof2007b, Reference Schoof2012).
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 1 > 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 (a1 and a3) as stable grounding-line positions, while the steady-state grounding-line position on the retrograde part of the bed (a2) 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 (b1) 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, Reference Haseloff and Sergienko2018, for details):
Equation (22a) is plotted as a dashed yellow line in Figures 3a4. As it qualitatively agrees with the other solutions, we can use (22a) to better understand the behaviour shown in Figure 3a1–a3: (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 (Reference Schoof2007b) (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:
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 3b4.
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
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)):
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 (Reference Schoof2012), 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, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010) 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):
τ 0 is the grounding-line stress of an unconfined ice shelf, given by (11c) 3, and q(x) is the steady-state flux on the ice shelf:
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:
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−1). We start our calculations from the stable steady state shown in Figure 3a1 and keep the ice shelf length constant (L 0 = 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
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.
Influence of relative melt location x r = (x m − x g)/L 0 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:
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 (γ1 = γ2 = γ3 = 0). We then vary the melt rate parameters γ1, γ2 or γ3, 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 a4–c4 (but not in the profiles), even though these correspond to freeze-on, and the models are not applicable in this parameter range.
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 a1–c3) 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 a4–c4. 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 a1–c1) and different calving laws (compare for example panels b1–b3). 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 a3–c3) 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 a1 and a2). The same is true for the slope-dependent melt rate when the melt rate parameter γ3 exceeds the critical value γc (panels c1 and c2).
γ1, γ2 and γ3 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 3a4–c4. 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−1 to $\bar {m}\approx -1$
m a−1 (ignoring positive values, as these are only included for completeness).
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 a1–a2 and c1–c2, 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 γ3. 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 (Reference Schoof2007b), Haseloff and Sergienko (Reference Haseloff and Sergienko2018) and Sergienko and Wingham (Reference Sergienko and Wingham2022), 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 5a4, c4).
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 (Reference Schoof2012). 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, Reference Schoof, Davis and Popa2017).
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, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010), 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) 3), 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 γ3, which regulates the strength of melting. If γ3 < γ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, Reference Jenkins1991; Little and others, Reference Little, Goldberg, Gnanadesikan and Oppenheimer2012; Sergienko and others, Reference Sergienko, Goldberg and Little2013; Slater and others, Reference Slater, Nienow, Goldberg, Cowton and Sole2017). 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, Reference Sergienko, Goldberg and Little2013). 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) 3, 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, Reference Van der Veen1996). 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, Reference Gladish, Holland, Holland and Price2012; Sergienko, Reference Sergienko2013; Langley and others, Reference Langley2014).
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, Reference Reese, Gudmundsson, Levermann and Winkelmann2018a; Zhang and others, Reference Zhang, Price, Hoffman, Perego and Asay-Davis2020), 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 (Reference Schoof2007a), Haseloff and Sergienko (Reference Haseloff and Sergienko2018) and Sergienko and Wingham (Reference Sergienko and Wingham2022). 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, Reference Gomez, Mitrovica, Huybers and Clark2010; Robel and others, Reference Robel, Schoof and Tziperman2014, Reference Robel, Schoof and Tziperman2016; Brondex and others, Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017; Åkesson and others, Reference Åkesson, Nisancioglu and Nick2018; Reese and others, Reference Reese, Winkelmann and Gudmundsson2018b; Sergienko and Wingham, Reference Sergienko and Wingham2019, Reference Sergienko and Wingham2022).
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.
Acknowledgements
MH acknowledges financial support through a VC Senior Fellowship at Northumbria University. OS is supported by award NA18OAR4320123 from the National Oceanic and Atmospheric Administration, USDepartment of Commerce. The statements, findings, conclusions and recommendations are those of the authors and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration, or the US Department of Commerce. We thank an anonymous reviewer and T. Zhang whose comments helped to improve and clarify the manuscript.
APPENDIX A. Numerical confirmation
Figure 7 shows a comparison of equilibrium grounding-line positions for varying widths obtained with the different methods listed in Section 2.3.2. The grey area indicates a retrograde bed. Note that we can only find stable steady states with the numerical method, while we can calculate stable and unstable steady state with the analytical and semi-analytical methods described above. The semi-analytical and numerical results agree well with each other. The analytic results are a good match for the case of a constant ice-shelf length and a constant calving front position, but deviate from the case with a constant ice-shelf front. This is most likely due to the approximation of the ice thickness at the calving front with (19a)–(19c) which is valid in the limits of no and strong buttressing only.
Comparison of grounding line positions for different ice-shelf widths obtained with the different methods (i)–(iii) outlined in the text. Analytic results are obtained from solution of (10) with (14), semi-analytic results are obtained from solution of (10) with (22a)–(22c) and numerical results are obtained from solution of (1)–(4) with Comsol. The grey-shaded area marks the upward-sloping part of the bed. Panel a: grounding-line positions for different ice-shelf widths with prescribed ice-shelf length. Panel b: grounding-line positions for different ice-shelf widths with prescribed calving front position. Panel c: grounding-line positions for different ice-shelf widths with prescribed ice thickness at the calving front. Panel d: bed elevation.
The grounding-line positions obtained for large values of the ice-shelf width W correspond to the unconfined, unbuttressed limit. Generally, as the width W decreases, buttressing increases, leading to equilibrium grounding-line positions that are downstream of the unconfined limit. However, the shape of the W–x g relationship is markedly different for different calving laws, as expected from other studies (Schoof and others, Reference Schoof, Davis and Popa2017; Haseloff and Sergienko, Reference Haseloff and Sergienko2018).
APPENDIX B. Analytic ice-shelf solutions
B.1. Unconfined ice shelves with slope-dependent melting
To derive an analytic solution, it is convenient to non-dimensionalize (1b), (3) 2, (4c) and (17) by setting x = L 0 x* + x g, h = h gh*, $u = q_{\rm g} h_{\rm g}^{-1} u^{\ast }$
and $\gamma _3 = q_{\rm g} h_{\rm g}^{-2}\gamma _3^{\ast }$
. We neglect asterisks rightaway in equations and immediately integrate (1b) with (4c) (see e.g. MacAyeal and Barcilon, Reference MacAyeal and Barcilon1988; Schoof, Reference Schoof2007b), leaving the simplified ice-shelf system for ε = 0:
with
the ratio between the scale of the extensional stress at the grounding line and the scale of the driving stress at the grounding line and the initial conditions (17):
Assuming dh/dx < 0, we can integrate (B1) 2 by use of (B3) which provides the ice velocity as a function of ice thickness
Using the product rule we can rewrite (B1) as an implicit integral equation for h:
which can be simplified to
Integration gives
and re-dimensionalization gives
The left-hand side of this equation is non-linear in h/h g (see Fig. 8a), and solutions of this equation require x = L s/L 0 < 1 below a critical γ3 = γc. To find γc we determine the ice thickness h = h crit that minimizes (B4) by taking the derivative with respect to h. This gives
Re-dimensionalization and substituting h crit into (B5) gives
For n = 3, it is possible to write solutions to (B5) in closed form, and the ice thickness at the ice-shelf front is given by:
We compare these analytic solutions for the ice thickness at the ice-shelf front h(x = 0) = h c with numerical solutions in Figure 8, panels b and c. Note that we plot the ‘apparent’ ice thickness at the shelf front for solutions with ε > 0 in panel b.
Examples for unbuttressed profiles with slope-dependent melting $\dot m = \gamma _3 ( h_{\rm g}-h) \vert {{\rm d} {h}}/{{\rm d} {x}}\vert \times ( 1 + \varepsilon ^{2} ( {{\rm d} {h}}/{{\rm d} {x}}) ^{2}) ^{-1/2}$
. If γ3 < γc, the ice thickness at the ice-shelf edge goes to zero over a short distance, as the melt rate becomes singular when u = −γ3(h g − h).
B.2. Grounding-line backstress for location-dependent melting $\dot m = f( x)$
To find solution for (23a), we closely follow Haseloff and Sergienko (Reference Haseloff and Sergienko2018), but assume a variable melt rate $\dot m = f( x)$
. It is convenient to non-dimensionalize (1b), (3) 2, (4c) and (17) by setting x = L sx* + x g, h = h gh*, $u = q_{\rm g} h_{\rm g}^{-1} u^{\ast }$
and $\dot m = q_{\rm g} L_{\rm s}^{-1}\dot {m}^{\ast }$
. We neglect asterisks rightaway and obtain the non-dimensional system of ice-shelf equations
with boundary conditions
where subscripts indicate derivatives. We have introduced the non-dimensional groups η and β
and will assume η ≪ β ~ 1, so that we can simplify (B8a):
Integration of (B8b) with (B8c) 1 is straightforward:
Using $M( x) = \int _0^{x} f( x')\; {\rm d} x'$
in (B9), we obtain:
For the flux to remain positive over the length of the ice shelf, we require (1 + M(x)) > 0 for 0 < x ≤ 1. Requiring this, we can integrate once more
where $h_1 = [ \eta ^{n} \beta ( 1 + M( 1) ) ^{1/n + 1}/2 ] ^{n/( n + 1) ^{2}}$
is the ice thickness at the ice-shelf front (see Haseloff and Sergienko, Reference Haseloff and Sergienko2018, for details). From (B12a) with (B10) we can obtain an expression for the ice-shelf velocities in the main ice-shelf body:
This velocity has to match the reduction in buttressing experienced at the grounding line (again, see Haseloff and Sergienko, Reference Haseloff and Sergienko2018, for details):
so that we obtain
Note that $1 + M^{\ast }( x^{\ast }) = [ q_{\rm g} + \int _{x_{\rm g}}^{x} f( x') {\rm d} x'] /q_{\rm g} = q( x) /q_{\rm g}$
. Re-dimensionalizing gives equation (23a).
APPENDIX C. Linear stability analysis
To perform a linear stability analysis, we use the time-dependent version of the reduced ice-sheet model (8) written in terms of h and q = uh
where subscripts indicate partial derivatives and we have introduced
The stress-continuity condition (13) with the backstress for a location-dependent melt rate $\dot m( x)$
(23a) becomes
where
and h satisfies the floatation condition
at the grounding line.
From the derivative of the floatation condition with respect to t we can also determine the rate of grounding-line migration:
by substituting h t from (C1b), and subsequently q x from (13) into (C3) (note that we leave the buttressing term in its general form):
To determine the stability of grounding-line positions obtained from (C2), we consider small time-dependent perturbations around steady-state solutions
where ɛ is a small parameter. At leading order, the linearized perturbation problem recovers the steady-state outer problem (8) for $\hat h$
and $\hat q$
. The linearized perturbation problem of O(ɛ) is
The perturbed flotation condition is
and the boundary conditions at the origin are
Equation (C2) is of the form
To arrive to a linearized, perturbed version of (C2) we therefore start by recognizing that
with $\hat h_x = {{\rm d} {\hat h}}/{{\rm d} {x}}$
and $\hat q_x = {{\rm d} {\hat q}}/{{\rm d} {x}}$
and
Introducing
the O(ɛ)-term on the left side of (C9) is
Before we can repeat similar steps for the right-hand side of (C9), we have to make a decision about the form of the calving law, which enters through ξ:
Complex calving laws might depend on the ice thickness, flux, and grounding-line position, i.e. x c = x c(x g, h g, q g), see for example (6). If we allow for a general calving law like this, then (C9) is of the form
with
and
We also used
C.1 Sturm–Liouville form and general stability condition
To show that (C6)–(C8) with (C13) is a Sturm–Liouville problem for $\tilde h$
, we assume the solution is separable, and can be written as $\tilde h( x,\; \, t) = \tilde h( x) {\rm e}^{\Lambda t}$
, where Λ is an eigenvalue, whose sign determines the stability. From (C1b)
and from (C7)
With a little algebra (largely following Sergienko and Wingham, Reference Sergienko and Wingham2022), we can thus rewrite (C6b) as
with
It is straightforward to show that (C17) can be put into Sturm–Liouville form (again following Schoof, Reference Schoof2012; Sergienko and Wingham, Reference Sergienko and Wingham2022)
with μ an integrating factor and
For the boundary condition at the grounding line, we substitute (C16) into (C13) and divide by $\tilde h$
, giving
As in Sergienko and Wingham (Reference Sergienko and Wingham2022), the functions μ, P and R are continuous and μP is a positive function, and the boundary conditions (C8) and (C19) can be written as homogeneous function of $\tilde h$
and $\tilde h_x$
. Importantly, as Sergienko and Wingham (Reference Sergienko and Wingham2022) elaborate on, it is possible to show that for a Sturm–Liouville problem as the above, the eigenfunction corresponding to the largest eigenvalue does not change sign, provided [A 1 + B(c 1 + c 3∂x c/∂q)] > 0 (the parameter A 3 is positive by definition). A 1 is positive as long as $( {1\over n} + 1) \gamma _{\rm w}\hat q^{1/n}h^{m + 1} + ( m + 1) \gamma _b\hat q^{m} \hat h^{1/n} + \hat h^{1/n + m + 1}b_x> 0$
, B > 0 by definition, and the sign of ∂x c/∂q depends on the calving law. If we reasonably assume that the steady solution is to remain grounded upstream of the grounding line, then $\hat h_x + b_x/( 1-\delta ) < 0$
(Schoof, Reference Schoof2012).
Under these circumstances, Eqn (C19) determines the stability of a grounding-line position determined from (13):
where the third condition $\hat \zeta > 0$
has to be satisfied for steady-state solutions to exist. Note that the choice of the calving law enters through ∂x c/∂h g, ∂x c/∂x g and ∂x c/∂q g which makes it possible for different calving choices to alter the stability of a steady state from stable to unstable.
C.2 Relationship with flux gradient
For comparison with the existing literature (Schoof, Reference Schoof2012; Sergienko and Wingham, Reference Sergienko and Wingham2022) it is convenient to express the stability condition as a function of the gradient of the grounding-line flux dq g/dx g. The stress condition (C2) provides the grounding-line flux for Θ given by (23a):
Taking the derivative with respect to x g and using the definitions of A 1–A 4 (C12), B and B 1(C14), evaluated at x = x g, gives
with
Reordering gives:
with
We can identify dh g/dx g = −b x/(1 − δ), $\hat q_x = q_x = \dot a$
, $\hat q = q_{\rm g}$
, $\hat h = h_{\rm g}$
, b xx = db x/dx g, $\hat q_{xx} = {{\rm d} {\dot a}}/{{\rm d} {x_{\rm g}}}$
so that the stability condition (C19) is linked to the flux gradient through
and we can write (C20) as
The same analysis can be applied to stability conditions derived by Sergienko and Wingham (Reference Sergienko and Wingham2022) and reduce them to the form (C26).
Open access
