Diverse behaviors of marine ice sheets in response to temporal variability of the atmospheric and basal conditions

Abstract

The observed retreat of the grounding line of the present-day ice sheets and the simulated grounding-line retreat of ice sheets under changing climate conditions are often interpreted as indications of marine ice-sheet instability (MISI), driven by a positive feedback between the ice discharge and conditions at the grounding line. However, the arguments that support this feedback are valid only for steady-state conditions. Here, we assess how unconfined marine ice sheets may behave if atmospheric conditions and basal conditions evolve with time. We find that the behavior of grounding lines can exhibit a range from unstoppable advance and retreat to irregular oscillation irrespective of the stability of the corresponding steady-state configurations obtained with time-invariant conditions. Our results show that numerical simulations with a parameterization of the ice flux through the grounding line used in large-scale ice-sheet models produce markedly different results from simulations without the parameterization. Our analysis demonstrates that the grounding-line migration can be driven by the temporal variability in the atmospheric and basal conditions and not by MISI, which assumes unchanging conditions. Instead, the grounding-line advance or retreat is determined by interactions between ice flow, basal processes and environmental conditions throughout the length of a marine ice sheet in addition to the circumstances at its grounding line.

1. Introduction

The contributions of marine ice sheets to sea level are controlled by the dynamics of their grounding lines. Typically, migration of grounding lines on bedrock that slopes toward interior of the ice sheet is thought to be caused by marine ice-sheet instability (MISI) – a hypothesis proposed by Weertman (1974). According to Weertman's hypothesis, the stability of a steady-state marine ice sheet is determined by the bed slope at the location of the grounding line: if the slope is ‘retrograde’, i.e. the bed slopes toward the interior, the ice sheet is inherently unstable; if the slope is ‘prograde’, i.e. the bed slopes away from the interior, the ice sheet is unconditionally stable. As the West Antarctic ice sheet, and many parts of the East Antarctic and Greenland ice sheets rest on beds with retrograde slopes, the behavior of their grounding lines is described ‘unstable’ (e.g. Shepherd and others, 2018a) as a corollary of Weertman's result.

The original hypothesis (Weertman, 1974) and its consequent analysis (Schoof, 2007a, 2007b, 2012) define stability as a property of steady states of ice sheets, i.e. all their environmental conditions and internal properties (e.g. basal sliding) do not change with time. While many studies broadened the definition of MISI and term ‘instability’ as any positive feedback between the grounding-line retreat and increase of ice flux (or discharge) through the grounding line (Pattyn and Morlighem, 2020), they use similar arguments to those first proposed by Weertman (1974) and Schoof (2007b) for steady states. These arguments have been used to explain the observed retreats of the grounding lines of Antarctic and Greenland ice sheets (Rignot, 1998; Shepherd and others, 2018a; Khan and others, 2020). Similarly, the grounding-line retreat produced in simulations of the future ice-sheet behavior under projected climate conditions changing in time has also been interpreted as an indication of MISI (Cornford and others, 2015; DeConto and Pollard, 2016).

Recent studies that considered steady-state configurations for laterally confined marine ice sheets (Gudmundsson and others, 2012; Kowal and others, 2016; Haseloff and Sergienko, 2018; Pegler, 2018; Reese and others, 2018; Sergienko and Wingham, 2022; Sergienko, 2022a), non-negligible bed topography (Sergienko and Wingham, 2022) and the regime of low basal stress (Sergienko and Wingham, 2019) have demonstrated that the bed slope alone does not necessarily determine stability of steady-state marine ice sheets, and in particular configurations that can be stable and unstable with their grounding lines located on either prograde or retrograde beds (Gudmundsson, 2013; Haseloff and Sergienko, 2018, 2022; Sergienko and Wingham, 2019, 2022). A study by Sergienko and Haseloff (2023) that considered a laterally confined marine ice sheet that experiences temporally variable submarine melting has found that the grounding line can intermittently advance and retreat as well as retreat in an unstoppable manner, even though a steady state obtained with time-averaged submarine melt rates is stable.

Here, we use the same 1-D model of an unconfined marine ice sheet resting on a smooth bed topography which has been used to establish stability conditions of steady-state marine ice sheets (Schoof, 2007a, 2007b, 2012), and subject it to time-evolving environmental conditions with the goal to investigate the marine ice-sheet response to changing conditions. In a laterally unconfined configuration, the problem reduces to the grounded part only, with submarine melting having no effects on the grounding line. We initialize time-variant simulations with two kinds of steady-state configurations (Fig. 1) both of which conform to the MISI hypothesis, i.e. the one shown with a green line, whose grounding line is located on the prograde slope, is stable when subject to small perturbations from its steady-state position; and the second, shown with a blue line, whose grounding line is located on the retrograde slope, and is unstable when subject to small perturbations.

Figure 1. Ice-sheet configurations: x _d, the ice divide location; x _g, the grounding-line location; x _c, the calving front location; b(x), the bed elevation; sea level lies at zero elevation (dot-dash blue line). The green line indicates a stable steady-state configuration; the blue line indicates an unstable steady-state configuration. The black line indicates bed topography.

Our results show that when accumulation rate (external) or basal sliding (internal) conditions change with time, marine ice sheets could persist or disappear irrespective of the stability of the steady states obtained with the time-averaged conditions. We illustrate with time-variant examples that the ice-sheet mass balance is different from that in a steady state; the partitioning between its terms is not obvious, and in consequence the grounding-line migration need not result in a sustained advance or retreat on retrograde beds or stable behavior on prograde beds.

