the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The stability of presentday Antarctic grounding lines – Part 1: No indication of marine ice sheet instability in the current geometry
Emily A. Hill
Benoît Urruty
Ronja Reese
Julius Garbe
Olivier Gagliardini
Gaël Durand
Fabien GilletChaulet
G. Hilmar Gudmundsson
Ricarda Winkelmann
Mondher Chekki
David Chandler
Petra M. Langebroek
Theoretical and numerical work has shown that under certain circumstances grounding lines of marinetype ice sheets can enter phases of irreversible advance and retreat driven by the marine ice sheet instability (MISI). Instances of such irreversible retreat have been found in several simulations of the Antarctic Ice Sheet. However, it has not been assessed whether the Antarctic grounding lines are already undergoing MISI in their current position. Here, we conduct a systematic numerical stability analysis using three stateoftheart ice sheet models: Úa, Elmer/Ice, and the Parallel Ice Sheet Model (PISM). For the first two models, we construct steadystate initial configurations whereby the simulated grounding lines remain at the observed presentday positions through time. The third model, PISM, uses a spinup procedure and historical forcing such that its transient state is close to the observed one. To assess the stability of these simulated states, we apply shortterm perturbations to submarine melting. Our results show that the grounding lines around Antarctica migrate slightly away from their initial position while the perturbation is applied, and they revert once the perturbation is removed. This indicates that presentday retreat of Antarctic grounding lines is not yet irreversible or selfsustained. However, our accompanying paper (Part 2, Reese et al., 2023a) shows that if the grounding lines retreated further inland, under presentday climate forcing, it may lead to the eventual irreversible collapse of some marine regions of West Antarctica.
 Article
(5129 KB)  Fulltext XML

Supplement
(21610 KB)  BibTeX
 EndNote
Ice loss from the Antarctic Ice Sheet has increased in recent decades (Otosaka et al., 2023), and changes in ice discharge from Antarctica entail one of the greatest uncertainties in future projections of global sea level rise (Robel et al., 2019; Pattyn and Morlighem, 2020; IPCC, 2021). The grounding lines, i.e. the zones where the grounded ice sheet becomes so thin that it floats, are a key indicator of the health of the ice sheet. Accelerated retreat of the grounding lines could be an indication of the onset of a collapse of large marine regions of the ice sheet. Such a collapse (Mengel and Levermann, 2014; Feldmann and Levermann, 2015) could commit global sea levels to rise by several metres over the coming centuries to millennia (e.g. DeConto and Pollard, 2016; Golledge et al., 2015; Ritz et al., 2015; Cornford et al., 2015).
The term marine ice sheet instability (MISI) is often used to describe the potential for a selfreinforcing, positive feedback to be at play, where the retreat of the Antarctic grounding lines is internally driven (e.g. Pattyn and Morlighem, 2020). In particular, MISI is often discussed in relation to the susceptibility of the West Antarctic Ice Sheet to collapse due to this positive feedback mechanism. The existence of MISI means that a shift in the position of the grounding line can cause it to cross a critical threshold (or “tipping point”), beyond which the system is driven towards a different steady state. For marine, laterally uniform ice sheets with constant conditions, it has been shown that ice flux across the grounding line increases with thickness. Hence, retreat on a retrograde (inland) sloping bed into deeper water, and thus regions of greater ice thickness, promotes a positive feedback in which retreat continues, unabated, inland. In this case, no stable steadystate grounding lines exist on a retrograde sloping bed (Weertman, 1974; Schoof, 2007, 2012). However, in the case of laterally confined ice shelves that buttress the inland grounded ice, this feedback mechanism becomes more complex. Indeed, in the presence of buttressing ice shelves, stable steadystate groundingline positions can exist on a retrograde bed slope (Gudmundsson et al., 2012; Pegler, 2018; Haseloff and Sergienko, 2018). Most ice shelves around Antarctica provide such buttressing and have an important impact on inland ice dynamics (Fürst et al., 2016; Reese et al., 2018b). In addition to the ice shelf lateral confinement, nonnegligible bed topography found, for instance, under Thwaites Glacier and very weak beds, such as under Siple Coast ice streams, complicate stability conditions (Sergienko and Wingham, 2019, 2022).
Previous numerical modelling studies have shown the potential for tipping points related to MISI in the Antarctic Ice Sheet. For example Rosier et al. (2021) showed that three tipping points exists for the Pine Island Glacier. The retreat of Pine Island Glacier was found to be irreversible, because once initiated, reversing the perturbation to prethreshold conditions was not sufficient to halt or reverse it, and instead the forcing had to be reversed past its initial value to recover the initial state (Rosier et al., 2021). Furthermore, in their simulations of the whole Antarctic Ice Sheet, Garbe et al. (2020) found that retreat of West Antarctic grounding lines could be initiated by around 1–2 ^{∘}C of global warming above preindustrial, while the recovery of these grounding lines to their modern positions requires temperatures that are at least −1 ^{∘}C below the preindustrial average.
Previous studies have suggested that presentday retreat in regions of the West Antarctic Ice Sheet could mean that irreversible retreat has begun (Joughin et al., 2014; Rignot et al., 2014). However, to date there has not yet been a systematic analysis to assess whether irreversible retreat of Antarctic grounding lines is already underway. In this paper we use a systematic modelling approach to assess whether, under steady climate conditions, the groundingline positions of the Antarctic Ice Sheet are reversible with respect to a smallamplitude perturbation away from their current positions. Our modelling approach is outlined in detail at the beginning of Sect. 2. Briefly, we perform numerical experiments using three stateoftheart ice sheet models – Elmer/Ice (Gagliardini et al., 2013), Úa (Gudmundsson, 2020), and the Parallel Ice Sheet Model (PISM; Bueler and Brown, 2009; Winkelmann et al., 2011) – by applying a small but numerically significant perturbation to our initial model states, which all closely replicate the current geometry and velocity of the Antarctic Ice Sheet. If we find the grounding lines to either revert back to their former position (if the state is steady) or stay within the vicinity (if drifting through time), then this would indicate that current retreat of Antarctic grounding lines is unlikely to be due to an ongoing positive feedback mechanism, i.e. related to MISI.
A further question regarding the stability of the Antarctic grounding lines is the following: could the currently observed retreat, driven by presentday climate conditions, eventually commit the grounding lines to undergo irreversible retreat? This is specifically addressed in our accompanying paper (Part 2, Reese et al., 2023a), where longterm model simulations are used to assess whether presentday climate forcing has the potential to eventually lead to a collapse of major marine basins.
The paper is structured as follows: in the following section (Sect. 2) we present the initialization of our three ice sheet models and the perturbation experiments. The results are presented for the entire ice sheet and individual drainage basins in Sect. 3 and discussed further in Sect. 4.
To determine the stability regime of the Antarctic Ice Sheet in its current geometry, we use three different ice sheet models and two complementary methodologies. The ice sheet models used are Elmer/Ice, Úa, and PISM. Modelling experiments were conduced independently by three modelling teams, and each team tailored the details of the methodologies to the needs of their own models. The two methodologies applied are (i) stability analysis of steady ice sheet states and (ii) trajectory analysis of transient ice sheet states.
In the steadystate stability analysis, we construct a steady Antarctic Ice Sheet state (ice volume changes through time are approximately equal to zero) as close as possible to the currently observed ice sheet geometry and surface velocities. To determine whether this steady state of the ice sheet is stable or unstable, a smallamplitude perturbation was applied. If the ice sheet reverted back to the steady state after the perturbation was removed, it was considered to be stable with respect to the perturbation applied. This experimental procedure is motivated by the definition of asymptotic stability; i.e. a steadystate solution f_{e} is stable if there exists a δ>0 such that if $\Vert {f}_{\mathrm{e}}f\left(t\right)\Vert <\mathit{\delta}$, then f(t)→f_{e} when $t\to +\mathrm{\infty}$. This definition requires the construction of a steady state, i.e. f_{e}, to which the perturbation is applied, and an analysis of the convergence, or nonconvergence, over time of the perturbed state towards that steady state.
In the trajectory analysis, which is motivated by an orbital stability analysis, we perturbed a given state (steady or not) and analysed if the trajectory of the perturbed solution stayed within a given distance from the trajectory of the unperturbed solution. Formally, for the state f(t_{0}) and the resulting perturbed state g(t_{0}) for which $\Vert f\left({t}_{\mathrm{0}}\right)g\left({t}_{\mathrm{0}}\right)\Vert ={\mathit{\delta}}_{\mathrm{0}}$, we determined whether $\Vert f\left(t\right)g\left(t\right)\Vert \le {\mathit{\delta}}_{\mathrm{0}}$ for t>t_{0}. Hence, we observe whether the grounding line of the perturbed state stays with time at the same, or smaller, initial distance from the grounding line of the unperturbed state. Note that here we are only interested in assessing the local stability of the system around its current geometry and hence only observe the evolution of the state for a limited period after the applied perturbation.
An overview of our entire methodology is shown in Fig. 1. The steadystate stability experiments were performed using Elmer/Ice and Úa and can be summarized as follows:

The models are initialized using an inversion procedure to replicate the observed presentday geometry and ice velocities (Sect. 2.2.1).

A modification is made to the mass balance term such that if held constant in time, the presentday geometry is as close as numerically possible to a steady ice sheet state (Sect. 2.2.2).

To test whether this state is stable or unstable, the model states are perturbed (by increasing subshelf melt) for a short period of time to cause a small but numerically significant shift in the position of the grounding line (Sect. 2.4).