Using simplified assumptions of negligible bed slope and accumulation rates in the vicinity of the grounding line, Schoof (2007b) has derived an expression for the steady-state ice flux as a function of the ice thickness at the grounding line. Due to its simplicity, this expression is used in a variety of applications (e.g. simplified conceptual models, Robel and others, 2018; analysis of ice-sheet wide observations, Slater and Straneo, 2022). It is also used as a parameterization in several large-scale ice-sheet models (e.g. DeConto and Pollard, 2016; Pattyn, 2017; Quiquet and others, 2018). However, the expression is a statement that at the grounding line, the internal deformation equals the ice advection, and, as we illustrate, the imbalance of these terms contributes to the rate of grounding-line motion. The results of simulations with and without the parameterization of the ice flux at the grounding line are significantly different, demonstrating its unsuitability for time-variant conditions.

The article is organized as follows. In Section 2 we provide a description of the model and numerical methods. Section 3.1 demonstrates the effects of time variability in the surface accumulation and Section 3.2 in the sliding conditions. In Section 3.3 we examine the performance of the ice-flux parameterization and provide physical interpretations of the results. We give our conclusions in Section 4. Readers with less interest in the mathematical and numerical aspects can proceed directly to Sections 3 and 4.

2. Methods

2.1 Model description

The model is the same as one used to investigate steady-state configurations of marine ice sheets (Schoof, 2007a, 2007b). Here, we provide its brief description. Flow of an unconfined ice stream into an unconfined ice shelf (Fig. 1) can be described by vertically integrated momentum balance under assumptions of negligible vertical shear appropriate for ice stream and ice shelf flow (MacAyeal, 1989), that is,

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

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

where u(x) is the depth-averaged ice velocity, h(x) is the ice thickness, b(x) is the bed elevation (negative below sea level and positive above sea level), A is the ice stiffness parameter (assumed to be constant), n is an exponent of Glen's flow law (n = 3), g is the acceleration due to gravity, τ _b is basal shear, $g^\prime$ is the reduced gravity defined as

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

where

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

is the buoyancy parameter, ρ and ρ _w are the densities of ice and water, respectively. x _d is the location of the ice divide, x _c is the location of the calving front and x _g is the location of the grounding line. The basal shear is assumed to follow a power law:

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

where C is the sliding parameter and m = 1/n is the sliding exponent.

The mass balance is

(5)$$h_t + \left(uh\right)_x = \left\{\matrix{& \dot a & \qquad & 0\leq x\leq x_{\rm g},\; \cr & \dot m & \qquad & x_{\rm g}< x\leq x_{\rm c},\; }\right.$$

where $\dot a$ is the net accumulation/ablation rate, usually referred to as the surface mass balance (SMB) of the ice stream, and $\dot m$ is the net accumulation and submarine melting rate of the ice shelf.

The boundary conditions at the divide x _d and the calving front x _c are

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

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

At the grounding line x _g the continuity conditions

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

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

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

(where τ = 2 A ^−(1/n)h|u _x|^1/(n−1)u _x is the longitudinal stress) and the flotation condition

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

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

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

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

In circumstances where ice shelves are unconfined, the momentum balance of the ice shelf (1b) can be integrated with the boundary condition at the calving front, x _c, (6b) and the continuity conditions (7), and the problem can be reduced to the ice-stream part only with the boundary conditions at the grounding line – the flotation condition (8) and the stress condition:

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

The rate of the grounding-line migration can be obtained by taking the total time derivative of the flotation condition (8) and rearranging terms:

(11)$$\dot x_{\rm g} = -{h_t + {b_t}/{( 1-\delta) }\over h_x + {b_x}/{( 1-\delta) }},\; $$

where b _t is the rate of change of the bed elevation that can be due to subglacial morphological processes (e.g. erosion or sediment deposition), or due to glacial isostatic adjustment, or due to changes in sea level. Here, we do not take into account these processes and assume b _t = 0. Using the mass-balance equation, this expression becomes (19).

2.2 Numerical implementation

We solve numerically the system of equations describing the evolution of the grounded part of the marine ice sheet and its flow. The system includes the momentum, Eqn (1a), and the mass, Eqn (5), balances with the boundary conditions (6a), (8) and (10), and is solved using the finite-element solver Comsol^TM (COMSOL, 2024). In all simulations, the grid resolution is spatially variable: it is 200 m through 95% of the length of the domain, and 1 m in the 5% closest to the grounding-line position. The initial steady-state configuration is obtained using a minimization procedure based on the bound optimization by quadratic approximation optimization algorithm (Powell, 2009). The time-variant simulations are performed on domains with a moving boundary, the grounding line. This is done using an arbitrary Lagrangian–Eulerian method (Donea and others, 2017). This boundary moves with a prescribed velocity, expression (19).

We use the following model setup. The bed topography is described by $\displaystyle {b( x) = b_0 + b_{\rm a} \cos {( 2\pi x}/{L}) }$, with b ₀ = −500 m, b _a = 250 and L = 1000 km, $\dot a_{ss} = 0.1$ m a⁻¹, the sliding law parameters C ₀ = 7.6 × 10⁶ Pa s^1/3 m^−1/3 and m (chosen to be m = 1/n), and ice stiffness parameter A = 1.35 × 10⁻²⁵ Pa⁻³ s⁻¹ (which corresponds to $T_{\rm ice}\approx -20^\circ$C). With the chosen bed elevation, we consider the largest possible extent of the ice sheet to be 1000 km.

All steady-state configurations used in this study conform with the MISI hypothesis.

2.3 Model experiments

To examine the marine ice-sheet behavior in response to the time-varying accumulation rate and time-varying basal conditions, we perform two sets of time-variant experiments. The first scenario aims to mimic the effect of changing climate conditions – atmospheric temperature and hence the SMB. The second scenario aims to mimic possible changes in basal conditions internal to the ice sheet, caused by, for instance, subglacial processes.

2.3.1 Stochastically varying atmospheric conditions

This set of experiments aims to resemble the effects of changing climate conditions on the dynamics of marine ice sheets. This is done by varying SMB in time and with the distance along the ice sheet. The $\dot {a}$ is determined by atmospheric conditions. If the atmospheric temperature is below the freezing point, snow mass accumulates on the surface; as the temperature approaches and exceeds the freezing point, mass is lost through ablation due to sublimation and melting. The atmospheric temperature decreases as elevation increases, and even under climate warming the higher elevations may experience net accumulation, whereas lower elevations may experience net ablation. Thus atmospheric temperature, which is controlled by the climate conditions, can be used as a proxy for the SMB. Here, we use an empirical relationship between $\dot a$ and atmospheric temperature at the ice-sheet surface derived by Sergienko (2022b) who analyzed the results of regional climate model simulations for the Antarctic and Greenland ice sheets for projected climate conditions under a scenario in which emissions continue to rise throughout the 21st century (IPCC, 2013).

This empirical expression relates $\dot a$ to the atmospheric temperature at the ice-sheet surface T _S:

(12)$$\dot a( T_{\rm S}) ( x,\; t) = a_1 \exp{\left[- {( T_{\rm S}( x,\; t) -T_0) ^2\over 2\sigma^2}\right]}-a_2\exp{\left[-2{T_{\rm S}( x,\; t) -T_0\over T_0}\right]},\; $$

where a ₁ = 2.4 m a⁻¹, a ₂ = 0.8 m a⁻¹, $T_0 = -15^\circ$C and $\sigma = 6^\circ$C are empirical parameters and T _S(x,  t) is the temperature at the surface elevation S, which is

(13)$$T_{\rm S}( x,\; t) = T^{\rm sl}( t) -\Gamma S( x,\; t) ,\; $$

where $\Gamma = 9.8^\circ$C km⁻¹ is the lapse rate, assumed to be adiabatic in this study, and T ^sl(t) is the temperature at sea level.

A number of previous numerical studies investigating the response of grounding lines to variability in climate forcing using realistic (Hoffman and others, 2019; Robel and others, 2019) and idealized configurations (Christian and others, 2022; Felikson and others, 2022) have demonstrated that variability in the climate forcing causes the grounding line to behave differently from that resulting from time-invariant forcing. Ice-core records indicate that the climate of polar regions exhibits variability on a variety of temporal scales (Jouzel and others, 2007a, 2007b, 2013), which range from hundreds of thousands of years governed by orbital cycles (Milanković, 1941) to decades and years governed by variability in atmospheric and oceanic circulations such as the Southern Annular Mode and El Ninõ-Southern Oscillation (Kim and others, 2020). These ice-core records are also characterized by noise at all temporal scales. To capture the observed variability, we choose to represent the temporal evolution of $\dot {a}$ as a noise function of time with decadal and centennial correlation times, i.e. we assume that T ^sl(t) varies with time according to

(14)$$T_{\rm sl}( t) = T^{\rm sl}_{0} + T^{\rm sl}_{10} N\left({t\over T_{10}}\right) + T^{\rm sl}_{100} N\left({t\over T_{100}}\right),\; $$

where $T^{\rm sl}_{0}$ is a steady-state value of atmospheric temperature at sea level, which was used to compute steady-state configurations of the ice sheet that are used as initial conditions for time-variant simulations; $T^{\rm sl}_{10} = 1.25^\circ$C is the amplitude of the decadal variability and $T^{\rm sl}_{100} = 2.5^\circ$C is the amplitude of the centennial variability, respectively; N(t) is a noise function with a uniform distribution and zero mean value, T ₁₀ = 10 years is the decadal and T ₁₀₀ = 100 years is the centennial correlation timescale. The choice of these timescales and respective magnitudes are motivated by analyses of ice-core records (e.g. Kobashi and others, 2010; Thomas and others, 2013). We restrict our model to decadal and centennial timescales because introducing longer, millennial scales would require simulations in excess of 100 ka, that are run here. For all experiments we perform five simulations with different seeds in the noise functions, which results in 30 experiments in total.

2.3.2 Periodic variability of basal conditions

Our simulations with time-evolving basal conditions aim to capture the consequences of subglacial processes on the ice flow in the ice-sheet interior. Inferences of basal conditions beneath both Antarctic and Greenland ice sheets, made from radar observations (Schroeder and others, 2013) and using inverse method techniques (Sergienko and others, 2008; Morlighem and others, 2013; Sergienko and Hindmarsh, 2013; Sergienko and others, 2014) indicate that these conditions are highly heterogeneous and can vary by many orders of magnitudes. This variability is attributed to a wide range of processes operating on the wide range of timescales – from the rapid flow of subglacial water (Wingham and others, 2006) to the formation of subglacial landforms (King and others, 2007). In the absence of direct or indirect estimates of the characteristic timescales of such processes, we choose to investigate the effects of changing basal conditions by imposing periodic variability on the sliding parameter with periods ranging from 25 ka to 400 years. As we use the same model used by Schoof (2007a, 2007b, 2012) the sliding law is in the form of (4). While all other parameters remain constant, the sliding parameter C evolves with time periodically:

(15)$$C( x,\; t) = 10^{\alpha}C_0,\; \quad\alpha = {k_t{x_{0}-x\over x_{\rm g}}}\sin{2\pi t\over T} ,\; $$

where C ₀ is the steady-state value of the sliding parameter that was used to compute the corresponding steady-state configurations (we use the same value used by Schoof, 2007b); k _t is the amplitude of the order of magnitude variability, x ₀ is a ‘catchment extent’ that affects the grounding line downstream of it; T is the period of cyclic variability. We have chosen this model to include a variation in sliding as a function of position in addition to a variation in time. The values of k _t that we use are such that the value of C produced by Eqn (15) and the corresponding basal shear stress are within the range of values obtained for the present-day ice sheets using inversion techniques (e.g. Sergienko and others, 2008; Morlighem and others, 2013). In contrast to the experiments with time-evolving SMB, we do not consider stochastic variability due to lack of knowledge of any such characteristics. To focus on the effects of temporal variability in basal conditions we keep all other parameters constant in space and time and use $\dot a = 0.1$ m a⁻¹.

The design of these experiments reflects the current state of the knowledge: much more is known about the temporal variability of atmospheric conditions than of basal conditions. Consequently, the first scenario is guided by the results of analyses of ice-core records (e.g. Kobashi and others, 2010; Thomas and others, 2013). However, there are no direct observations of the temporal variability of basal conditions; consequently the second scenario is highly idealized. In both sets of experiments all other parameters remain constant in space and time.

2.3.3 Experiments with the steady-state grounding-line flux formula

Additionally, we perform experiments described above with a parameterization of the grounding-line stress condition which is based on the widely used expression of the steady-state ice flux at the grounding line obtained by Schoof (2007b). In a steady state $\dot x_{\rm g} = 0$, and if the bed slope b _x and the accumulation rate $\dot a$ at the grounding line can be neglected, the internal deformation at the grounding line and ice advection balance each other. For these circumstances, Schoof (2007b) formulated an approximate expression for the ice flux at the grounding line:

(16)$$q_{{\rm g}S} = \left({A\left(\rho g\right)^{n + 1}\delta^n\over 4^nC}\right)^{{1}/{( m + 1) }}h^{{( m + n + 3) }/{( m + 1) }}.$$

We repeat simulations with the time-variant accumulation rate and basal sliding using the ice-flux expression (16) as a boundary condition. The ice velocity at the grounding line

(17)$$u_{\rm g} = \left({A\left(\rho g\right)^{n + 1}\delta^n\over 4^nC}\right)^{{1}/{( m + 1) }}h^{{( n + 2) }/{( m + 1) }},\; \quad x = x_{\rm g},\; $$

is used as a boundary condition instead of stress condition (10). All other parameters and conditions are identical to numerical simulations described above. We compare q _gS to the ice flux q _g computed in simulations with stress condition (10) at the grounding line, and which is given simply by

(18)$$q_{\rm g} = uH,\; \quad x = x_{\rm g}.$$

2.4 Model analysis

In order to understand what governs the motion of its grounding line, we analyze the rate of the grounding-line migration $\dot x_{\rm g}$:

(19)$$\dot{x}_{\rm g} = ( uh_x + u_xh-\dot a ) / \left(h_x + {b_x\over 1-\delta}\right).$$

For brevity, we denote the denominator $\displaystyle {D = h_x + {b_x}/{( 1-\delta ) }}$. In this expression, the three terms of the numerator are all contributions to rate of change of height h _t at the grounding line, due respectively to the advection of ice from upstream, the internal deformation at the grounding line and the accumulation at the grounding line; the denominator translates this rate to the corresponding horizontal velocity of the grounding line. The last two terms are determined by the local conditions at the grounding line. The accumulation term is determined solely by conditions at the grounding line, and, if the flow enters an unconfined ice shelf, this is true too of the ice deformation term, as in this case it is balanced by the pressure deficit. In contrast, the first term is determined by the ice flow along the length of the ice stream, and reflects the integrated effects of the accumulation, changes of the ice thickness and basal conditions of the grounded part of the marine ice sheet. This expression indicates that in a steady state ($\dot x_{\rm g} = 0$) the accumulation, ice advection and internal deformation at the grounding line balance each other. Generally, however, the grounding line migrates due to imbalance of these terms.

We also analyze the integrated form of the grounded ice sheet mass balance (5), i.e.

(20)$$\int_{x_{\rm d}}^{x_{\rm g}} {\rm d}x\left[h_t + ( uh) _x\right] = \int_{x_{\rm d}}^{x_{\rm g}}{\rm d}x \dot a.$$

Taking into account the boundary conditions at the ice divide x _d and recognizing that $\left.uh\right \vert _{x = x_{\rm g}} = q( x_{\rm g})$, the above expression can be written as

(21)$$q( x_{\rm g}) = \int_{x_{\rm d}}^{x_{\rm g}}\dot a( x) {\rm d}x-\int_{x_{\rm d}}^{x_{\rm g}}h_t( x) {\rm d}x.$$

In our analysis we use the form (21) of the integrated mass balance of the grounded part of the ice sheet.

3. Results

3.1 Time-evolving SMB

In response to temporal variations in the accumulation, the simulated ice sheets exhibit diverse dynamic behaviors, which are illustrated in Figure 2. In this figure, simulations that are initialized from positions on the prograde bed (illustrated by the green ice sheet in Fig. 1) are shown in the left-hand side panels; those initialized on the retrograde bed (illustrated by the blue ice sheet in Fig. 1) are shown in the right-hand side panels. The panels are arranged vertically according to their various ‘modes’ of behavior, which depend on the value of temperature at sea level $T_{0}^{\rm sl}$ in Eqn (14). Figures 2a, b illustrate retreats, in case 2a after a long duration of oscillatory behavior of retreat and growth, in case 2b from retrograde positions to a prograde positions; Figures 2c, d illustrate oscillatory behavior; while Figures 2e, f illustrate unstopped growth to the edge of the model domain. The duration of each of the plots is chosen to illustrate the character of their behavior. In the cases shown in Figures 2c, d, we extended the simulations to 100 ka (not shown) to confirm that the grounding-line behavior does not change on longer timescales than those shown in the figure.

Figure 2. Grounding-line response to variable accumulation. (a–f) Grounding-line positions x _g(t). All simulations were initialized with respective steady-state configurations and were performed with the respective values of $T^{\rm sl}_0$ (Eqns (12)–(13)). Left panels correspond to the initial configurations with the grounding-line positions on the prograde slope; right panels correspond to the initial configurations with the grounding-line positions on the retrograde slope. Colors represent simulations with different seeds in the noise function. The blue rectangle in panel a marks the 2000 years interval shown in Figure 3. The red boxes outlining panels a and d indicate simulations that are repeated with the ice-flux parameterization and described in Section 3.3.

In Figure 2, the sea-level temperatures $T_{0}^{\rm sl}$ that determine the initial steady states are given in the panels. There is no simple monotonic relationship between $T_{0}^{\rm sl}$ and the horizontal extent of the ice sheet. This is due to several factors that include the possibility of multiple steady-state configurations for the same set of parameters; the highly non-linear dependence of the ice-sheet thickness and the horizontal extent on $\dot a$; and the highly non-linear dependence of $\dot a$ on the surface temperature, which is a function of the ice-sheet surface elevation (Eqns (12) and (13)).

Among the behaviors shown in Figure 2 are those in accordance with the MISI hypothesis. In Figure 2b, the grounding lines move from their initial position on the retrograde bed slope to stable positions on the prograde bed slope; in Figure 2c, the grounding lines oscillate around a stable position on the prograde bed slope; and in Figure 2f the grounding lines continuously advance from its initial position on the retrograde bed slope (instability allowing for unstopped advance as much as retreat). However, and equally, there are three counter cases. Figure 2a shows ultimate extinctions from initial positions on the prograde bed slope; Figure 2d shows oscillations about stable locations on the retrograde bed slope and Figure 2e shows unstopped advances from an initial position on the prograde bed slope.

To get insight into what governs the behavior of the grounding line, we analyze the rate of the grounding-line migration $\dot x_{\rm g}$ for 2000 years of one simulation (the blue box in Fig. 2a). As Figure 3 illustrates, all the terms of the right-hand side of Eqn (19) have similar magnitudes. In addition to the immediate effect at the grounding line of the variability of the SMB (the term $\displaystyle {-\dot a / D}$), it appears in a more muted fashion in the ice advected from upstream (the term $\displaystyle {uh_x / D}$). The resulting rate of the grounding-line migration (the dark green line) is the imbalance between all these effects. As a result, the magnitude of the rate of the grounding-line migration is substantially smaller than the magnitudes of any of the individual terms.

Figure 3. Rate of the grounding-line migration the terms of Eqn (19) for the simulation described by the dark blue line in Figure 2a (the grounding line is on a prograde slope), during the 2000 year interval marked by the blue rectangle in panel a. Here, $\displaystyle {D = h_x + {b_x}/{( 1-\delta ) }}$.

The net effect of the three terms in Eqn (19) has no simple connection to the local conditions at the grounding line. The sign of the rate of the grounding-line migration $\dot x_{\rm g}$, Eqn (19), determines whether the grounding line advances (positive) or retreats (negative). In Eqn (19), the ice-thickness gradient h _x as well as ice velocity u depend on the ice flux q at x _g, which in turn depends on the integrated $\dot a$ and the rate of the ice-thickness change h _t throughout the extent of the ice sheet. As a result, the rate of the grounding-line migration $\dot x_{\rm g}$ depends on the size of the ice sheet, that is, the grounding-line position itself, in a complex, non-linear way.

In the example shown in Figure 2a, while the grounding line remains on the prograde slope (Fig. 1) throughout the course of the simulations, it also exhibits a long-term retreat, and, after ~20–60 ka, depending on the simulation, the ice sheet vanishes. (A similar retreat from a prograde slope was observed in stochastic simulations with the presence of peaks in the bed topography; Christian and others, 2022.) One might suppose that the disappearance of the ice sheet results from a negative surface elevation feedback in which the lowering of the ice-sheet surface results in the increased surface ablation that leads to further surface lowering and eventual contraction of the ice sheet. This feedback has been used to explain ice-sheet collapse under steady-state climate conditions (Garbe and others, 2020). However, a detailed examination of this simulation shows the collapse to be more complicated than a simple elevation-SMB feedback.

As Figure 4 illustrates, there is no simple connection between the loss of the ice-sheet surface area through which it gains mass and the disappearance of an ice sheet. For the simulations shown by the dark blue line in Figure 2a, the mass gain through the ice-sheet surface (blue line in Fig. 4a) for the most part exceeds the ice loss through the grounding line (orange line in Fig. 4) throughout the ice-sheet lifetime. It is only when ice advection toward the grounding line (orange line in Fig. 3b) significantly reduces that the ice sheet completely disappears. As Figures 3 and 4b illustrate, in the 2 ka period of the grounding-line advance and retreat (the green line in Fig. 3 shows the rate of the grounding-line migration) the ice flux through the grounding line (the orange line in Fig. 4b) does not change greatly; however, the integrated mass gain (the blue line in Fig. 4b) experiences significant variations in its magnitude and also sign. It is the rate of the ice-thickness change that balances these variations in the integrated mass gain (the red line in Fig. 4b). All three terms of the integrated mass balance have similar magnitudes and are equally important in determining both the instantaneous and long-term ice-sheet mass balance.

Figure 4. Time series of various terms of the integrated mass balance Eqn (21). Panel a shows the terms for the dark blue line in Figure 2a (the grounding line is on a prograde slope); panel b shows the terms during 2 ka period outlined by the dark rectangle in panel a.

The behavior of the grounding lines in other cases shown in Figure 2 can be understood using the same analysis described above for the case of Figure 2a. Ultimately, it is the imbalance between the advection of ice from upstream, the internal deformation and the accumulation at the grounding line, together with the geometric conditions at the grounding line, such as the bed slope and the ice-thickness gradient that determines whether the grounding line advances or retreats and at what rate.

3.2 Time-evolving basal conditions

Depending on the choice of parameters that determine temporal variability of basal sliding, Eqn (15), with all other parameters remained constant at their steady-state values, the ice sheets and their grounding lines exhibit a wide variety of behaviors including oscillation, retreat and advance. For example, we illustrate in Figure 5 evolutions with a 25 ka period of variability in the sliding parameter. In Figures 5a, b, the grounding line oscillates between limiting positions on the prograde and retrograde slopes with the same period regardless of the initial steady-state configuration. (Simulations (not shown) were run for 2 Ma to confirm the oscillatory behavior.) The bed slope alone is insufficient to explain this oscillatory behavior. While the grounding line indeed retreats from the retrograde slope, it continues to retreat on the prograde slope until it reaches its limiting position, and then re-advances far into the parts of the bed with retrograde slopes (Supplementary information, Movie 1).

Figure 5. Grounding-line response to time-variable sliding coefficients. (a–f) The grounding-line position x _g(t). All simulations were initialized with respective steady-state configurations and performed with the following parameters in Eqn (15): panels a–b: k _t = 2.8; x ₀ = 0.6 x _g, T = 25 ka; panel c: k _t = 9; x ₀ = 0.3 x _g, T = 400 years; panel d: k _t = 6; x ₀ = 0.59 x _g, T = 8.5 ka; panel e: k _t = −4, x ₀ = 0.2x _g, T = 20 ka; panel f: k _t = 3, x ₀ = 0.3x _g, T = 20 ka; in all simulations C ₀ = 7.6 × 10⁶ Pa m^−1/3 s^1/3, $\dot a = 0.1$ m a⁻¹. The blue rectangle in panel b marks the 25 ka interval shown in Figure 6. The red boxes outlining panels a and b indicate simulations that are repeated with the ice-flux parameterization and described in Section 3.3.

The grounding-line behavior during one cycle (marked by the thick blue line in Fig. 5b) is as follows. After advancing to its most downstream position, the grounding line rapidly retreats, and then slowly re-advances from its most upstream position. The two phases – retreat and advance – are not symmetric. The rate of the grounding-line retreat (Fig. 6, the green line) reaches its maximum magnitude ~690 m a⁻¹, and then slows down until it reaches its limiting upstream position. The magnitude of the rate of the grounding-line advance is an order of magnitude smaller than the magnitude of its retreat rate; its maximum ~70 m a⁻¹. The term $\displaystyle {-\dot a / D}$ is substantially smaller than the other two terms in Eqn (19). Consequently, the behavior of the grounding line (its advance and retreat) is primarily controlled by ice advection, deformation and changes in the ice thickness gradient caused by changes in the sliding conditions. The temporal evolution of the basal friction is such that the retreat from the most downstream position coincides with low basal shear near the grounding line, and the re-advance of the grounding line from its most upstream position coincides with the increase of the basal shear. Simulations with shorter periods and slightly different values of other parameters in Eqn (15) result in an unstopped retreat of the grounding lines starting from steady-state configurations on the prograde and retrograde parts of the bed (Figs 5c, d). As Figures 5e, f illustrate, the grounding lines can advance in an unstopped manner from the prograde and retrograde steady-state positions.

Figure 6. Individual terms contributing to the rate of the grounding-line migration in Eqn (19) for the simulation described by the thick blue line in Figure 5b (the grounding line oscillates between retrograde and prograde parts of the bed). Here, $\displaystyle {D = h_x + {b_x}/{( 1-\delta ) }}$.

In circumstances where $\dot a$ is constant, but the sliding properties vary in time, the temporal variability of the ice flux through the grounding line and the rate of the ice-thickness change integrated through the length of the ice sheet mimic each other (Fig. 7). Irrespective of the long-term behavior (i.e. either the grounding line exhibits regular oscillations shown in Figs 5a or 5b or unstopped advance shown in Fig. 5f), the terms of the integrated mass balance have similar magnitudes as shown in Figures 7a, b, respectively. The instantaneous balance of these terms is not informative about the long-term behavior of the ice sheet.

Figure 7. Time series of various terms of the integrated mass balance Eqn (21) for the simulations shown in Figure 5a (the grounding line oscillates between retrograde and prograde parts of the bed) and Figure 5f (the grounding line advances in an unstoppable manner).

3.3 Grounding-line behavior with the ice-flux parameterization

Under the assumptions of negligible bed slope b _x and the SMB $\dot a$ at the grounding line, Schoof (2007b) derived an expression for the ice flux (Eqn (16)) for the steady-state conditions. Due to its simplicity, it has been widely used in various applications in place of the exact description of the longitudinal stress at the grounding line. As we have noted, this expression equates the ice advection and ice deformation at the grounding line. However, it is the imbalance of these terms that contributes to the motion of the grounding line in Eqn (19), and it is not apparent to us that Eqn (16) is suitable in the time-variant case.

Previous studies (e.g. Gudmundsson, 2013; Reese and others, 2018) have demonstrated that this parameterization is not suitable for marine ice sheets whose ice shelves are laterally confined and experience buttressing. Here, we consider a configuration with unbuttressed ice shelves, for which expression (16) was derived. To assess its performance, we undertake simulations with the same time-variant SMB and basal sliding that resulted in the grounding-line behavior shown in respectively Figures 2a, d and 5a, b. The results of these simulations are shown in Figure 8. In general, we find that Eqn (16) results in dynamic evolutions that are markedly different in kind to those performed with the exact boundary conditions for the longitudinal stress at the grounding line. For example, in the simulations with time-variant SMB, which are shown in Figures 8a, b, the use of Eqn (16) (the red dashed lines) replaces either an irreversible retreat on the prograde slope with oscillatory behavior (Fig. 8a), or replaces on the retrograde slope an oscillating grounding-line behavior with an unstoppable advance (Fig. 8b). In simulations with the time-variant basal sliding, which are shown in Figures 8c, d, unstopped advances replace oscillatory behavior.

Figure 8. Effects of the ice-flux parameterization on the grounding-line migration. (a, b) Simulations with a time-variant SMB; (c, d) simulations with a time-variant sliding parameter. The red lines are simulations using Eqn (16); the green and blue lines are simulations with the exact treatment of the longitudinal stress at the grounding line (the lines are the same as in Figs 2a, d (marked by the red rectangles) and in Figs 5a, b). Simulations using Eqn (16) are truncated at the point when the ice sheet reaches the edge of the domain at 1000 km.

These markedly different evolutions arise because although the longitudinal stress condition relates the velocity gradient to the thickness at the grounding line, Eqn (16) insists that it is the ice flux that is determined solely by the thickness at the grounding line (i.e. no other characteristics such as the bed slope or the rate of the ice-thickness change affect it). To illustrate the difference that occurs in the grounding line flux, we compare the ice flux at the grounding line obtained with the exact treatment of the longitudinal stress, Eqn (18), to that computed with Eqn (16). As Figure 9 illustrates, these two fluxes are substantially different. In the case of the time-variant SMB (Fig. 9a), which corresponds to the grounding line exhibiting an unstopped retreat on the prograde slope (Fig. 2a), the ice flux computed with Eqn (16) both under- and over-estimates the simulated flux by factors ranging from four to more than ten. This is a result of the expression equating the ice advection and ice deformation at the grounding line. During the interval of grounding-line retreat in our simulations of Figure 2a, the SMB at the grounding line, and as a consequence the rate of the ice thickness change, experiences a broad range of values and cannot be neglected if one is to form a time-variant expression for the ice flux at the grounding line (Sergienko and Wingham, 2022). In the case of a time-variant sliding parameter, Eqn (16) under- and over-estimates the ice flux by ~30% (Fig. 9b). The discrepancy between the two fluxes is due to the contributions of the rate of the ice thickness change to the time-variant ice flux at the grounding line, and its dependence on the bed slope, whose effects become more pronounced for smaller values of the sliding parameter (Sergienko and Wingham, 2022, 2019).

Figure 9. Performance of ice-flux parameterization: the ratio of the ice flux computed in time-variant simulations with the exact treatment of the longitudinal stress at the grounding line to the ice flux computed with Eqn (16), vs ice thickness. (a) The case shown in Figure 2a (the grounding line oscillates between retrograde and prograde parts of the bed due to time-variant SMB); (b) the case shown in Figure 5b (the grounding line oscillates between retrograde and prograde parts of the bed due to time variant basal sliding).

4. Conclusions

Our results show that, once temporal variability of the external or internal conditions is accounted for, the same model (Schoof, 2007a, 2007b, 2012) that exhibits under constant conditions the irreversible retreat of the MISI hypothesis exhibits a diverse range of the grounding-line behavior – an unstoppable advance or retreat or irregular limited advance and retreat – regardless of the stability of a steady-state configuration achieved with constant conditions. Such behavior cannot be explained by a simple model of ice-sheet instability. This is because grounding-line migration is generally determined by the interplay between processes both at the grounding line and throughout the interior of the ice sheet, in addition to the geometric properties of the bed at the grounding line.

The model we employ is a very simple description of the ice dynamics: it lacks any description of lateral variability or lateral shear in either the sheet or the shelf, either of which may impact the dynamic behavior (e.g. Sergienko, 2012, 2022a; Gudmundsson, 2013; Schoof and others, 2017; Haseloff and Sergienko, 2018). Equally, it is the same model employed by Schoof (2007a, 2007b, 2012) to demonstrate instability in small perturbations from the steady state, and from which instability in more complex situations has been inferred. The models we used to capture the effects of time-variant SMB and time-variant basal conditions within the context of this simple model are asymmetric in their complexity, which is a reflection of our relative understanding of these processes. SMB is strongly dependent on temperature and contains considerable stochastic variability, and we have accommodated these effects within our model. Very little is known about the centennial to millennial variation in basal shear stress. We do not claim any particular virtue for our particular choice, beyond that it allows us to show the consequences on grounding-line migration that can emerge when the bed stress is time-variant. The detailed behaviors of the grounding line is sensitive to the choice of model parameters, particularly the sea-level temperature, but the variety of behaviors we illustrate is a common feature of the model. They reflect the variety of grounding-line behavior in the generally time-variant situation.

As the results of simulations with the time-variant SMB show, even if the grounding-line migration is caused by only stochastic variability in the climate conditions (here encapsulated in the variability of the SMB), this interplay can give rise to long-term trends in grounding-line behavior. Conversely, changes in the external conditions need not cause an immediate response of the grounding line, because other processes (deformation and sliding) also control its dynamics. These examples also illustrate that grounding-line migration depends on the history of changing environmental conditions, even if these changes are random in time. Consequently, the short-term grounding-line behavior (e.g. over several decades) may not indicate a response to the immediate environmental conditions; equally, it need not indicate a long-term behavior of the ice sheet and a grounding line. These results have direct implications for the interpretation of the behavior of the present-day ice sheets. The ice-sheet wide observations spanning the satellite era, which are a few decades long (Shepherd and others, 2018b), may be too short to make conclusive statements about the long-term behavior of their marine parts.

The results of simulations with the time-variant basal sliding illustrate that the grounding line can respond to changes in the basal conditions in the interior of ice sheets far away from the grounding line. Previous conceptual studies used similar mechanisms – changes in ice-sheet basal conditions – to explain the long-time variability of West Antarctic ice sheet (MacAyeal, 1992a, 1993). Although inferences about the spatial variability of the present-day basal conditions from surface observations have been performed routinely (MacAyeal, 1992b; Joughin and others, 2004; Sergienko and others, 2008; Brinkerhoff and others, 2021), nothing is known about their long-term temporal evolution. Current modeling projections of the future behavior of the present-day ice sheets are based on the assumption that basal conditions remain constant in time (Cornford and others, 2015; Seroussi and others, 2020). However, the results presented here illustrate that long-term changes in the basal conditions might cause an increase in the short-term (e.g. decadal) grounding-line migration rate that is an order of magnitude larger than the longer-term average. Thus, there is an urgent need to find ways to determine the temporal evolution of basal conditions in order to make reliable projections of the ice-sheet behavior in changing environmental conditions.

Our analysis of the integrated mass balance demonstrates that in time-variant conditions all its terms may have similar magnitudes and play an equal role in determining the behavior of the marine ice sheet. In circumstances where the surface accumulation varies in time, grounding-line retreat does not always lead to the reduction in the mass gain that happens under steady-state conditions. Additionally, in time-variant conditions, the rate of the ice-thickness change integrated through the horizontal extent of the ice sheet, which is zero in steady-state conditions, plays a significant role in the integrated ice-sheet mass balance.

We have also examined in the time-variant setting the use of the boundary condition due to Schoof (2007b), which equates the ice advection with the ice deformation at the grounding line. This is a reasonable approximation in the steady state (Sergienko and Wingham, 2022). However, in the general time-variant case, the grounding-line motion depends on small differences between the effects of advection and deformation. The ice flux computed in the time-variant simulations with the exact treatment of the longitudinal stress at the grounding line is substantially different from that obtained with this parameterization of the grounding-line ice flux in terms of the ice thickness. With an increasing number of climate models that use the large-scale ice-sheet models (Sadai and others, 2020; Pelletier and others, 2022; Park and others, 2023) it is necessary to recognize limitations of this ice-flux parameterization on the simulated behavior of marine ice sheets. In the time-variant case, the longitudinal stress at the grounding line requires a careful treatment.

Taking together, our results indicate that arguments and expressions developed for ice sheets in steady states are limited only to steady-state conditions. Studies of ice sheets experiencing temporally variable conditions require new, dedicated approaches.