The perturbation is reversed, and we use the temporal evolution of the flux and groundingline positions to determine whether the constructed steady state is asymptotically stable (Sect. 3).
To ensure that our conclusions are not conditioned on the steadystate assumptions made in our first set of numerical experiments, we conducted a series of trajectory analysis experiments where the perturbation is applied to drifting grounding lines. In this case the model states are not in steady state. In PISM we initialize the ice sheet model as follows: (1) the model is initialized to preindustrial conditions using a spinup procedure, and (2) the model is run using historical forcing from ISMIP6 to obtain a transient state consistent with the observed presentday trend in mass loss (see Sect. 2.3 for further details). We then perturb the model states and determine whether the grounding lines of the perturbed and unperturbed transient model states remain close to each other or whether their trajectories start to diverge over time. As we are here interested in the stability of the grounding lines in their current position, we only follow the groundingline evolution for a limited period of time. It is possible that if we had followed the groundingline evolution for extended periods of time, the trajectories of the perturbed and unperturbed states would have started to diverge. However, that would no longer be a statement about the local stability regime of the grounding lines in the current ice sheet geometry.
2.1 Common approach
We first initialize all three models using as many common aspects of our models as possible to create initial states based on the observed geometry of the Antarctic Ice Sheet. These common inputs and approaches are summarized in Table 1. Namely, we use observations of bedrock topography, ice surface elevation, and ice thickness from BedMachine Antarctica (Version 2; Morlighem et al., 2020). To replicate the current ice flow, we take surface velocities from a recent snapshot in time (2015–2016) from the MEaSUREs Annual Ice Velocity Maps dataset (Version 1; Mouginot et al., 2019), which has a resolution of 1 km and good coverage across the entire ice sheet.
Across the surface of the ice sheet, all models initially apply a constantintime surface mass balance (later modified for Elmer/Ice and Úa), which is the output from the regional atmospheric climate model RACMO2.3p2 averaged from January 1995 to December 2014 (van Wessem et al., 2018). Melting at the bottom of ice shelves is a key control on the dynamics of the Antarctic Ice Sheet and is the focus of the perturbation experiments. All three models parameterize subshelf melt rates using their respective implementations of the Potsdam Ice shelf Cavity mOdel (PICO; Reese et al., 2018a). We ensure that the PICO geometry is the same as that in Reese et al. (2018a). PICO parameters were selected to reflect the sensitivity of subshelf melt rates in the Amundsen Sea Embayment (ASE) sector and Filchner–Ronne Ice Shelf to ocean temperature changes. Further details on this tuning of PICO parameters can be found in Reese et al. (2023a). This results in the parameter for vertical heat exchange ${\mathit{\gamma}}_{T}^{*}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ m s^{−1} and the parameter for the overturning strength C=2 Sv m^{3} kg^{−1}. While PICO is not a perfect representation of presentday melt rates, it can track the groundingline movement and provides melting rates for newly ungrounded regions.
(Joughin et al., 2019)(AsayDavis et al., 2016)(Schoof and Hindmarsh, 2010; Bueler and van Pelt, 2015)2.2 Steady states in Elmer/Ice and Úa
Elmer/Ice and Úa are both finite element models that solve the vertically integrated ice dynamics equations using the shallow shelf approximation (SSA; MacAyeal, 1989). Both models have been used extensively to solve ice flow problems in Antarctica and have participated in a number of model intercomparison experiments, e.g. Pattyn et al. (2012); Cornford et al. (2020). Both models first create an unstructured finite element mesh, where element sizes were refined in fastflowing regions and close to the grounding line. The initial ice sheet model states were then obtained using a twostep approach outlined in the following sections.
2.2.1 Model inversions
Using presentday observations of ice sheet geometry and velocity, the models are first initialized using a model inversion to estimate sliding law and rheology parameters using the adjoint method (MacAyeal, 1993). Both models impose a pressuredependent sliding law, which usually leads to more physically representative sliding close to the grounding line as compared to the Weertman sliding law (Brondex et al., 2019). The optimal fields for the basal slipperiness and viscosity parameters are found by minimizing a cost function, which is the sum of misfit and regularization terms. The main misfit term is the difference between observed and modelled velocities. Both models also apply an additional penalty on the rates of thickness change, to reduce nonphysical ice flux divergence anomalies. We regularize the inverse solutions using Tikhonov regularization terms that enforce smoothness of the inferred parameters (basal slipperiness and ice rate factor). The regularization weights are determined using an Lcurve analysis. By design, at the end of the inversion the surface velocities are in very good agreement with observations. See Appendix A1 and A2 for details on the inversion.
2.2.2 Mass balance modification
Observations clearly show that the presentday geometry of the ice sheet is not in steady state; i.e. ice thickness changes through time $\mathrm{d}h/\mathrm{d}t\ne \mathrm{0}$, and the current mass balance is $\mathrm{d}h/\mathrm{d}t={\dot{m}}_{\mathrm{pd}}\mathbf{\nabla}\cdot \left(\mathit{v}h\right)$ where v is the ice velocity and ${\dot{m}}_{\mathrm{pd}}$ is the total presentday mass balance (surface + base). However, given that stability is strictly speaking a property of steady states, in order to test the stability of the current grounding lines, we require a steadystate configuration of the ice sheet. Ideally, this steady state would be achieved through the inversion step. However, the penalty we apply to rates of thickness change is not sufficient to bring the model into a steady state, and when the model is run forward in time it drifts away from the presentday geometry and ice velocity. Hence, we create a modified mass balance term ${\dot{m}}_{\mathrm{mod}}$ such that $\mathrm{d}h/\mathrm{d}t=\mathrm{0}$ or as close as numerically possible. To do this, we follow a similar approach to others (Price et al., 2011; Goelzer et al., 2013) in which we subtract the modelled rates of ice thickness change ${\left(\mathrm{d}h/\mathrm{d}t\right}_{\mathrm{relax}}$ needed to keep the model in balance from ${\dot{m}}_{\mathrm{pd}}$ as follows:
Thereby, the presentday surface mass balance field ${\dot{b}}_{\mathrm{RACMO}}$ is the 1995–2014 averaged surface mass balance provided by RACMO2.3p2 (van Wessem et al., 2018) and ${\dot{b}}_{\mathrm{PICO}}$ the subshelf melt rates provided by PICO using 1975–2012 averaged ocean forcing from Schmidtko et al. (2014). By construction, the modification term ${\left(\right)}_{\mathrm{d}h/\mathrm{d}t}\mathrm{relax}$ reduces the mass balance to $\mathbf{\nabla}\cdot \left(\mathit{v}h\right)={\dot{m}}_{\mathrm{mod}}$ with $\mathrm{d}h/\mathrm{d}t\approx \mathrm{0}$. It is calculated in slightly different ways for the two models, the details of which are in Appendix A1 and A2. In summary, Elmer/Ice calculates ${\left(\mathrm{d}h/\mathrm{d}t\right}_{\mathrm{relax}}$ after a 20year relaxation period and then inputs this into a 1year forward relaxation, which forms the beginning of the perturbation experiments. Úa uses a semiiterative approach which takes ${\left(\mathrm{d}h/\mathrm{d}t\right}_{\mathrm{inv}}$ calculated after the inversion as input to Eq. (1) for a 1year relaxation period. The ${\dot{m}}_{\mathrm{mod}}$ field is then recalculated using ${\left(\right)}_{\mathrm{d}h/\mathrm{d}t}\mathrm{relax}$ at 1 year, and this is the beginning of the perturbation experiments. The modified mass balance fields and terms in Eq. (1) are shown in Figs. S1 to S3 in the Supplement.
We find that applying this modified mass balance term over the entire ice sheet (both grounded and floating portions) brings the ice sheet models to a steady state. However, by creating our modified mass balance field, we have created a surface mass balance field that deviates from climate conditions in RACMO. Broadly speaking, these fields represent the spatial pattern across the ice sheet; both models show higher positive mass balance around the coastal margins and more precipitation around the West Antarctic coast of the ASE and Bellingshausen Sea regions (Figs. S1 and S2). Additionally, the ${\left(\right)}_{\mathrm{d}h/\mathrm{d}t}\mathrm{relax}$ fields are regionally similar, where in both models the modification needed to obtain a steady state is largest in West Antarctica, in particular in the ASE sector. This reflects the current flux imbalance in West Antarctica and ongoing retreat and mass loss from this region. While the spatial patterns are broadly the same, e.g. large modifications across the ASE, we note that the distribution and amplitude of the modification to the mass balance vary greatly between the two models, likely due to subtle differences in our model initialization. In Úa the modification is higher in magnitude and locally spatially heterogeneous; i.e. both positive and negative values occur in close proximity to one another (see Fig. S3). By comparison, the modification in Elmer/Ice is smoother and generally smaller in magnitude. Crucially, it is assumed that modifying the mass balance in this manner has not altered the location of the critical thresholds in the groundingline position for the real ice sheet and that perturbing these steady states can allow us to learn about the stability of the presentday grounding lines in their current position. Given the differences between our two modified mass balance fields, and our results are consistent for both models (and PISM does not modify the mass balance; Sect. 2.3), our results are unlikely to be an artefact of this field having a specific shape. We also conducted a series of numerical experiments where we allowed the grounding lines to drift instead of forcing them to be stationary. This was done by applying either limited or no additional modification to the surface mass balance (Fig. S3). In all cases we found that once the perturbation was removed, grounding lines started to drift at the same rate as the (unperturbed) reference run. For unstable states one would expect the trajectories to diverge with time.
2.2.3 Initial states
The initial ice sheet model states created by Elmer/Ice and Úa closely replicate the observed ice thickness and surface velocities, particularly in the location of fastflowing ice streams (see Figs. S5 and S6). Importantly, the positions of the grounding lines in our models closely replicate the currently observed position of the grounding lines (Fig. 2) from BedMachine (Morlighem et al., 2020). In addition, the good agreement of observed and modelled groundingline positions for all models can be seen in our individual glacier profile figures presented in the results (Fig. 4). We calculate the error in the grounded ice area (a proxy for groundingline change), by differencing the simulated and observed grounded ice areas. To obtain a relative displacement of the grounding line itself, we normalize this area change by the simulated length of the grounding line. The resulting groundingline position errors are 84 m for Elmer/Ice and 540 m for Úa. During a steadystate control simulation for both models, there is less than 4.2 mm of change in sealevelequivalent volume over 100 years. Furthermore, there is very little drift in the position of the grounding line from the initial location, deviating by only 12.6 (Elmer/Ice) and 2.1 m (Úa) over 100 years, which is very small with respect to grid resolution (1 km close to the grounding line).
2.3 Transient state in PISM
To support the findings of this study further, we generate an initial state using PISM. Due to the spinup procedure of PISM, when it arrives at a steady ice sheet state, by definition that state will be stable. Hence, we would not need to additionally perturb such a state to determine whether it is stable or unstable. Instead, we conduct the trajectory analysis described above whereby we perturb a state including a presentday trend in mass changes. This allows us to examine whether current groundingline positions stay within the vicinity of the unperturbed control run. More details on the physics modelled in PISM and a comparison to the other models are given in Appendix A3.
2.3.1 Spinup to preindustrial conditions and historical run
The strategy adopted to build an initial state with PISM relies on spinup methods, with the approach taken here detailed in Reese et al. (2023a) and summarized in Fig. 1: we first create a thermal equilibrium with fixed geometry, followed by an ensemble of equilibrium simulations with full dynamics run under constant ∼ 1850 climate conditions (here equilibrium is defined as the integrated ice sheet mass balance being close to zero).
Starting from the 1850 initial configurations, historic simulations are run from 1850 to 2015. We use ISMIP6 historic forcing for the atmosphere and the ocean from 1850 to 2015. Climatologies for the preindustrial equilibrium state were created such that when adding the historic anomalies the atmosphere and ocean forcings between 1995 and 2014 match the presentday observations. Observed presentday velocities do not directly enter the initialization, but instead they are used, amongst other criteria, indirectly to determine optimal parameters in the initial state related to basal sliding and subshelf melting. We select the state that best replicates presentday ice thickness, groundingline position, mass loss, and ice surface velocities for the “best” PICO parameters (“AIS5” in Reese et al., 2023a).
2.3.2 Initial state
Despite the different initialization procedure of PISM, the resultant ice sheet reference state with respect to geometry and ice velocities is in good agreement with observations and the states obtained in the other models. PISM is capable of locating most, if not all, ice streams in their correct positions and accurately replicating the observed ice thickness and surface velocities after the spinup (see Figs. S5 and S6).
Due to the spinup procedure and the coarser grid resolution of PISM, there are some areas with greater deviation in the initial groundingline position. The overall error in the groundingline position for PISM is 11.4 km, which, while higher than the other two models, is very close to the horizontal grid size used in the model. Overall there is good agreement between modelled and observed groundingline positions (Fig. 2), including in regions of particular interest, such as Thwaites Glacier. The model groundingline positions are calculated using the flotation criterion.
During the historic simulation (1850–2015), the change in sealevelequivalent volume is −1.88 mm. While this is lower than observed rates of mass loss over the last decades, due to mass gains from snowfall, the pattern of thickness change is comparable to observations (see Reese et al., 2023a). During the control run the presentday climate conditions (2015) are held constant, so the climate forcing itself is steady, but the ice sheet state itself is not in a steady state (ice thickness changes through time are not equal to zero).
2.4 Experimental design
We have generated three ice sheet model initial states which closely replicate the current geometry, surface velocity, and groundingline positions of the Antarctic Ice Sheet (Table 2). We account for any remaining discrepancies between individual models by running control simulations alongside the perturbation experiments, and the results presented are then with respect to these control runs.
Morlighem et al. (2020)Morlighem et al. (2020)Morlighem et al. (2020)Morlighem et al. (2020)Morlighem et al. (2020)Morlighem et al. (2020)Gardner et al. (2018)van Wessem et al. (2018)During the control simulations the mass balance in Elmer/Ice and Úa is
where $\frac{\mathrm{d}h}{\mathrm{d}t}\approx \mathrm{0}$ and ${\dot{m}}_{\mathrm{mod}}$ are the same as defined for the previous models in Eq. (1). In PISM the mass balance is
and $\mathrm{d}h/\mathrm{d}t\ne \mathrm{0}$, broadly comparable to presentday rates of ice thickness change.
We apply a perturbation to the current position of the grounding lines by increasing the farfield ocean temperature that drives the melt rates calculated by PICO. Note that, in general, it would be possible to perturb the system using a number of different control parameters in our models. Given the important role that ice shelves play on inland ice sheet dynamics and groundingline position, via buttressing forces, we here choose to perturb the system by applying a shift in subshelf melt rates.
We increase the input ocean temperature, which is assumed to be representative of the conditions at depth on the continental shelf, by $\mathrm{\Delta}T\in \mathit{\{}\mathrm{1},\mathrm{3},\mathrm{5}\mathit{\}}$ ^{∘}C all around the Antarctic Ice Sheet. This perturbation in temperature is applied for 20 years to create a numerically significant groundingline retreat. By applying different increases in ocean temperature, we are able to test the robustness of this perturbation and found that a 5 ^{∘}C perturbation over 20 years remained small; i.e. it resulted in a small but obvious deviation in the position of the grounding line. While +5 ^{∘}C appears to be an unrealistically high magnitude change, we want to stress that this perturbation is not designed to be realistic and is only applied over a few decades. Instead, we can think of our small perturbation as a small movement of the grounding line away from its current position, on the order of a few grid cells or ice thicknesses at the grounding line. Indeed, the +5 ^{∘}C perturbation leads to approximately 10 km of retreat along the profiles in the ASE region (see Results) and along other profiles (see Figs. S16 and S17). We note that a recovery for the +5 ^{∘}C perturbation also implies that a smaller perturbation would also be reversible and indicates that the steady state is stable.
Across the floating ice shelves, during the 20year perturbation, the mass balance ${\dot{m}}_{\mathrm{\Delta}T}$ is the anomaly in subshelf melting calculated using PICO $\mathit{\delta}{\dot{b}}_{\mathrm{PICO}}$ applied to the modified mass balance fields in Elmer/Ice and Úa as
Across the grounded areas no perturbation is applied, and the mass balance is unchanged (Eq. 1).
In PISM the mass balance field has not been modified from presentday conditions, and unlike the other two models the initial temperatures T in PICO are perturbed directly by ΔT
We compare the integrated subshelf melt perturbation applied in all three models and find that it is comparable for all experiments (see Fig. S10b).
After 20 years we remove the perturbation and allow a recovery phase of the simulation for a further 80 years (Fig. 1). We extend some simulations by 400 years to test the robustness of the groundingline evolution over longer timescales. Throughout the forwardintime experiments, the extent of the ice shelves remains unchanged; however, all models impose a minimum ice shelf thickness, which is sufficiently thin to represent ice that has been removed and provides no buttressing.
During the recovery phase (80–480 years after the end of the perturbation), we examine the temporal evolution of the integrated ice flux across the grounding line. We choose groundingline flux as our metric of the system state because it is found to recover faster to a perturbation than grounded area or volume, due to the long timescales needed for the ice to thicken and readvance. In the steadystate initial states of Úa and Elmer/Ice, by design, the ice flux across the grounding line balances the surface accumulation upstream. Increased subshelf melt in our simulations leads to a sharp increase in ice flux across the grounding line, which is assumed to be due to a loss of buttressing as a result of ice shelf thinning. If the flux across the grounding line returns to its initial value after the perturbation is removed, this indicates that the ice sheet reverts to a steady state with a balance between surface accumulation in the grounded regions and groundingline flux (note that surface accumulation is altered a little by the groundingline movement). Hence, it is assumed that a return of the groundingline flux indicates the grounding line has either reverted back to its initial position or has begun to readvance towards its former position. When the grounding line does not retreat further, it means that it has found a new stable position very close of the previous one. If the flux were to increase away from its initial value, the grounding line is unstable. To support this we also examine the trend in groundedline position after the perturbation is removed, which is calculated as the change in grounded area for a constant groundingline length.
To further analyse the response of the grounding lines after the perturbation is removed, we calculate the efolding relaxation time, i.e. the time taken for the flux to decrease by a factor of e (Euler's number; ≈2.17). To do this we fit an exponential decay function in the form
to the change in flux during the recovery period of 80 years ΔQ(t), where ΔQ_{pert} is the change in flux at the end of the 20year perturbation relative to the initial, unperturbed flux. t is time after the perturbation, and τ is the recovery timescale. We repeat this for all three models and perturbation experiments, and these exponential curves can be seen in Figs. S7 to S9.
In the following sections we present the results of our perturbation experiments of presentday Antarctic grounding lines as described in Sect. 2.4. Figure 3 shows the integrated ice flux across the entire Antarctic grounding lines in each model during the perturbation experiments. We also present the results integrated across the 27 basins used in the Ice Sheet Mass Balance Intercomparison Exercise (IMBIE) (Fig. 2; Zwally et al., 2012), excluding basins 7, 8, and 25, which contain only small ice shelves. Additional figures showing groundingline position change and volume above flotation can be found in the Supplement (Figs. S10 to S15).
3.1 Antarctic Ice Sheet
On the Antarcticwide scale, all models and all perturbations show a similar trend: a strong increase in ice flux, reaching a maximum at the end of the 20year perturbation period, followed by an exponential decrease for the remaining 80 years (Fig. 3). The magnitude of the flux response to the melt perturbations is comparable between all ice sheet models, in particular for the 1 ^{∘}C temperature experiment, in which ice flux increased by approximately 300 Gt yr^{−1}. In the higher temperature scenarios (3 and 5 ^{∘}C), the flux responses diverge slightly from one another; PISM and Úa remain similar, while Elmer/Ice shows a stronger increase. This is likely due to subtle differences in the imposed basal melt perturbation in each model, which also diverge with the magnitude of the perturbation (see Fig. S10b).
The calculated efolding flux response time reveals that the Antarcticwide recovery time is in good agreement for all models, ranging between 10 and 20 years, and is largely independent of the magnitude of the perturbation (Fig. 3, bottom panel). During the 80year recovery period, the flux decreases rapidly but is not fully recovered. To check whether it is able to recover eventually, we extended the relaxation period in the 5 ^{∘}C simulations in all models to 480 years and found that the Antarcticwide flux returns to within 3.5 % of its original value (see Fig. S13). Alongside the ice flux evolution, the total retreat of the grounding line and ice volume show similar trends: rapid retreat and reduction in ice volume during the perturbation, after which retreat rates subside and grounding lines begin a slow recovery. Furthermore, for the Antarcticwide signal there is no indication of accelerated retreat (see Fig. S10). The recovery of the ice flux, alongside a short, 2decade efolding time, strongly indicates that the majority of Antarctic grounding lines are reversible. Also, ice flux evolution for individual basins (in all models) appears to show an exponential decrease after the perturbation is removed (Fig. 3). However, some basins recover more quickly than others. In the remainder of this section we explore the response of individual basins in more detail.
3.2 Amundsen Sea Embayment sector
The Amundsen Sea Embayment (ASE) sector encompasses the Pine Island, Thwaites, and Getz drainage basins (basin numbers 20 to 22 in Fig. 2). For all three of these basins, our results show a rapid exponential decay in the ice flux (efolding times <20 years) after the 20year perturbation (Fig. 3). We extract profiles along four marine glaciers in the ASE sector (Fig. 4) to show the response of a vertical transect along the ice shelves at four critical points of the 5 ^{∘}C experiment. These flow lines were directly interpolated from the original model grids. We note some differences in the bed geometry due to the different resolution and interpolation methods used by each model.
At Pine Island Glacier, the initial groundingline position in Úa is close to observations, whereas in Elmer/Ice and PISM the initial groundingline positions are located downstream at a topographic ridge (Fig. 4a). Grounding lines in Elmer/Ice and Úa retreat approximately 3 to 7 km (Fig. S12) across sections of retrograde sloping bed topography during the perturbation. In PISM the grounding line retreats along a section of prograde slope to the top of the ridge. For all models the 5 ^{∘}C perturbation causes the ice shelf to disappear entirely. After the perturbation is removed, grounding lines in all models advance back to their initial positions and the ice shelves regrow to their initial thickness. Groundingline flux (Fig. 3) increases during the 20year perturbation and sharply decreases thereafter. At year 50 the flux in Úa and PISM has decreased either below the control simulation (PISM) or within 0.7 % (Úa) of the initial flux. In Elmer/Ice the recovery is slower (11year efolding time) than the other models, and the flux remains 10 % higher than the initial flux at year 100.
Our results for Thwaites Glacier are similar to Pine Island (Fig. 4b). In the entire Thwaites basin (including Dotson and Crosson ice shelves) groundingline flux increases between 50 and 170 Gt yr^{−1} during the perturbation, and then it decreases rapidly (response times 8–15 years) in all models after the perturbation is removed. At Thwaites Glacier, the initial grounding lines in all three models are in close agreement with one another (Fig. 4b). For Elmer/Ice and Úa this is slightly downstream of the observed groundingline position at a topographic ridge, whereas PISM is located at the currently observed position. The retreat during the perturbation is not across a section of reverse bed slope, but instead it remains downstream of a second topographic ridge (located 70 km along the profile). For PISM and Úa, the ice shelf thins strongly during the perturbation, whereas in Elmer/Ice, the ice shelf is relatively thick in the initial state and does not disappear during perturbation. After the end of the perturbation, the ice shelf grows back to its initial shape in all models. In Úa the grounding lines readvance fully to their initial positions (within 100 years), and PISM readvances to the position arrived at in the control simulation after 500 years. The recovery in Elmer/Ice is slower than Úa, which corresponds to the longer recovery timescale for the entire Thwaites basin (Fig. 3).
Similar behaviour is also observed for Crosson and Dotson ice shelves for Elmer/Ice and Úa (Fig. 4c–d). The ice shelf thins and the grounding line retreats across a prograde slope during the perturbation. The grounding lines advance towards their initial position after the perturbation is removed, and in Úa it returns fully to its former position after 500 years. In contrast, the Crosson and Dotson grounding lines in PISM show signs of ongoing retreat. Crucially, the grounding lines in PISM start in slightly retreated positions compared to those of observations and the two other models. After the perturbation (20–100 years) the ice shelves thicken and the grounding lines readvance. However, after 100 years, the ice surface lowers and the grounding lines retreat to a location further inland of the perturbed position (at 20 years). Importantly, the groundingline positions at year 500 are very similar to their location in the control simulation, showing that they are reversible with respect to the trend in the control simulations. Longterm (10 000year) simulations under presentday climate forcing also show that the Crosson and Dotson grounding lines will remain (in the absence of additional forcing) at this slightly retreated position in the future (see Fig. 4 in Reese et al., 2023a).
For additional basins in this region of West Antarctica, the Getz (Fig. 3) and those in the Bellingshausen Sea sector (basins 23 and 24 in Fig. S11), the results are similar: the ice flux tends toward its initial value, the grounding lines readvance to their former positions (Fig. S12), and the response times are less than 20 years (Fig. 3). In summary, despite strong ice shelf thinning and retreat of the grounding lines in the ASE sector of the ice sheet during our perturbation experiments, we find that the grounding lines return to their initial positions after the perturbation is removed.
3.3 East Antarctica
Our experiments for basins in East Antarctica show a similarly strong signal of recovery to those in the ASE region. In particular, Wilkes, Totten, and Amery basins all display a rapid exponential decay in flux (efolding times <20 years) after the perturbation is removed (Fig. 3). The flux in these basins appears to be reversible in Úa and PISM within 100 years and within 500 years in Elmer/Ice (Fig. S15). Grounding lines in these basins also show a strong signal of readvance after the perturbation is removed (Fig. S12). This is also evident in selected glacier profiles at Totten and Cook glaciers (Fig. S17b–c). At both glaciers, the ice shelves disappear under the strongest perturbation (5 ^{∘}C), but at year 100 the ice shelves have almost recovered their former ice thickness and groundingline positions. At Cook Glacier the grounding lines have fully returned to their initial positions after 500 years in Úa and PISM, and in Elmer/Ice it tends towards the initial position. At Totten Glacier, all models start close to the observed groundingline position and retreat approx. 3 to 6 km inland (Fig. S12). After the perturbation is removed, Úa fully recovers by 500 years and Elmer/Ice has regrounded at the observed position (Fig. S17). In PISM the grounding line has not reverted to its initial position after 500 years but has instead approached the same location in the control simulation and can therefore still be considered reversible.
At Amery Ice Shelf, all models show that the ice shelves thicken, and the ice flux has almost fully recovered within 100 years, but the grounding line has not returned to its former position. However, it appears to have found a new stable groundingline position within a short (approx. <10 km) distance inland, and it does not experience accelerated retreat thereafter (Fig. S17a). We note that a number of smaller basins surrounding Totten, Wilkes, and Amery (basins 12, 15, and 16) show a similar recovery of the ice flux within approximately <20 years and readvance of the grounding lines.
In Dronning Maud Land (basins 4, 5, and 6), the signal of recovery is similar; all three ice sheet models show an increase in ice flux and retreat of the grounding line during the perturbation and a decreasing trend in ice flux after the perturbation is removed (Figs. S11 and S12). In general, the recovery of the groundingline position is slower in these basins (efolding times >15 years; Fig. 3). However, the extended relaxation period shows that the flux in the 5 ^{∘}C experiment recovers to within 3 % of its initial value for all models and all three basins (4–6). Concurrently, all grounding lines slowly readvance to their former positions and show no signs of accelerated retreat.
3.4 Filchner–Ronne and Ross ice shelf sectors
Elsewhere, the large basins draining into the Filchner–Ronne and Ross ice shelves show more complicated behaviour. In general, we do not see any strong signs of instability; i.e. the flux decays towards its initial value after the perturbation is removed, similar to the previously discussed basins. However, for the Filchner–Ronne basin in particular, the ice flux in all models remains 20 %–40 % away from its initial value after the 80year relaxation period in the 5 ^{∘}C experiment, and all models show response times of approximately 20 years (Fig. 3). In addition, all models show some signs of further retreat after the perturbation is removed (between 20–100 years in Fig. S12) but at a reduced rate. However, we find that in our extended simulations (for 5 ^{∘}C) the grounding lines either (1) tend to readvance towards their initial positions, (2) show a readvance of a few to several kilometres before appearing to reach a new steady position (Fig. S14), or (3) settle on a slightly retreated position a short distance inland. In no case do we see further retreat inland by the end of the simulations after 500 years. For both the Filchner–Ronne and Ross basins, the flux decreases to within 10 % (Elmer/Ice) or 6 % (PISM and Úa) of its initial value after 500 years.
We extract profiles at major outlets feeding into the large Filchner–Ronne and Ross ice shelves (see Fig. S16). For both profiles feeding the Ross Ice Shelf and the Ronne Ice Shelf, all models show signs of thinning and some additional retreat after the perturbation is removed (between 20 and 100 years). However, during the extended relaxation period (100 to 500 years), the ice shelves have begun to thicken again and the grounding lines do not retreat further inland. In the Filchner profile we find similar behaviour, with some additional retreat after the end of the perturbation (20 years). The grounding lines in PISM and Elmer/Ice do not retreat further after 100 years but remain in the same (slightly inland) position. In Úa the grounding line begins to readvance and has almost returned to its initial (and full) ice shelf extent after 500 years. Overall, there are no signs of accelerated retreat in these basins, and instead the grounding lines of the larger ice shelves may need longer to recover from a perturbation.
The perturbation experiments described above all showed that with time the grounding lines reverted back towards the state prior to the application of the perturbation. Where we ensured that the initial state was steady, as done in experiments using Elmer/Ice and Úa, the grounding lines migrated back to their original steadystate locations. Hence, these steady states are stable. In other experiments, where the grounding lines were allowed to migrate freely following the initialization procedure, the groundingline migration rates reverted back to those of the reference (i.e. unperturbed) runs. In no experiments conducted with the three icesheet models did we see any indication of selfenhanced irreversible retreat from the currently observed locations of the grounding lines of the Antarctic Ice Sheet.
Previous studies have suggested that the retrograde bed slopes of Pine Island and Thwaites glaciers may mean they could undergo phases of internally driven (MISI) retreat (e.g. Joughin et al., 2010; Rignot et al., 2014; Joughin et al., 2014; Favier et al., 2014; Mouginot et al., 2014; Seroussi et al., 2017; Milillo et al., 2022). Our results suggest that the current grounding lines in the ASE sector of West Antarctica are reversible in response to a small deviation from their current position. While previous work has indeed shown that the bed topography further inland could cause tipping points to be crossed at Pine Island Glacier (Rosier et al., 2021), our results show that it is unlikely that the current retreat of Pine Island is due to MISI. Importantly, once the grounding line retreats further towards the critical regions identified in Rosier et al. (2021), internal instability is likely to be found.
We find similar results for Thwaites Glacier where an entire removal of the ice shelf (in the strongest perturbation) is not sufficient to destabilize the current groundingline position, as has been suggested by previous studies (e.g. Wild et al., 2022). Several modelling studies have also shown that once the grounding lines retreat further inland under future increases in ocean forcing (and in the absence of any changes in surface mass balance from present day), it is possible that they will enter phases of accelerated retreat (Joughin et al., 2014; Seroussi et al., 2017). This accelerated retreat could be an indication of the onset of MISI. Note however that the projected groundingline evolution of Thwaites Glacier in response to future ocean changes is associated with large uncertainties (Nias et al., 2019; Robel et al., 2019). To assess whether presentday Antarctic grounding lines are committed to enter irreversible retreat under current climate forcing, we conduct multimillennialscale simulations using an ensemble of presentday configurations in PISM (including the state used here); see our accompanying paper (Reese et al., 2023a). Indeed, we find that current climate conditions can force grounding lines in the ASE sector to eventually enter irreversible retreat but not within 300 years of constant presentday climate. This supports previous suggestions that presentday climate conditions may be sufficient to trigger rapid groundingline retreat in West Antarctica in the long term (Golledge et al., 2021; Garbe et al., 2020; Joughin et al., 2014).
While our experiments, performed with three different ice sheet models, have shown consistent results, in the following paragraphs we note several caveats and potential sources of uncertainty. First, we parameterize the subshelf melt rates using the PICO subshelf melt model. While any parameterization of subshelf melt has associated uncertainties and may not capture certain physical processes in subshelf cavities, the subshelf melt distribution is not considered to greatly influence our results for two main reasons. First, our results are consistent across a wide range of temperature perturbations, including the strongest perturbation (+5 ^{∘}C) in which several ice shelves disappear. Second, our results are consistent across all three models despite the different implementations of PICO, which lead to subtly different melt rate distributions. In addition, we could have perturbed the grounding line using different model parameters, e.g. relating to ice viscosity or basal sliding. Indeed, we modified the basal sliding using Úa, and the results were consistent with those presented here; grounding lines readvanced after the perturbation was removed.
Secondly, the results of our experiments are dependent on accurate representations of the ice thickness and bedrock topography. One approach to assess the impact of ice thickness and bedrock topography on our results would be to conduct additional perturbation experiments with modifications to the topography within some bounds. In part, we have achieved this by using three different ice sheet models with different grid resolutions, where in particular the bed topography around the grounding lines in PISM is smoother than the other models (see Fig. 4). Grid resolution itself is a source of uncertainty in ice sheet model simulations, in particular around the grounding line (Durand et al., 2009; Pattyn et al., 2013), but again this is captured by our experiments with different models and grids. In addition, we repeated some experiments in Úa in which we modified the mesh resolution around the grounding line and found that our results were not affected by a finer resolution (500 m) at the grounding line.
An additional source of uncertainty in our results is the absence of dynamic calving front positions in our models. In Elmer/Ice and Úa the ice front is fixed and the ice shelves maintain a numerically required minimum thickness, due to the numerical challenges of removing ice entirely from the domain. It is possible that this thin layer of ice promotes ice shelf regrowth at a quicker rate than a true readvance of the ice front, due to the external forcing, surface, and basal mass balance, which are still applied across areas that have reached the minimum ice thickness. However, we note that the experiments conducted with PISM show similar recovery timescales, where in this case ice shelf cells are converted to ocean cells when the ice thickness decreases to zero. In addition, as the calving fronts are kept fixed, they are unable to advance beyond our prescribed initial position, which could also alter the stability of the grounding lines. We do not expect this to be relevant, because previous studies have shown that downstream regions (close to the calving fronts) of ice shelves provide little to no buttressing to upstream flow (Fürst et al., 2016; Reese et al., 2018b). We acknowledge that an evolving calving front may produce different results, especially as recent work has shown that (in the presence of buttressing) iceberg calving could impact the stability regime of grounding lines (Haseloff and Sergienko, 2022). Future studies should seek to include calving in such stability analyses.
Alongside calving, ice shelf damage can have an important impact on the magnitude of groundingline retreat due to weakening in the shear zones (Lhermitte et al., 2020). Hence, it is possible that, with damage included, the retreat of the grounding line in our perturbation experiments would be larger. However, damage is not yet sufficiently well parameterized for inclusion in our models, and implementing timedependent damage in particular is an ongoing (computational) challenge. An additional source of uncertainty is that we do not incorporate the surface melt elevation feedback during our simulations (Sergienko, 2022). Another point to make is that in our study here we only consider steady climate conditions in our experiments. Studies have shown that variability in the climate can cause noiseinduced tipping, which causes a system to transgress towards a qualitatively different state before the actual tipping point is crossed (e.g. Ashwin et al., 2012). Future work would benefit from incorporating timevarying climate conditions to explore the possibility of noiseinduced tipping in the real Antarctic Ice Sheet.
Our key finding is that the grounding lines of Antarctica, including Pine Island and Thwaites glaciers, show no indication of marine ice sheet instability (MISI) in their current state. We arrive at this conclusion by showing the existence of stable steadystate model configurations of the Antarctic Ice Sheet which closely resemble the observed presentday configuration as well as by showing reversibility of transient model configurations of the Antarctic Ice Sheet, under steady climate conditions. Our results indicate that MISI is not causing the presentday groundingline migration.
There is a general consensus within the ice sheet modelling community that the West Antarctic Ice Sheet is susceptible to MISI. Here, we have argued that the currently observed changes of the Antarctic Ice Sheet are not a manifestation of ongoing positive feedback related to MISI. While our experiments suggest that an internal instability threshold has not yet been crossed in Antarctica, future retreat driven by changes in climate conditions could force the grounding lines to cross a tipping point, after which retreat becomes driven by MISI. Further work is needed to quantify the amplitude and duration of forcing required for the Antarctic Ice Sheet to enter a phase of a largescale irreversible collapse involving groundingline retreat over hundreds of kilometres and a concomitant sea level rise of potentially several metres.
A1 Elmer/Ice
A1.1 Model description
Elmer/Ice (https://elmerice.elmerfem.org/, last access: 25 July 2023) is an opensource finite element software for ice sheets, glaciers, and ice flow modelling built on the multiphysics finite element model suite Elmer (Ruokolainen et al., 2023). Elmer/Ice is a very general and flexible tool and as such has been used for a large diversity of applications (181 publications since 2004). The main features and capabilities of Elmer/Ice have been described in Gagliardini et al. (2013) and in the associated publications (https://elmerice.elmerfem.org/publications, last access: 25 July 2023). Here, only the main characteristics relevant for our analysis will be presented. Regarding the physics included in Elmer/Ice, the ice flow velocity is computed solving the shallow shelf approximation (SSA) assuming an isotropic rheology following Glen's flow law. The initial viscosity field is computed using the 3D ice temperature field given by Van Liefferinge and Pattyn (2013) and using the values given in Cuffey and Paterson (2010) for the activation energies and prefactors used to compute the temperaturedependent rate factor A in Glen's flow law. This initial viscosity is then modified using inverse methods. Regarding the boundary conditions, the ice front of ice shelves is assumed to be fixed (i.e. the sum of calving flux and frontal melt flux equals the ice flux at the front). For the grounded part, the slidinglaw parameter field is inferred using inverse methods (see Sect. 2.2). For the floating part of the basal boundary, basal melt from the ocean is applied using the PICO model (Reese et al., 2018a), and no melt is applied to partially floating mesh elements. The grounding line is determined using a flotation criterion, and a subgrid scheme is applied for the basal drag over partially floating elements: SEP3 in Seroussi et al. (2014). Regarding the mesh, we use an anisotropic mesh adaptation scheme that uses the observed ice velocities and thickness. The mesh is preferentially refined along the directions of highest curvature of these two fields with an additional criterion function of the distance to the grounding line. The resulting mesh contains 545 837 nodes and 1 070 444 linear elements, and the size varies from 1 km close to the grounding line to 50 km in the interior. The mesh is held constant during the transient simulations.
Details concerning the model initialization, relaxation, and the basal sliding law used in the transient simulations are given below.
A1.2 Model inversion
To initialize the model, we use an inverse method to estimate model parameters that control the basal slipperiness and the ice rate factor by minimizing the misfit between observed (u_{obs}) and modelled (u_{mod}) velocities. We optimize β in a linear sliding law that relates the basal shear stress, t_{b}, to the basal sliding velocity, v_{b}:
This ensures that C_{eff}=10^{β} is positive. For the ice rheology, the vertically averaged effective viscosity used in the SSA is written as
where n is the Glen exponent, I_{D} is the second invariant of the strainrate tensor D defined as ${I}_{D}^{\mathrm{2}}=\mathrm{2}\mathit{D}:\mathit{D}$, and ${\overline{\mathit{\mu}}}_{\mathrm{ini}}$ is the initial vertically averaged ice rigidity computed using the 3D ice temperature field given by Van Liefferinge and Pattyn (2013) as explained above; to ensure that the viscosity remains positive, we optimize the parameter η starting from an initial guess of one (i.e. no modification from the viscosity initially predicted from the temperature field).
We apply a standard inverse methodology in which a cost function J(β,η), which is the sum of misfit (I) and regularization (R) terms, is minimized. The gradients of J with respect to β and η are determined in a computationally efficient way using the adjoint method. The misfit I and regularization R terms are defined as
As the velocity observation grids might be incomplete or have a better spatial resolution than the finite element mesh in the ice sheet interior, the difference between the model and the observations in I (Eq. A3) is evaluated at the N_{obs} locations where observations are present. The model velocities are interpolated using the natural finite element interpolation functions. For the regularization term R (Eq. A4), the first term computes the misfit between the modelled flux divergence ∇⋅(u_{mod}h) and the apparent point mass balance $\dot{a}$ integrated over the model domain Ω. Due to uncertainties in other model parameters that are not controlled during the inversion (e.g. the bed elevation), it is not possible in general to match both the velocities and the apparent point mass balance. So this term acts more as a regularization term that penalizes the highest ice flux divergence anomalies. The remaining anomalies are then dampened during a relaxation period (see next section). Here, for the apparent point mass balance we use the parameterization by DeConto and Pollard (2016) for the basal mass balance and RACMO for the surface mass balance and neglect the thickness rate of changes. The two other terms impose a smoothness constraint for the two inferred parameters β and η.
The regularization weights used in this study are ${\mathit{\lambda}}_{\mathrm{1}}=\mathrm{3.162}\times {\mathrm{10}}^{\mathrm{5}}$, ${\mathit{\lambda}}_{\mathrm{2}}=\mathrm{1.259}\times {\mathrm{10}}^{\mathrm{2}}$, and ${\mathit{\lambda}}_{\mathrm{3}}=\mathrm{7.943}\times {\mathrm{10}}^{\mathrm{4}}$. Following the principle of the Lcurve analysis, they have been empirically chosen from a large set of inversions, as those that give a good compromise with a low misfit and small regularization terms.
As the velocity dataset used for the common initial state is spatially incomplete, the inversion is first performed with a mosaic that aggregates observations from 2007 to 2018 (Mouginot et al., 2019) and thus has a nearly complete spatial coverage, but this comes at the expense of the accuracy in areas where velocities have largely changed. The minimization is then continued for 100 iterations using the 2015–2016 dataset to get closer to those observations while staying close to the initially inverted values in areas where observations are missing.
A1.3 Relaxation
The role of the relaxation is to reduce the inconsistency between input data and inverted data when we switch from a diagnostic to a prognostic simulation (GilletChaulet et al., 2012), which results in an unreasonably high surface elevation rate of change.
The model relaxation in Elmer/Ice is divided into three different steps. During the first 5 years, the relaxation is applied with the linear sliding law used for the inversion and the inverted friction coefficient. In floating areas, the friction parameter cannot be inverted, but it is set to a fixed value of $C=\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}\mathrm{a}$ to allow some friction if the grounding line were to advance. Then, to use a more realistic friction law, the regularized law from Joughin et al. (2019) (see section below) is applied for the following 15 years of relaxation. The conversion from the linear friction law to the regularized one is also described in Sect. A1.4. The last part of the relaxation is applying the mass balance modification as explained in Sect. 2.2. It consists of 1 year of balanced flux relaxation, in which the mass balance term is defined from Eq. (1) and the term ${\left(\right)}_{\frac{\mathrm{d}h}{\mathrm{d}t}}\mathrm{relax}$ is defined as the thickness rate of change from the last time step of the previous relaxation step. This step allows the model to reach a nearsteady state.
A1.4 Regularized Coulomb sliding law
For basal sliding, we adopt the regularized Coulomb sliding law proposed by Joughin et al. (2019):
which depends on the two parameters C_{s,m} and u_{0} and where
with h_{af} the height of ice above flotation and h_{T} a threshold height. Following Joughin et al. (2019), we adopt m=3, u_{0}=300 m a^{−1}, and h_{T}=75 m.
This sliding law (Eq. A5) exhibits two asymptotic behaviours: a Weertman regime ${\mathit{t}}_{\mathrm{b}}={C}_{\mathrm{s},m}/{u}_{\mathrm{0}}^{\mathrm{1}/m}\Vert {\mathit{v}}_{\mathrm{b}}{\Vert}^{\mathrm{1}/m\mathrm{1}}{\mathit{v}}_{\mathrm{b}}$ for $\Vert {\mathit{v}}_{\mathrm{b}}\Vert \ll {u}_{\mathrm{0}}$ and a Coulomb regime ${\mathit{t}}_{\mathrm{b}}=\mathit{\lambda}{C}_{\mathrm{s},m}{\mathit{v}}_{\mathrm{b}}/\Vert {\mathit{v}}_{\mathrm{b}}\Vert $ for $\Vert {\mathit{v}}_{\mathrm{b}}\Vert \gg {u}_{\mathrm{0}}$. It does not include a direct dependency on effective pressure N, whose role is subsumed in the model parameters. In pressuredependent sliding laws like that used in Úa (Eq. A10), u_{0} depends on the effective pressure. However, assuming a perfect hydrological connection between the subglacial drainage system and the ocean to compute N usually restricts the Coulomb regime to a small area close to the grounding line where the ice is close to flotation. Note that, in this particular case, we have $N=\mathit{\rho}g{h}_{\mathrm{af}}=\mathit{\rho}g{h}_{\mathrm{T}}\mathit{\lambda}$ such that both Eqs. (A5) and (A10) have the same dependency to water pressure in the vicinity of the grounding line. As the slidinglaw parameter C_{s,m} is determined through an inversion, it should include the dependency to N so that keeping it constant through time implicitly assumes that N does not change. Because this assumption is certainly not valid as the ice column approach flotation, the factor λ imposes a linear correction to the friction when h_{af} drops below the threshold h_{T} so that the friction decreases smoothly toward zero at the grounding line.
With the inversion being done assuming a linear Weertman sliding law (Eq. A1), C_{s,m} is inferred from the inverted effective slidinglaw coefficient C_{eff} such that
In floating areas, where C_{s,m} is not constrained by the inversion, we set a constant value of ${C}_{\mathrm{s},m}=\mathrm{10}$ kPa.
A2 Úa
A2.1 Model description
Úa is an opensource finite element ice flow model (Gudmundsson et al., 2012; Gudmundsson, 2020). The model is based on a vertically integrated formulation of the momentum equations and can be used to simulate the flow of large ice sheets such as the Antarctic and Greenland ice sheets, ice caps, and mountain glaciers. Úa solves the ice dynamics equations using the shallowshelf approximation (SSA) (MacAyeal, 1989) and Glen's flow law (Glen, 1955). The location of the grounding line is determined using the flotation condition. During forward transient experiments Úa allows for fully implicit time integration, and the nonlinear system is solved using the Newton–Raphson method. A minimum ice thickness constraint is used to ensure that ice thicknesses remain positive.
To initialize the Antarcticwide model we take the ice extent from BedMachine Antarctica (Morlighem et al., 2020), and within this boundary we generated a finite element mesh with 194 193 nodes and 385 097 linear elements using the Mesh2D Delaunaybased unstructured mesh generator (Engwirda, 2015). Element sizes were refined based on effective strain rates and distance to the grounding line and have a maximum size of 226 km in the very interior of the ice sheet, a mean size of 3.8 km, a median size of 1.57 km, and a minimum size of 0.68 km. Element sizes close to the grounding line are 1 km in size. We then linearly interpolated ice surface, thickness, and bed topography from BedMachine Antarctica (Morlighem et al., 2020) onto the model mesh. The surface boundary condition is stressfree, allowing the surface to respond freely to changes in surface velocity and surface mass balance. Surface mass balance is initially prescribed using output from RACMO2.3p2 (van Wessem et al., 2018), and subshelf melt is parameterized using an implementation of the PICO model (Reese et al., 2018a). Parameters of the basal slipperiness coefficient C in the sliding law and ice rate factor A in the flow law are determined using an inverse approach described in the following section.
A2.2 Model inversion
To initialize the model we use a data assimilation approach in which we estimate unknown parameters C and A by minimizing the misfit between observed (u_{obs}) and modelled (u_{mod}) velocities. Observed ice velocities (Mouginot et al., 2019) were linearly interpolated onto the model mesh. Úa uses a standard inverse methodology in which a cost function J, which is the sum of a misfit (I) and regularization (R) term, is minimized. The gradients of J with respect to A and C are determined in a computationally efficient way using the adjoint method and Tikhonovtype regularization. The misfit I and regularization R terms are defined as
where $\mathcal{A}=\int \mathrm{d}\mathcal{A}$ is the area of the model domain, ϵ_{obs} denotes measurement errors, and $\widehat{p}$ denotes the prior values for model parameters ($\widehat{A}$ and $\widehat{C}$). We use a uniform prior $\widehat{A}\approx \mathrm{1.15}\times {\mathrm{10}}^{\mathrm{8}}\phantom{\rule{0.125em}{0ex}}{\mathrm{kPa}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{yr}}^{\mathrm{1}}$ equivalent to an ice temperature of approx. −10 ^{∘}C using an Arrhenius temperature relation. For $\widehat{C}$ we estimate the prior as $\widehat{C}={u}_{\mathrm{obs}}/{\mathit{\tau}}_{\mathrm{b}}^{m}$ with and τ_{b}=80 kPa and m=3. The value of C beneath the ice shelves does not deviate from this prior value. Tikhonov regularization parameters γ_{s} and γ_{a} control the slope and amplitude of the gradients in A and C. Optimum values were determined using Lcurve analysis and are equal to γ_{s}=10 000 and γ_{a}=1. The inversion was run for 10 000 iterations, after which the cost function had converged.
Here we invert the model to estimate C using a commonly used sliding law (Eq. A10) that relates the basal traction t_{b} to the horizontal components of the bed tangential basal sliding velocity v_{b}. This was proposed by AsayDavis et al. (2016) (Eq. 11 in that paper) and used in a recent intercomparison experiment (Cornford et al., 2020). In our notation it reads
where μ_{k} is a coefficient of kinetic friction, and 𝒢 is the grounding–floating mask, with 𝒢=1, where the ice is grounded and 𝒢=0 otherwise. Here β^{2} is defined as
where C is the slipperiness coefficient, and we set m=3. When calculating the effective pressure, N, we assume a perfect hydrological connection with the ocean, i.e. $N=\mathcal{G}\phantom{\rule{0.125em}{0ex}}\mathit{\rho}g\phantom{\rule{0.125em}{0ex}}(h{h}_{\mathrm{f}})=\mathcal{G}\phantom{\rule{0.125em}{0ex}}\mathit{\rho}g\phantom{\rule{0.125em}{0ex}}{h}_{\mathrm{af}}$, where h_{f} is the flotation ice thickness (maximum ice thickness possible for a given ocean water column thickness, H, where $H=SB$ and S is the ocean surface and B is the bedrock) or ${h}_{\mathrm{f}}=H{\mathit{\rho}}_{\mathrm{w}}/\mathit{\rho}$. This ensures that the basal drag approaches zero as the grounding line is approached from above, i.e. t_{b}=0 at the grounding line.
The sliding law (Eq. A10) combines basal drag as calculated separately by the Coulomb and Weertman sliding laws through reciprocal weighting and, thus, represents a combination of those two sliding laws. Equation (A10) gives the limits:
A2.3 Relaxation
Similar to Elmer/Ice as described above (Sect. A1.3), Úa requires a period of relaxation after the inversion to dampen ice flux anomalies. Here, we use a twostep, semiiterative approach to apply a modification to the mass balance term in Eq. (1). First, we take rates of thickness change that are calculated at the end of the inversion ${\left(\mathrm{d}h/\mathrm{d}t\right}_{\mathrm{inv}}$ and use these as input to Eq. (1). Then, we relax the model forward in time for 1 year using this initial modification to $\dot{m}$, after which we take rates of thickness change ${\left(\mathrm{d}h/\mathrm{d}t\right}_{\mathrm{relax}}$ and use these as the second modification to the mass balance. This modified mass balance field is sufficient to bring the model into a steady state, in which ice volume changes are approximately zero. This modified mass balance field then remains fixed for all of the remaining control and perturbation simulations.
A3 PISM
We here extend on the model description of PISM. For more details on the spinup procedure, please see Reese et al. (2023a).
A3.1 Model description
The Parallel Ice Sheet Model (PISM; https://www.pism.io, last access: 25 July 2023; Bueler and Brown, 2009; Winkelmann et al., 2011) is an opensource ice dynamics model developed at the University of Alaska, Fairbanks, and the Potsdam Institute for Climate Impact Research. Ice flow velocities are obtained from a superposition of the SSA and the shallow ice approximation (SIA) velocity fields. PISM is thermomechanically coupled and solves the enthalpy evolution on a threedimensional grid. The ice rate factor is calculated taking into account the ice enthalpy and following the Glen–Paterson–Budd–Lliboutry–Duval flow law (Lliboutry and Duval, 1985), assuming a Glen exponent of n=3 and a SSA flow enhancement factor of E_{SSA}=1. The SIA flow enhancement factor is varied in the parameter ensemble. For the run presented in this article, we use a value of E_{SIA}=2.
Basal sliding is parameterized in the form of a generalized powerlaw formulation (Schoof and Hindmarsh, 2010):
where t_{b} is the basal shear stress, τ_{c} is the yield stress for basal till (see below), v_{b} is the SSA basal sliding velocity, and u_{0} is a threshold velocity. The slidinglaw exponent $q=\mathrm{1}/m$ can vary between 0 (purely plastic Coulomb sliding) and 1 (linear relationship between basal velocity and shear stress with coefficient ${\mathit{\tau}}_{c}/{u}_{\mathrm{0}}$). For the experiments presented here we adopt values of ${u}_{\mathrm{0}}=\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{yr}}^{\mathrm{1}}$ and q=0.75, respectively.
Basal shear stress in the vicinity of the grounding line is linearly interpolated on a subgrid scale between adjacent grounded and floating grid cells according to the height above buoyancy (Feldmann et al., 2014), which allows the grounding line to evolve freely. Note that subshelf melt is not interpolated across the grounding line and not applied in partially floating grid cells, as usually done in PISM. The yield stress τ_{c} is a function of parameterized till material properties (heuristic till friction angle ϕ) and effective till pressure N_{till} (“Mohr–Coulomb criterion”; Bueler and van Pelt, 2015):
where N_{till} is a function of the ice overburden pressure and the modelled effective amount of water in the till layer. No connection to the ocean is assumed in the calculation of the till water content; however, the till is assumed to be fully saturated when in contact with the ocean (in grid cells with floating ice or icefree ocean), which means that freshly grounded grid cells are usually slippery.
In the simulations presented here, for consistency with the other models, glacial isostatic adjustment is not considered, and we only calve floating ice that extends beyond the presentday extent of Antarctic ice shelves.
In contrast to Úa and Elmer/Ice, PISM simulations use a finite difference scheme to solve the momentum balance on a regular grid of 8 km horizontal resolution. While Seroussi et al. (2014) report that a horizontal resolution of 2 km is required to accurately represent groundingline dynamics, Feldmann et al. (2014) find that, by using a subgrid interpolation of friction, groundingline reversibility in PISM is also captured at coarser (Δx>10 km) resolutions. While a higher horizontal resolution is desirable, we here employ this interpolation to be able to run PISM over millennial timescales at 8 km horizontal resolution. Despite the abovementioned differences, we find that PISM results are in line with results from Elmer/ice and Úa that employ finer resolution around the grounding lines.
Elmer/Ice code is publicly available through GitHub (https://github.com/ElmerCSC/elmerfem, last access: 1 August 2023; Gagliardini et al., 2013) and an archived version of the Elmer code is available through Zenodo at https://doi.org/10.5281/zenodo.7892180 (Ruokolainen et al., 2023). All the simulations were performed with version 9.0 (Rev: 242e4bb) of Elmer/Ice. PISM code is publicly available at https://github.com/pism/pism (last access: 1 August 2023) and an archived version of the PISM code is available through Zenodo at https://doi.org/10.5281/zenodo.8101891 (Reese et al., 2023b). The Úa code is publicly available through GitHub (https://github.com/GHilmarG/UaSource/, last access: 25 July 2023), and an archived version of the model can be found at https://doi.org/10.5281/zenodo.3706623 (Gudmundsson, 2020).
Model output data are publicly available on Zenodo at https://doi.org/10.5281/zenodo.8086403 (Hill et al., 2023).
The supplement related to this article is available online at: https://doi.org/10.5194/tc1737392023supplement.
The experiments presented in this paper were collectively designed by members of the TiPACCs work package 2. BU, RR, and EAH set up and ran the simulations for Elmer/Ice, PISM, and Úa, respectively, and wrote the paper. JG created Figs. 1, 2, and 4. All authors contributed to the writing and discussion of ideas.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work is part of the TiPACCs project, which receives funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 820575. The Elmer/Ice computations presented in this paper were performed using HPC resources of CINES under the allocation 2021A0060106066 made by GENCI. We further acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research (BMBF), and Land Brandenburg for supporting this project by providing resources on the highperformance computer system at the Potsdam Institute for Climate Impact Research. The development of PISM is supported by NSF grant nos. PLR1644277 and PLR1914668 and NASA grant nos. NNX17AG65G and 20CRYO20200052. We acknowledge the use of the Northumbria University HPC facility Oswald to perform the Úa simulations. We like to thank Christian Schoof for providing helpful input on stability analysis. We thank the editor, Johannes Fürst, for handling our article and one anonymous reviewer and Alexander Robinson for their constructive comments that greatly improved the article.
This research has been supported by Horizon 2020 (TiPACCs (grant no. 820575)).
This paper was edited by Johannes J. Fürst and reviewed by Alexander Robinson and one anonymous referee.
AsayDavis, X. S., Cornford, S. L., Durand, G., GaltonFenzi, B. K., Gladstone, R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.: Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), Geosci. Model Dev., 9, 2471–2497, https://doi.org/10.5194/gmd924712016, 2016. a, b
Ashwin, P., Wieczorek, S., Vitolo, R., and Cox, P.: Tipping points in open systems: bifurcation, noiseinduced and ratedependent examples in the climate system, Philo. T. Roy. Soc. A, 370, 1166–1184, https://doi.org/10.1098/rsta.2011.0306, 2012. a
Brondex, J., GilletChaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177–195, https://doi.org/10.5194/tc131772019, 2019. a
Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res.Sol. Ea., 114, 1–21, https://doi.org/10.1029/2008JF001179, 2009. a, b
Bueler, E. and van Pelt, W.: Massconserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd816132015, 2015. a, b
Cornford, S. L., Martin, D. F., Payne, A. J., Ng, E. G., Le Brocq, A. M., Gladstone, R. M., Edwards, T. L., Shannon, S. R., Agosta, C., van den Broeke, M. R., Hellmer, H. H., Krinner, G., Ligtenberg, S. R. M., Timmermann, R., and Vaughan, D. G.: Centuryscale simulations of the response of the West Antarctic Ice Sheet to a warming climate, The Cryosphere, 9, 1579–1600, https://doi.org/10.5194/tc915792015, 2015. a
Cornford, S. L., Seroussi, H., AsayDavis, X. S., Gudmundsson, G. H., Arthern, R., Borstad, C., Christmann, J., Dias dos Santos, T., Feldmann, J., Goldberg, D., Hoffman, M. J., Humbert, A., Kleiner, T., Leguy, G., Lipscomb, W. H., Merino, N., Durand, G., Morlighem, M., Pollard, D., Rückamp, M., Williams, C. R., and Yu, H.: Results of the third Marine Ice Sheet Model Intercomparison Project (MISMIP+), The Cryosphere, 14, 2283–2301, https://doi.org/10.5194/tc1422832020, 2020. a, b
Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Academic Press, ISBN 9780123694614, 2010. a
DeConto, R. M. and Pollard, D.: Contribution of Antarctica to past and future sealevel rise, Nature, 531, 591–597, https://doi.org/10.1038/nature17145, 2016. a, b
Durand, G., Gagliardini, O., De Fleurian, B., Zwinger, T., and Le Meur, E.: Marine ice sheet dynamics: Hysteresis and neutral equilibrium, J. Geophys. Res.Earth, 114, F03009, https://doi.org/10.1029/2008JF001170, 2009. a
Engwirda, D.: Locally optimal Delaunayrefinement and optimisationbased mesh generation, PhD Thesis, School of Mathematics and Statistics, University of Sydney, https://ses.library.usyd.edu.au/handle/2123/13148 (last access: 25 July 2023), 2015. a
Favier, L., Durand, G., Cornford, S., Gudmundsson, G., Gagliardini, O., GilletChaulet, F., Zwinger, T., Payne, A., and Le Brocq, A.: Retreat of Pine Island Glacier controlled by marine icesheet instability, Nat. Clim. Change, 4, 117–121, https://doi.org/10.1038/nclimate2094, 2014. a
Feldmann, J. and Levermann, A.: Collapse of the West Antarctic Ice Sheet after local destabilization of the Amundsen Basin, P. Natl. Acad. Sci. USA, 112, 14191–14196, https://doi.org/10.1073/pnas.1512482112, 2015. a
Feldmann, J., Albrecht, T., Khroulev, C., Pattyn, F., and Levermann, A.: Resolutiondependent performance of grounding line motion in a shallow model compared with a fullStokes model according to the MISMIP3d intercomparison, J. Glaciol., 60, 353–360, https://doi.org/10.3189/2014JoG13J093, 2014. a, b
Fürst, J. J., Durand, G., GilletChaulet, F., Tavard, L., Rankl, M., Braun, M., and Gagliardini, O.: The safety band of Antarctic ice shelves, Nat. Clim. Change, 6, 479–482, https://doi.org/10.1038/nclimate2912, 2016. a, b
Gagliardini, O., Zwinger, T., GilletChaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a newgeneration ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd612992013, 2013. a, b, c
Garbe, J., Albrecht, T., Levermann, A., Donges, J. F., and Winkelmann, R.: The hysteresis of the Antarctic ice sheet, Nature, 585, 538–544, https://doi.org/10.1038/s4158602027275, 2020. a, b
Gardner, A. S., Moholdt, G., Scambos, T., Fahnstock, M., Ligtenberg, S., van den Broeke, M., and Nilsson, J.: Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years, The Cryosphere, 12, 521–547, https://doi.org/10.5194/tc125212018, 2018. a, b
GilletChaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sealevel rise from a newgeneration icesheet model, The Cryosphere, 6, 1561–1576, https://doi.org/10.5194/tc615612012, 2012. a
Glen, J. W.: The creep of polycrystalline ice, P. Roy. Soc. Lond. A, 228, 519–538, https://doi.org/10.1098/rspa.1955.0066, 1955. a
Goelzer, H., Huybrechts, P., Fürst, J. J., Nick, F. M., Andersen, M. L., Edwards, T. L., Fettweis, X., Payne, A. J., and Shannon, S.: Sensitivity of Greenland ice sheet projections to model formulations, J. Glaciol., 59, 733–749, https://doi.org/10.3189/2013JoG12J182, 2013. a
Golledge, N. R., Kowalewski, D. E., Naish, T. R., Levy, R. H., Fogwill, C. J., and Gasson, E. G.: The multimillennial Antarctic commitment to future sealevel rise, Nature, 526, 421–425, https://doi.org/10.1038/nature15706, 2015. a
Golledge, N. R., Clark, P. U., He, F., Dutton, A., Turney, C., Fogwill, C., Naish, T., Levy, R. H., McKay, R. M., Lowry, D. P., Bertler, N. A. N., Dunbar, G. B., and Carlson, A. E.: Retreat of the Antarctic Ice Sheet during the Last Interglaciation and implications for future change, Geophys. Res. Lett., 48, e2021GL094513, https://doi.org/10.1029/2021GL094513, 2021. a
Gudmundsson, G. H., Krug, J., Durand, G., Favier, L., and Gagliardini, O.: The stability of grounding lines on retrograde slopes, The Cryosphere, 6, 1497–1505, https://doi.org/10.5194/tc614972012, 2012. a, b
Gudmundsson, H.: GHilmarG/UaSource: Ua2019b, Zenodo [code], https://doi.org/10.5281/zenodo.3706624, 2020. a, b, c
Haseloff, M. and Sergienko, O. V.: The effect of buttressing on grounding line dynamics, J. Glaciol., 64, 417–431, https://doi.org/10.1017/jog.2018.30, 2018. a
Haseloff, M. and Sergienko, O. V.: Effects of calving and submarine melting on steady states and stability of buttressed marine ice sheets, J. Glaciol., 68, 1149–1166, https://doi.org/10.1017/jog.2022.29, 2022. a
Hill, E., Urruty, B., Reese, R., Garbe, J., Gagliardini, O., GilletChaulet, F., Gudmundsson, H., Winkelmann, R., Chekki, M., Chandler, D., and Langebroek, P.: Datasets of TC2022104, Zenodo [data set], https://doi.org/10.5281/zenodo.8086404, 2023. a
IPCC: Summary for Policymakers, in: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: MassonDelmotte, V., Zhai, P., Pirani, A., Connors, S., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J., Maycock, T., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, https://doi.org/10.1017/9781009157896.001, 2021. a
Joughin, I., Smith, B. E., and Holland, D. M.: Sensitivity of 21st century sea level to oceaninduced thinning of Pine Island Glacier, Antarctica, Geophys. Res. Lett., 37, L20502, https://doi.org/10.1029/2010GL044819, 2010. a
Joughin, I., Smith, B. E., and Medley, B.: Marine ice sheet collapse potentially under way for the Thwaites Glacier Basin, West Antarctica, Science, 344, 735–738, https://doi.org/10.1126/science.1249055, 2014. a, b, c, d
Joughin, I., Smith, B. E., and Schoof, C. G.: Regularized Coulomb friction laws for ice sheet sliding: Application to Pine Island Glacier, Antarctica, Geophys. Res. Lett., 46, 4764–4771, https://doi.org/10.1029/2019GL082526, 2019. a, b, c, d
Lhermitte, S., Sun, S., Shuman, C., Wouters, B., Pattyn, F., Wuite, J., Berthier, E., and Nagler, T.: Damage accelerates ice shelf instability and mass loss in Amundsen Sea Embayment, P. Natl. Acad. Sci. USA, 117, 24735–24741, https://doi.org/10.1073/pnas.1912890117, 2020. a
Lliboutry, L. and Duval, P.: Various isotropic and anisotropic ices found in glaciers and polar ice caps and their corresponding rheologies, Ann. Geophys., 3, 207–224, 1985. a
MacAyeal, D. R.: Largescale ice flow over a viscous basal sediment: theory and application to ice stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, https://doi.org/10.1029/jb094ib04p04071, 1989. a, b
MacAyeal, D. R.: A tutorial on the use of control methods in icesheet modeling, J. Glaciol., 39, 91–98, https://doi.org/10.3189/S0022143000015744, 1993. a
Mengel, M. and Levermann, A.: Ice plug prevents irreversible discharge from East Antarctica, Nat. Clim. Change, 4, 451–455, https://doi.org/10.1038/nclimate2226, 2014. a
Milillo, P., Rignot, E., Rizzoli, P., Scheuchl, B., Mouginot, J., BuesoBello, J., PratsIraola, P., and Dini, L.: Rapid glacier retreat rates observed in West Antarctica, Nat. Geosci., 15, 48–53, https://doi.org/10.1038/s4156102100877z, 2022. a
Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell, P., Goel, V., Greenbaum, J. S., Gudmundsson, H., Guo, J., Helm, V., Hofstede, C., Howat, I., Humbert, A., Jokat, W., Karlsson, N. B., Lee, W. S., Matsuoka, K., Millan, R., Mouginot, J., Paden, J., Pattyn, F., Roberts, J., Rosier, S., Ruppel, A., Seroussi, H., Smith, E. C., Steinhage, D., Sun, B., Broeke, M. R. van den, Ommen, T. D. van, Wessem, M. van, and Young, D. A.: Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nat. Geosci., 13, 132–137, https://doi.org/10.1038/s4156101905108, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m
Mouginot, J., Rignot, E., and Scheuchl, B.: Sustained increase in ice discharge from the Amundsen Sea Embayment, West Antarctica, from 1973 to 2013, Geophys. Res. Lett., 41, 1576–1584, https://doi.org/10.1002/2013GL059069, 2014. a
Mouginot, J., Rignot, E., and Scheuchl, B.: Continentwide, interferometric SAR phase, mapping of Antarctic ice velocity, Geophys. Res. Lett., 46, 9710–9718, https://doi.org/10.1029/2019GL083826, 2019. a, b, c
Nias, I. J., Cornford, S. L., Edwards, T. L., Gourmelen, N., and Payne, A. J.: Assessing uncertainty in the dynamical ice response to ocean warming in the Amundsen Sea Embayment, West Antarctica, Geophys. Res. Lett., 46, 11253–11260, https://doi.org/10.1029/2019GL084941, 2019. a
Otosaka, I. N., Shepherd, A., Ivins, E. R., Schlegel, N.J., Amory, C., van den Broeke, M. R., Horwath, M., Joughin, I., King, M. D., Krinner, G., Nowicki, S., Payne, A. J., Rignot, E., Scambos, T., Simon, K. M., Smith, B. E., Sørensen, L. S., Velicogna, I., Whitehouse, P. L., A, G., Agosta, C., Ahlstrøm, A. P., Blazquez, A., Colgan, W., Engdahl, M. E., Fettweis, X., Forsberg, R., Gallée, H., Gardner, A., Gilbert, L., Gourmelen, N., Groh, A., Gunter, B. C., Harig, C., Helm, V., Khan, S. A., Kittel, C., Konrad, H., Langen, P. L., Lecavalier, B. S., Liang, C.C., Loomis, B. D., McMillan, M., Melini, D., Mernild, S. H., Mottram, R., Mouginot, J., Nilsson, J., Noël, B., Pattle, M. E., Peltier, W. R., Pie, N., Roca, M., Sasgen, I., Save, H. V., Seo, K.W., Scheuchl, B., Schrama, E. J. O., Schröder, L., Simonsen, S. B., Slater, T., Spada, G., Sutterley, T. C., Vishwakarma, B. D., van Wessem, J. M., Wiese, D., van der Wal, W., and Wouters, B.: Mass balance of the Greenland and Antarctic ice sheets from 1992 to 2020, Earth Syst. Sci. Data, 15, 1597–1616, https://doi.org/10.5194/essd1515972023, 2023. a
Pattyn, F. and Morlighem, M.: The uncertain future of the Antarctic Ice Sheet, Science, 367, 1331–1335, https://doi.org/10.1126/science.aaz5487, 2020. a, b
Pattyn, F., Schoof, C., Perichon, L., Hindmarsh, R. C. A., Bueler, E., de Fleurian, B., Durand, G., Gagliardini, O., Gladstone, R., Goldberg, D., Gudmundsson, G. H., Huybrechts, P., Lee, V., Nick, F. M., Payne, A. J., Pollard, D., Rybak, O., Saito, F., and Vieli, A.: Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP, The Cryosphere, 6, 573–588, https://doi.org/10.5194/tc65732012, 2012. a
Pattyn, F., Perichon, L., Durand, G., Favier, L., Gagliardini, O., Hindmarsh, R. C. A., Zwinger, T., Albrecht, T., Cornford, S., Docquier, D., Fürst, J. J., Goldberg, D., Gudmundsson, G. H., Humbert, A., Hütten, M., Huybrechts, P., Jouvet, G., Kleiner, T., Larour, E., Martin, D., Morlighem, M., Payne, A. J., Pollard, D., Rückamp, M., Rybak, O., Seroussi, H., Thoma, M., and Wilkens, N.: Grounding Groundingline migration in planview marine icesheet models: results of the ice2sea MISMIP3d intercomparison, J. Glaciol., 59, 410–422, https://doi.org/10.3189/2013JoG12J129, 2013. a
Pegler, S. S.: Suppression of marine ice sheet instability, J. Fluid Mech., 857, 648–680, https://doi.org/10.1017/jfm.2018.742, 2018. a
Price, S. F., Payne, A. J., Howat, I. M., and Smith, B. E.: Committed sealevel rise for the next century from Greenland ice sheet dynamics during the past decade, P. Natl. Acad. Sci. USA, 108, 8978–8983, https://doi.org/10.1073/pnas.1017313108, 2011. a
Reese, R., Albrecht, T., Mengel, M., AsayDavis, X., and Winkelmann, R.: Antarctic subshelf melt rates via PICO, The Cryosphere, 12, 1969–1985, https://doi.org/10.5194/tc1219692018, 2018a. a, b, c, d
Reese, R., Gudmundsson, G. H., Levermann, A., and Winkelmann, R.: The far reach of iceshelf thinning in Antarctica, Nat. Clim. Change, 8, 53–57, https://doi.org/10.1038/s415580170020x, 2018b. a, b
Reese, R., Garbe, J., Hill, E. A., Urruty, B., Naughten, K. A., Gagliardini, O., Durand, G., GilletChaulet, F., Gudmundsson, G. H., Chandler, D., Langebroek, P. M., and Winkelmann, R.: The stability of presentday Antarctic grounding lines – Part 2: Onset of irreversible retreat of Amundsen Sea glaciers under current climate on centennial timescales cannot be excluded, The Cryosphere, 17, 3761–3783, https://doi.org/10.5194/tc1737612023, 2023a. a, b, c, d, e, f, g, h, i, j
Reese, R., Garbe, J., Hill, E., Urruty, B., Naughten, K., Gagliardini, O., Durand, G., GilletChaulet, F., Gudmundsson, G. H., Chandler, D., Langebroek, P. M., Winkelmann, R., and other PISM authors: Data and code for publication “The stability of presentday Antarctic grounding lines – Part B”, Zenodo [code and data set], https://doi.org/10.5281/zenodo.8101891, 2023b. a
Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H., and Scheuchl, B.: Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011, Geophys. Res. Lett., 41, 3502–3509, https://doi.org/10.1002/2014GL060140, 2014. a, b
Ritz, C., Edwards, T. L., Durand, G., Payne, A. J., Peyaud, V., and Hindmarsh, R. C.: Potential sealevel rise from Antarctic icesheet instability constrained by observations, Nature, 528, 115–118, https://doi.org/10.1038/nature16147, 2015. a
Robel, A. A., Seroussi, H., and Roe, G. H.: Marine ice sheet instability amplifies and skews uncertainty in projections of future sealevel rise, P. Natl. Acad. Sci. USA, 116, 14887–14892, https://doi.org/10.1073/pnas.1904822116, 2019. a, b
Rosier, S. H. R., Reese, R., Donges, J. F., De Rydt, J., Gudmundsson, G. H., and Winkelmann, R.: The tipping points and early warning indicators for Pine Island Glacier, West Antarctica, The Cryosphere, 15, 1501–1516, https://doi.org/10.5194/tc1515012021, 2021. a, b, c, d
Ruokolainen, J., Malinen, M., Råback, P., Zwinger, T., Takala, F., Kataja, J., GilletChaulet, F., Ilvonen, S., Gladstone, R., Byckling, M., Chekki, M., Gong, C., Ponomarev, P., van Dongen, E., Robertsen, F., Wheel, I., Cook, S., t7saeki, luzpaz, and Rich_B.: ElmerCSC/elmerfem: Elmer 9.0 (release9.0), [Code] Zenodo, https://doi.org/10.5281/zenodo.7892181, 2023. a, b
Schmidtko, S., Heywood, K. J., Thompson, A. F., and Aoki, S.: Multidecadal warming of Antarctic waters, Science, 346, 1227–1231, https://doi.org/10.1126/science.1256117, 2014. a
Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res.Earth, 112, F03S28, https://doi.org/10.1029/2006JF000664, 2007. a
Schoof, C.: Marine ice sheet stability, J. Fluid Mech., 698, 62–72, https://doi.org/10.1017/jfm.2012.43, 2012. a
Schoof, C. and Hindmarsh, R. C.: Thinfilm flows with wall slip: an asymptotic analysis of higher order glacier flow models, Q. J. Mech. Appl. Math., 63, 73–114, https://doi.org/10.1093/qjmam/hbp025, 2010. a, b
Sergienko, O. and Wingham, D.: Grounding line stability in a regime of low driving and basal stresses, J. Glaciol., 65, 833–849, https://doi.org/10.1017/jog.2019.53, 2019. a
Sergienko, O. V.: No general stability conditions for marine icesheet grounding lines in the presence of feedbacks, Nat. Commun., 13, 1–6, https://doi.org/10.1038/s41467022298923, 2022. a
Sergienko, O. V. and Wingham, D. J.: Bed topography and marine icesheet stability, J. Glaciol., 68, 124–138, https://doi.org/10.1017/jog.2021.79, 2022. a
Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087, https://doi.org/10.5194/tc820752014, 2014. a, b
Seroussi, H., Nakayama, Y., Larour, E., Menemenlis, D., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, 6191–6199, https://doi.org/10.1002/2017GL072910, 2017. a, b
Van Liefferinge, B. and Pattyn, F.: Using iceflow models to evaluate potential sites of million yearold ice in Antarctica, Clim. Past, 9, 2335–2345, https://doi.org/10.5194/cp923352013, 2013. a, b
van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Lenaerts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016), The Cryosphere, 12, 1479–1498, https://doi.org/10.5194/tc1214792018, 2018. a, b, c, d, e
Weertman, J.: Stability of the junction of an ice sheet and an ice shelf, J. Glaciol., 13, 3–11, https://doi.org/10.3189/S0022143000023327, 1974. a
Wild, C. T., Alley, K. E., Muto, A., Truffer, M., Scambos, T. A., and Pettit, E. C.: Weakening of the pinning point buttressing Thwaites Glacier, West Antarctica, The Cryosphere, 16, 397–417, https://doi.org/10.5194/tc163972022, 2022. a
Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISMPIK) – Part 1: Model description, The Cryosphere, 5, 715–726, https://doi.org/10.5194/tc57152011, 2011. a, b
Zwally, H. J., Giovinetto, M. B., Beckley, M. A., and Saba, J. L.: Antarctic and Greenland drainage systems, GSFC cryospheric sciences laboratory, https://earth.gsfc.nasa.gov/cryo/data/polaraltimetry/antarcticandgreenlanddrainagesystems (last access: 25 July 2023), 2012. a, b