Present-day mass loss rates are a precursor for West Antarctic Ice Sheet collapse

Observations of recent mass loss rates of the West Antarctic Ice Sheet (WAIS) raise concerns about its stability since a collapse would increase global sea levels by several meters. Future projections of these mass loss trends are often estimated using numerical ice sheet models, and recent studies have highlighted the need for models to be benchmarked against present-day observed mass change rates. Here, we present an improved initialization method that optimizes local agreement not only with observations of ice thickness and surface velocity but also with satellite-based estimates of mass change rates. This is achieved by a combination of tuned thermal forcing under the floating ice shelves and friction under the ice sheet. Starting from this improved present-day state, we generate an ensemble of future simulations of Antarctic mass change by varying model physical choices and parameter values while fixing the climate forcing at present-day values. The dynamical response shows slow grounding-line retreat over several centuries, followed by a phase of rapid mass loss over about 200 years with a consistent rate of ∼3 mm GMSL yr⁻¹ (global mean sea level). We find that, for all ensemble members, the Thwaites Glacier and Pine Island Glacier collapse. Our results imply that present-day ocean thermal forcing, if held constant over multiple centuries, may be sufficient to deglaciate large parts of the WAIS, raising global mean sea level by at least a meter.

1 Introduction

The projected Antarctic contribution to global mean sea level (GMSL) rise ranges from 0.03 m (SSP1-1.9 (Shared Socioeconomic Pathway), low end of the likely range) to 0.34 m (SSP5-8.5, high end of the likely range) in 2100 (Fox-Kemper et al., 2021). Dynamical processes are not expected to accelerate or decelerate ice mass change significantly before 2100 (van de Wal et al., 2022). After 2100, however, the GMSL driven by nonlinear processes is highly uncertain, possibly accelerating GMSL rise significantly (Fox-Kemper et al., 2021; Payne et al., 2021).

One such process is marine ice sheet instability (MISI; see Schoof, 2012, 2007, and Durand et al., 2009), which could drive a mostly dynamical irreversible retreat (i.e., the glacier will not return to its present grounding-line position when the forcing is removed) of marine-terminating glaciers. Robel et al. (2019) argued that the mere possibility of MISI amplifies uncertainty in sea level rise projections from ice sheet models. MISI can occur for parts of the ice sheet where outlet glaciers terminate on a retrograde-sloping, submerged bedrock, due to a positive feedback between the flux across the grounding line (GL; the transition line between a grounded ice sheet and a floating ice shelf) and the ice thickness at the GL. Thwaites Glacier (TG) and Pine Island Glacier (PIG), in the Amundsen Sea region of the West Antarctic Ice Sheet (WAIS), are located in such regions and may therefore be vulnerable to MISI.

Currently, the Thwaites GL retreats by up to 1.2 km yr⁻¹ on the eastern shelf and up to 0.9 km yr⁻¹ along the western margin (Milillo et al., 2019), and the GL of PIG retreated 31 km in 20 years (Rignot et al., 2014). Warm Circumpolar Deep Water (CDW) and related high basal melt rates have been observed under Pine Island Glacier (Jacobs et al., 2011a) since the early 1990s, and studies show that CDW intrusions have coincided observed with PIG speedup (Thoma et al., 2008; Jenkins et al., 2016). As a result, the Amundsen Sea Embayment (ASE; which contains TG and PIG) is currently the largest contributor to Antarctic Ice Sheet (AIS) mass loss (Smith et al., 2020; Rignot et al., 2019). Further retreat of PIG and TG would increase ice drainage from these basins into the ocean, eventually triggering the collapse (i.e., relatively fast, potentially irreversible ungrounding and ice loss) of a significant part of the WAIS over several centuries.

Several modeling studies have assessed the potential for ASE collapse. One study (Joughin et al., 2014) argued that under present-day melt rates, TG might already be on a trajectory toward accelerated retreat; moderate retreat in this century will likely be followed by a phase of rapid collapse beginning in the next 200 to 900 years. Another study (Favier et al., 2014) used three ice sheet models to show that PIG is now undergoing a climate-driven 40 km retreat, but they made no projections on longer timescales. Both studies mention MISI as the main driver of retreat, but the retreat might be reversed by sufficient ocean cooling (Favier et al., 2014). A number of recent studies suggest that TG and PIG are unstable under the current climate and could collapse on a timescale up to 2000 years (Garbe et al., 2020; Golledge et al., 2021; Reese et al., 2023; Coulon et al., 2024). Lipscomb et al. (2021) projected accelerated retreat leading to a collapse only when ocean thermal forcing increases by 1–2 K relative to a preindustrial or 20th-century equilibrium. However, there are few historical observations of Southern Ocean temperatures, so it is unknown how much the ASE has warmed in the past century and if the warming needed to drive such a retreat has already happened.

Studying ASE glacier stability requires century-scale model simulations, which contain significant uncertainties. Seroussi et al. (2024) found that potential AIS contributions to GMSL rise dramatically after 2100, with large variations across ice sheet models. According to Aschwanden et al. (2021), uncertainties in ice sheet modeling arise from four main sources: suboptimal model initialization, incomplete understanding of physical processes, numerical model uncertainty, and uncertainty in the climate forcing. Some ice sheet models are initialized using a long spin-up to the steady state. Since they are initialized to a stable ice sheet without model drift, they do not represent the observed mass change in recent decades. A good representation of recent mass changes is essential for reliable projections, as these changes are likely to continue or accelerate.

Broadly, two methods are used to initialize ice sheet models: data assimilation and spin-up. In the first, data assimilation methods (Larour et al., 2012; Gudmundsson et al., 2012; Cornford et al., 2015; Hoffman et al., 2018; Bradley and Arthern, 2021) capture conditions at a certain time. For example, based on the observed ice thickness, the ice surface velocities are calculated and compared to observations. Uncertain parameters are tuned to minimize a cost function based on the difference between observed and modeled quantities until the velocities converge to a state close to observations. This has the advantage that the resulting ice sheet matches the observed thickness and is ideally close to observed velocities at the start of a forward run, without necessarily being in a steady state. However, the resulting dynamic state might not replicate the observed mass change rates. The simulation of recent mass change can be improved where annual or subannual observations are available by doing a transient calibration as suggested by Goldberg et al. (2015). Also, the modeled mass change rates can be added to the cost function, as was done for example by Rosier et al. (2021) to minimize model drift. Alternatively, Bett et al. (2024) used the ice sheet model WAVI (Arthern et al., 2015) and added the difference between observed and modeled mass change rates to the cost function, using observed mass change rates from Smith et al. (2020). The WAVI inversion method yields spatially varying values of two free parameters relating to basal friction and ice viscosity. Since both ice velocities and mass change rates are targeted simultaneously through the cost function, it is likely that there is a trade-off between errors in these quantities and that the modeled mass change rates do not agree with observations in all locations. Neither Arthern et al. (2015) nor Bett et al. (2024) give a quantitative comparison of modeled and observed mass change rates. More recently, Rosier et al. (2024) used the misfit between observed and modeled mass change rates in their data assimilation inversion, similar to Rosier et al. (2021). They state explicitly that it is not possible to obtain a perfect fit for both ice surface velocities and observed mass change rates, since those two datasets are not mutually consistent.

The other initialization approach is the spin-up method (Winkelmann et al., 2011; Pollard and DeConto, 2012; Greve and Blatter, 2016; Quiquet et al., 2018; Lipscomb et al., 2019; Berends et al., 2021, 2022), which consists of a long run (typically for several thousand years) with paleoclimatic or preindustrial forcing during which the ice sheet can freely evolve. Optionally, uncertain parameter values can be tuned to nudge variables like ice thickness and surface velocities toward present-day observed values. This happens in the model runtime: at regular intervals, one or more free parameters are adjusted as necessary to decrease the misfit between the model and observations. The resulting ice sheet is ideally close to observations and in near equilibrium with little drift (i.e., with little mass change). The modeled ice sheet can also be tuned to historical estimates or observations and then advanced to the present day using forcing from global climate models (GCMs) to avoid an unrealistic steady state at the start of a forward run (Reese et al., 2023; Coulon et al., 2024). This has the advantage of obtaining the modeled mass change rates through a physical process (ocean warming or changes in the surface mass balance, SMB), but this adds two uncertainties: the ice sheet geometry from before the observational period and the forcing over the historical period.

This study describes an addition to the spin-up method to match observed mass change rates at the start of forward simulations without relying on uncertain historical ice sheet geometry and forcing. This method is similar to the mass balance modification mentioned in Hill et al. (2023), but instead of minimizing mass change rates, we aim to match the observed rates. These rates are derived from satellite-observed surface elevation changes (Smith et al., 2020), corrected for firn processes and glacial isostasy. We apply this method to the Community Ice Sheet Model (CISM; Lipscomb et al., 2021, 2019), allowing us to initialize the Antarctic Ice Sheet in its observed disequilibrium. We use this setup to conduct forward simulations without changing the climate forcing. In these simulations, we find that Thwaites Glacier, Pine Island Glacier, and eventually large parts of the WAIS collapse over many centuries. We run similar simulations using the Utrecht Finite Volume Ice-Sheet Model (UFEMISM) (Berends et al., 2021) to show that the method can be easily implemented across models and to make sure that WAIS collapse is not unique to CISM.

Section 2 describes the two ice sheet models, the initialization method, and the subsequent forward runs. In Sect. 3, we show that the new initialization procedure results in a modeled present-day state in CISM which is in very good agreement with observations (ice thickness RMSE of 34 m, ice surface velocity RMSE of 121.2 m yr⁻¹, GL position on average within 1.5 km from the observed GL) and importantly has the same disequilibrium as observed by Smith et al. (2020). In Sect. 4, we discuss the default forward run, starting from the initialized state, under sustained present-day atmosphere and (inverted) ocean forcing, purposely excluding further climate forcing. We repeat the default simulation with UFEMISM (Sect. 5) to show that the observed retreat of TG and PIG is consistent in the two models. In Sect. 6, we show that the CISM behavior is robust under a wide variety of physical approximations and parameter settings. The paper ends with a discussion in Sect. 7 and a conclusion in Sect. 8.

2 Methods

2.1 CISM

The Community Ice Sheet Model (CISM; Lipscomb et al., 2021, 2019) is a thermo-mechanical higher-order ice sheet model, which is part of the Community Earth System Model version 2 (CESM2; Danabasoglu et al., 2020). This study primarily uses CISM v2.1 and builds on earlier applications of CISM to Antarctic Ice Sheet retreat (Lipscomb et al., 2021; Berdahl et al., 2023). CISM is run at 4 km resolution with a vertically integrated higher-order approximation of the momentum balance (Goldberg, 2011; Robinson et al., 2022) and a regularized Coulomb sliding law (Zoet and Iverson, 2020) as the basal boundary condition. We apply a surface mass balance and surface air temperature climatology spanning 1979–2016 from RACMO version 2.3p2 (van Wessem et al., 2018). Basal melt rates are calculated using a quadratic relation with a thermal forcing observational dataset (Jourdain et al., 2020). Internal ice temperatures are allowed to evolve freely during the initialization and the forward runs. The calving front is maintained at its present-day position. This conservative approach has been chosen because calving physics and the impact of iceberg melange of ice shelves are still poorly understood, and overly retreating calving fronts enhance mass loss, which makes the results overly negative. Variable names and parameter values used in the following text are summarized in Tables 1 and 2.

Table 1 Variables used in this study.

Table 2 Parameters and their units and values used in this study.

2.1.1 Basal friction

In the default configuration, basal friction is parameterized using a basal sliding law based on laboratory experiments (Zoet and Iverson, 2020). This sliding law combines elements of ice sliding on a hard bed (e.g., bedrock) with sliding on deformable soft beds (e.g., saturated till). Its end cases are Weertman-style sliding, with basal velocities dependent on some power of the basal friction, and Coulomb sliding, with basal friction and velocities decoupled. The former mechanism usually dominates where the ice sheet is flowing on hard, non-deformable beds, while the latter describes fast-flowing outlet glaciers with softer deformable till such as Pine Island Glacier and Thwaites Glacier. In CISM, the basal sliding law is implemented as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{b}}=C_{\mathrm{c}}N{\left({\displaystyle \frac{u_{\mathrm{b}}}{u_{\mathrm{b}}+u_{\mathrm{0}}}}\right)}^{\frac{\mathrm{1}}{m}},\end{array}$$

where C_c is a unitless parameter controlling the strength of the Coulomb sliding. In CISM, C_c corresponds to tan ϕ of Zoet and Iverson (2020), in which ϕ is the friction angle, a material property of the subglacial till. Since C_c is poorly constrained by theory and observations, we use it as a spatially variable tuning parameter, which we tune in the following way using a nudging method (Lipscomb et al., 2021):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{c}}}{\mathrm{d}t}}=-C_{\mathrm{c}}\left[\left({\displaystyle \frac{H-H_{\mathrm{obs}}}{H_{\mathrm{0}}\mathit{\tau }}}\right)+{\displaystyle \frac{\mathrm{2}}{H_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}+{\displaystyle \frac{r}{\mathit{\tau }}}\ln {\displaystyle \frac{C_{\mathrm{c}}}{C_{\mathrm{r}}}}\right],\end{array}$$

in which H is the modeled ice thickness, H_obs is the observed ice thickness, and C_r is the relaxation target for C_c. H₀, τ, and r are empirical constants used to vary the relative magnitude of each term. The third term in brackets, which is not included in Lipscomb et al. (2021), provides a relaxation target C_r that helps prevent C_c from drifting toward extreme values (limited to 10⁻⁶ and 1). It is a function of elevation, with lower values at low elevation where soft marine sediments are likely more prevalent. We chose targets of 0.1 for bedrock below −700 m a.s.l. and 0.4 for above 700 m a.s.l., with a linear interpolation in between, similar to Aschwanden et al. (2013). Our values for H₀ and τ are 100 m and 100 years, respectively.

The effective pressure N in Eq. (1) is the pressure at the ice–bed interface, equal to the difference between the ice overburden pressure and the subglacial water pressure. The effective pressure is lowered near GLs to represent the connection of the subglacial hydrological system to the ocean (Leguy et al., 2014):

$$\begin{array}{}\text{(3)}& N={\mathit{\rho }}_{\mathrm{i}}gH{\left(\mathrm{1}-{\displaystyle \frac{H_{\mathrm{f}}}{H}}\right)}^p,\end{array}$$

where ρ_i is the density of ice; g is gravitational acceleration; H is the ice thickness; and p is an exponent between 0 and 1, with a larger value implying a stronger ocean connection. A value of p=0.5 is chosen in this study unless stated otherwise to include some hydrological connection (Leguy et al., 2021; Lipscomb et al., 2021). Different values are tested in the sensitivity analysis. The floatation thickness H_f is the height of an ice column resting on bedrock below sea level (b<0) at hydrostatic equilibrium:

$$\begin{array}{}\text{(4)}& H_{\mathrm{f}}=max\left(\mathrm{0},-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}b\right),\end{array}$$

where ρ_w is the density of seawater and b is the bedrock elevation.

2.1.2 Basal melt rates

Beneath grounded ice, CISM computes basal melt rates based on the net of geothermal, frictional, and conductive heat fluxes at the bed. In our runs, however, there is no coupling between basal water and sliding parameters. Beneath floating ice, we use the ISMIP6 local parameterization (Jourdain et al., 2020):

$$\begin{array}{}\text{(5)}& b_{\mathrm{mlt}}={\mathit{\gamma }}_{\mathrm{0}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}{c}_{p\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}}}\right)}^{\mathrm{2}}{\left(max\left[{\text{TF}}_{\mathrm{base}}+\mathit{\delta }T,\mathrm{0}\right]\right)}^{\mathrm{2}}.\end{array}$$

The max function ensures non-negative melt rates; negative rates imply refreezing (accretion), which is not included. We tune δT to optimize the ice thickness agreement with observations, similar to Eq. (2):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}\left(\mathit{\delta }T\right)}{\mathrm{d}t}}=-\mathit{\delta }T\left[\left({\displaystyle \frac{H-H_{\mathrm{obs}}}{H_{\mathrm{0}}{\mathit{\tau }}_{\mathrm{o}}}}\right)+{\displaystyle \frac{\mathrm{2}}{H_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}\right]+{\displaystyle \frac{\left(T_{\mathrm{r}}-\mathit{\delta }T\right)}{{\mathit{\tau }}_{\mathrm{o}}}},\end{array}$$

where T_r is a relaxation target for the temperature correction. Ideally, the correction is small, and therefore we set T_r=0 everywhere. Melt rates are sensitive to γ₀ in Eq. (5). We chose a default value of ${\mathit{\gamma }}_{\mathrm{0}}=\mathrm{3.0}\times {\mathrm{10}}^{\mathrm{4}}$ m yr⁻¹, for which the average δT in the TG and PIG basins at the end of the transient spin-up is close to 0. In other words, we assume that the observed ocean temperatures in this region are of the right magnitude to drive basal melt rates consistent with observations. Our value for γ₀ lies about halfway between the median MeanAnt and PIGL calibration values suggested by Jourdain et al. (2020). The depth dependency is introduced in the dataset TF_base from Jourdain et al. (2020) in two ways. First, the thermal forcing increases with depth due to the decrease in the pressure melting point. Second, in the dataset of Jourdain et al. (2020), there is warm CDW present in the ASE, leading to increased thermal forcing with depth. In a forward run, we keep δT constant at the value computed during the spin-up. When the grounding line retreats, we extrapolate basin-average values of δT (basins are defined following Zwally et al., 2015) to newly floating cells.

2.1.3 Grounding-line subgrid-scale parameterization

The GL is not explicitly modeled in CISM but follows from the hydrostatic balance. The modeled GL cuts through cells, with some cells partly grounded and partly floating. To prevent abrupt jumps in modeled quantities at the GL, we use a GL parameterization (Leguy et al., 2021). We define a floatation function as

$$\begin{array}{}\text{(7)}& f_{\mathrm{float}}=-b-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}H,\end{array}$$

which varies smoothly from negative for grounded ice to positive for floating ice. (For floating ice, f_float is the water column thickness in the cavity.) f_float is used to compute the floating fraction as a percentage of the grid cell area by bilinearly interpolating its value from cell vertices to the cell areas scaled to the cavity thickness (Leguy et al., 2021). The grounded and floating fraction of a cell are then used to scale basal friction and basal melting. For the basal melting, this is referred to as the partial melt parameterization (PMP). Other parameterizations tested in the sensitivity experiments of this study are the floatation condition melt parameterization (FCMP), in which full melt is assigned to a cell once the cell center is floating, and the no-melt parameterization (NMP), in which no melting is allowed in any cell that is partly grounded.

2.1.4 Calving

There are several calving laws in the literature (e.g., Yu et al., 2019; Wilner et al., 2023; Greene et al., 2022). However, there is no agreed-upon best approach to Antarctic calving (this holds for the Greenland Ice Sheet as well; see for example Benn et al., 2017), and most calving laws struggle to reproduce the observed calving front at multiple locations simultaneously without adjusting local parameters (Amaral et al., 2020). We therefore choose to apply a no-advance calving scheme, preventing the calving front from advancing beyond the observed present-day location. The calving front can retreat but only when the ice thins below a threshold thickness of 1 m. In most places, it does not move. This is a conservative approach for a retreating system; in practice, we would expect the calving front to move inland and therefore provide less buttressing for upstream ice. By taking a conservative approach, we rule out the possibility that a simplistic calving algorithm is (partly) responsible for the modeled retreat.

2.2 Initialization methods

Typically, CISM nudges basal friction parameters and ocean temperatures to optimize the agreement of ice thickness with observations, allowing a modest misfit (a misfit scale of H₀=100 m in Eqs. 2 and 6) to avoid overfitting but without considering the observed mass change rates. We refer to this method as the “equilibrium spin-up”. This results in a steady state, where the integrated SMB is in equilibrium with the integrated ice losses: the modeled calving and basal melt fluxes. There is no isostatic adjustment in the default experiments; we assume a static bedrock. The model converges within ∼10⁴ model years, which is typically assessed by considering several criteria like change in ice mass ice area and ice velocities. A converged model state should not drift significantly when run forward with fixed inverted quantities. We calculate the model drift in the modeled ice sheet after stopping the inversion (i.e., by continuing the simulation with the inverted fields kept constant) for 2000 years. We accept an initialization as “stable” once there is little to no instantaneous model drift and there is no TG and PIG collapse after 2000 years. Then, it is certain that a PIG and TG collapse is caused by applying present-day mass change rates and not influenced by model drift from the initialization.

The addition to the equilibrium spin-up, referred to as the transient spin-up, allows for the model to reach a geometry in agreement with observations, similar to the result of the equilibrium initialization, while also allowing for the ice to thin or thicken in agreement with satellite observations at the start of the forward run (i.e., at the present day). Our goal is to create the observed dynamic disequilibrium, so we want to include the observed mass change rates as a target. To hold the ice geometry near observed values during the spin-up, we apply a correction term equal to the negative of the observed mass change rates. That is, we artificially apply a thickening term to ice that is thinning and a thinning term to ice that is thickening. This procedure trains the model to produce mass change rates that are (nearly) identical to the observations during the subsequent forward run, when the correction term is removed.

The mass balances for the equilibrium initialization and the transient spin-up, respectively, are given by

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathit{\delta }H}{\mathit{\delta }t}}=-\mathrm{\nabla }F+B,\text{(9)}& {\displaystyle \frac{\mathit{\delta }H}{\mathit{\delta }t}}=-\mathrm{\nabla }F+B+C.\end{array}$$

In these equations, ∇F is the ice flux divergence, B is the sum of the basal and surface mass balances, and C is the correction term applied to hold the ice at the observed thickness during the spin-up. Ideally, the initialization procedure ends when the left-hand sides of Eqs. (8) and (9), $\frac{\mathit{\delta }H}{\mathit{\delta }t}$, are 0. Hence, at the end of the equilibrium initialization, which uses Eq. (8), ∇F is being balanced by B, and forward simulations start with $\frac{\mathit{\delta }H}{\mathit{\delta }t}\approx \mathrm{0}$. In Eq. (9), used for transient initialization, we end the initialization when $\frac{\mathit{\delta }H}{\mathit{\delta }t}\approx \mathrm{0}$. When we remove the C term, forward simulations start with $\frac{\mathit{\delta }H}{\mathit{\delta }t}\approx -C$, where the right-hand side is equal to the observed mass change rate, as desired.

The inversion methods are stopped when the modeled ice sheet has little to no model drift once run forward with the inverted friction and ocean temperatures and particularly no significant ice mass loss in the ASE basin. The modeled ice sheet should ideally have little instantaneous drift compared to the observed mass change rates and no thinning or collapsing glaciers during a continuation run of 2000 years with fixed inverted quantities and no mass change rates applied.

During a 2000-year model drift experiment, the model tends to grow to a slightly advanced state for PIG and the Dotson Ice Shelf (see Figs. S1, S2, and S3). This advance is mostly due to the local positive feedback between the ice shelf geometry and basal melt rates that amplifies any small residual drift. Model drift is largely reduced in a forward run with a time constant instead of geometry-based (e.g., Eq. 5) basal melt rate field (not shown). Initialization takes about 10⁴ model years.

2.3 Forward simulations

Forward simulations are continued without additional climate forcing to assess the impact of the present-day ocean thermal forcing on the evolution of the AIS and thereby the sea level contribution over time. The future simulation using the transient spin-up method is referred to from now on as the “default experiment”. As the thermal forcing at a given ocean depth does not change, the modeled future evolution of the ice sheet is the response to the present-day forcing.

In forward simulations, the ice sheet starts to retreat after the mass correction term in the spin-up is removed. This has implications for the inverted C_c and ocean temperatures in cells that change from grounded to floating or vice versa. If grid points initialized as grounded ice start floating, we extrapolate the basin-average ocean temperature corrections under the shelf from the end of the initialization. If initially floating grid points become grounded, we use an elevation-dependent parameterization to estimate the local C_c, similar to what was done in Aschwanden et al. (2013) and the same as what was used to calculate the relaxation target C_r.

The default forward simulation is compared against an extensive set of sensitivity experiments with modified physics and parameter choices. We designed these experiments conservatively, with many model choices that would tend to stabilize the modeled ice sheet and could potentially stop PIG and TG from collapsing. Section 6 and the Supplement describe the sensitivity experiments.

2.4 UFEMISM

We repeated the default experiment with the transient spin-up using the Utrecht Finite Volume Ice-Sheet Model (UFEMISM) (Berends et al., 2021). We use the same transient initialization procedure including the mass change rates from Smith et al. (2020) in the mass balance equation, as well as climate forcing as in the default CISM experiment, but apply a different momentum balance approximation, basal friction law, calving algorithm, and discretization scheme.

Figure 1 Initialized state of the Antarctic Ice Sheet using CISM. (a–b) Ice thickness difference between the final state of the transient initialization and observations for (a) all of Antarctica (RMSE = 34 m) and (b) the ASE region (RMSE = 19.3 m). The black line represents the modeled GL. The observed GL (green) is visible only where the modeled and observed GL deviate in position, for example at the PIG GL. Red and orange contours are drawn around the TG and PIG basins, respectively. (c) Surface ice velocity difference with respect to observations of Rignot et al. (2019) (RMSE = 156 m yr⁻¹). Positive differences indicate locations where CISM overestimates velocities. Observed surface velocities of 100 (outer) and 500 (inner) m yr⁻¹ are contoured by black dots to highlight fast-flowing regions with large velocities and gradients. For the ASE, the velocity RMSE is 99 m yr⁻¹. For ASE shelves only, the velocity RMSE is 216 m yr⁻¹. (d–f) Basal melt rates, with grey and white denoting zero basal melt rates and no data, respectively, of (d) CISM with an integrated flux of 106 Gt yr⁻¹, (e) Adusumilli et al. (2020) with an integrated flux of 94 Gt yr⁻¹, and (f) Rignot et al. (2013) with an integrated flux of 149 Gt yr⁻¹.

UFEMISM is a thermo-mechanical ice-sheet model designed for long-term simulations at variable resolution, e.g., simultaneously high resolution at the grounding line and lower resolution in the slow-moving interior (Berends et al., 2021). In this study, UFEMISM discretizes the hybrid shallow-ice approximation–shallow-shelf approximation (SIA–SSA) to the stress balance (Bueler and Brown, 2009) and other physical equations on a dynamic adaptive triangular mesh. Mesh resolution ranges from ∼4 km at the grounding lines of PIG and TG to ∼60 km for slow-moving inland ice. We apply a Budd-type sliding law (see Eq. 10) which contains spatially varying friction coefficients that aim to represent subglacial substrate conditions:

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{b}}=\tan \left(\mathit{\varphi }\right)u_{\mathrm{b}}^{\frac{\mathrm{1}}{m}}{N}^q.\end{array}$$

We invert for ϕ by nudging toward a target geometry (Bernales et al., 2017a; Pollard and DeConto, 2012). As part of this nudging method, UFEMISM simultaneously inverts for sub-shelf melt rates and the corresponding ocean temperatures needed to reproduce the target geometry over ice shelf sectors (Bernales et al., 2017b). Both calibrations are performed during model initialization, which brings the modeled ice sheet to a state of equilibrium with a prescribed set of boundary conditions. For this study, friction coefficients and ocean temperatures are inverted during the first 90 000 model years of the initialization, after which their values are kept fixed in time. The model is then run for another 10 000 years, which allows it to evolve unperturbed until equilibrium is attained. From 100 000 years, the default experiment with the mass change rates is started. We use both the PMP (default) and FCMP subgrid-scale melting parameterizations from Leguy et al. (2021), as described above. The calving front is conservatively fixed at its present-day position.

3 Results: transient present-day state

Figure 1 shows the modeled CISM state at the end of the transient initialization. The thickness error over the ASE region is low (RMSE = 19.3 m), and the GL closely follows observations, with an average 1.5 km difference (calculated as the average distance between the modeled GL position and the closest observed GL position). The modeled GL for PIG is shifted seaward by 5–10 km. The modeled basal melt rates and melt patterns beneath floating ice agree well with the values from Rignot et al. (2013) and Adusumilli et al. (2020), with an integrated melt flux within the range of the two datasets. The instantaneous and average model drift of this simulation is shown in Figs. S1 and S2, and the inverted quantities of both initializations are shown in Fig. S3. The initialized state in UFEMISM is shown in Fig. S4, the inverted quantities are in Fig. S5, and the result of the CISM equilibrium initialization is in Fig. S6.

Figure 2 Observed and modeled present-day mass change rates. (a) Ice thickness change rates in the observations of Smith et al. (2020) with the observed grounding line in green and the modeled grounding line at the end of the transient initialization in black, in the ASE region. (b–c) Modeled ice thickness change rates when starting a continuation run from the transient initialization after (b) 1 year and (c) 5 years.

The modeled ice surface velocities are in good agreement with observations, even though they are not a target for assimilation. The RMSE of the ice surface velocities is comparable to other ISMIP6 models (Seroussi et al., 2020), of which the range is 100–400 m yr⁻¹. Many of those models are optimized to match observed velocities and use a variety of spin-up and data assimilation methods. Modeled PIG and Thwaites velocities are too low midstream and too high in the lateral margins, possibly because of the lack of brittle processes in CISM. These regions are heavily crevassed (Lhermitte et al., 2020); taking crevasses into account would weaken the ice, speed up the main flowlines, and slow down the margins because the margins are then decoupled from the center flowline. The center flowline can accelerate and is not necessarily “dragging along” the slower moving margins, which would bring our modeled ice surface velocities closer to observations. The integrated ice fluxes across the GLs of TG and PIG are close to the observations of Morlighem et al. (2020) (Fig. S7). As the integrated ice flux is in good agreement with observations, we believe that the dynamical characteristics of the model with a very accurate GL position and ice thickness and ice thickness changes almost equal to observations are sufficient for our continuation experiments.

In the observations of Smith et al. (2020), the ASE shelves thin rapidly. Simulating these large negative mass change rates here requires higher ocean temperatures and lower basal friction under grounded ice compared to the equilibrium initialization. The greater friction during the equilibrium initialization leads to an underestimation of ice velocities (Fig. S6). As a result, the grounding-line flux is lower (see Fig. S7); therefore, lower basal melt rates and in general colder ocean conditions are required to match the observed ice shelf geometry.

The observed and modeled ice thickness change rates in the first years after the start of a default continuation run are shown in Fig. 2. The modeled mass change rate pattern over the first year (b) mimics the observations (a) in detail. After 5 years, the thinning rate of PIG decreases slightly, while the thinning rate of TG increases slightly. The overall pattern of the modeled rates changes little during the first time steps of a continuation simulation, highlighting that CISM is able to reproduce the observed rates robustly.

4 Results: future state

In the default experiment with present-day climate forcing, AIS mass loss is dominated by the collapse of TG (including the neighboring Smith, Pope, and Kohler glaciers) and PIG (Fig. 3). In the first 4 centuries, mass loss is comparable to the current observed rates (Fig. 3b). During this period, the total AIS mass loss averages 0.5 mm GMSL yr⁻¹, slightly less than the current combined loss from PIG and TG, because there is some East Antarctic Ice Sheet thickening, primarily in Dronning Maud Land (Smith et al., 2020), which continues in the forward simulation (Fig. 3a). Areas where little to no mass changes are observed in Smith et al. (2020), such as the Siple Coast and Victoria Land, remain in balance in our simulation.

After 500 years, the mass loss accelerates as both PIG and TG collapse. We can identify a shallow ridge in the TG bedrock profile approximately 45 km upstream of the present-day GL, which acts as the last pinning point before this accelerated collapse begins (see line AB in Fig. 3a and a zoomed-in image in Fig. 4). As soon as the grounding line passes over line AB in Figs. 3a and 4, the high ridge approximately 45 km upstream of the present-day grounding line, the collapse starts. This second ridge at about 100 km upstream in cross section CD has a height similar to ridge AB but is less extended in the cross-flow direction and surrounded by troughs (see Fig. 4).

Figure 3 Future state of the Antarctic Ice Sheet using CISM. (a) Ice thickness change after a 1000-year continuation experiment initialized with present-day ocean thermal forcing. The original modeled GL resulting from the spin-up is shown in green, and the GL position at t=1000 years is shown in black, only where it deviates from the spin-up GL. The thin grey lines show the modeled elevation with an interval of 1000 m at t=1000 years. (b) Ice mass above floatation for Thwaites Glacier (TG) and Pine Island Glacier (PIG) in blue, relative to the initial state shown in Fig. 1, and the corresponding sea level change equivalent (calculated according to Eq. 2 in Goelzer et al. 2020) of TG and PIG (combined, crosses) and the AIS (stars) in red. Present-day observed mass loss is extrapolated in time and labeled as extrapolated observations in black. (c) GL displacement over time and the bedrock profile along cross section CD. The last ridge before collapse is initiated, AB, is shown with a vertical black line. Once the GL passes this ridge, mass loss accelerates. A zoomed-in figure including the 2D bedrock profile is shown in Fig. S8.

Figure 4 ASE bedrock topography and grounding-line migration of the future simulation. The bedrock elevation (m) is shown with grey shading. Magma contours show the grounding-line position at different times during the default continuation experiment. Contour lines clumped together highlight areas where the grounding line stabilizes for long periods. The short white line (AB) shows the last prominent ridge before the accelerated collapse starts, and the long semi-horizontal white line (CD) shows the location of the cross section plotted in Fig. 3.

The TG GL recedes faster than that of PIG (Fig. 5b), in line with the mass loss above floatation (Fig. 3b), where the TG mass loss acceleration precedes the PIG acceleration. In model year 500, ice velocities in the main channel of TG exceed 4000 m yr⁻¹, about twice the current modeled values (Fig. 5e). Thinning is greatest near the new GL, which is on a retrograde bed and therefore deeper than the current GL. Ocean thermal forcing in the ASE increases with depth; therefore newly ungrounded grid cells have a higher basal melt rate and thin rapidly. After 750 years, TG and PIG have collapsed and the basins contain almost no grounded ice, except on a few small ridges with bedrock above sea level. A large, confined ice shelf has formed, with a thickness of several hundred meters near the GLs, decreasing to tens of meters near the (non-moving) calving front. At that time, the two basins have contributed about 1.2 m to global mean sea level. Hence, the retreat of TG and PIG consists of two phases: a gradual decrease for the first few centuries, followed by an accelerated collapse, which continues until the basins are empty. PIG retreat is accelerated by the simultaneous collapse of TG. However, each glacier collapses on its own when the observed mass change rates are applied only to a single basin (not shown).

Figure 5 Development of the collapse at three times. (a–c) Ice thickness difference with respect to the end of the initialization at (a) t=250, (b) t=500, and (c) t=750 years after the start of the default continuation experiment. The present-day GL from observations is shown in green, and the modeled GL is shown in black. The ice surface elevation lines are contours of 500 (dashed), 1000 (dash-dotted), and 2000 m (dotted). (d–f) Ice surface velocities at (d) t=250, (e) t=500, and (f) t=750 years, with the modeled GL shown in cyan and the observed GL in green.

During the continuation experiment, the collapse can be halted by turning off basal melting completely (yellow lines in Fig. 6a). We observe two features in these zero basal melting experiments: the ice sheet readvances more slowly than it collapsed, and the readvance is slower when basal melt rates are switched off later during the simulation. Setting the basal melt rates to 0 is equivalent to decreasing ocean temperatures up to 4 K near the GL and on average 2–3 K elsewhere. We also tested the effect of an instant cool-down of 1, 1.5, and 2 K at 250 years (brown lines) and 1.5, 2, 2.5, and 3 K after 500 years (red lines) during the simulation. None of the cool-down experiments were enough to regrow the ice sheet to the initial condition on the same timescales as the retreat. Most of the experiments merely postpone the total collapse of the TG and PIG basins. Only a full basal melt stop will at any moment halt the collapse. We also tested a percentage decrease in basal melt rates at certain time steps. Setting the basal melt rates to 75 % and 50 % of their present-day value (yellow and brown lines in Fig. 6b) does not stop but only delays the collapse. Setting melt rates of 25 % (red lines) of their present-day value results in eventual regrowth but only when applied in the first 200 years of the simulation (red lines in Fig. 6b).

Figure 6 The impact of decreasing basal melt rates on the ice retreat. (a) Evolution of ice volume (shown as a percentage relative to the initial volume (see Fig. 1) for PIG and TG combined and as sea level equivalent on the right y axis) over 1000 years in experiments with a sudden decrease in basal melt rates. Brown lines show where the thermal forcing (TF) is instantly decreased after 250 years by 1, 1.5, and 2 K, and the red lines show simulations where the thermal forcing is decreased after 500 years by 1.5, 2, 2.5, and 3 K. Yellow lines indicate runs where the basal melting is instantly switched off after 500, 600, 700, and 800 years. (b) Percentage of the basal melt rates applied: 25 % (red), 50 % (brown), and 75 % (yellow). The 50 % and 75 % basal melt rate experiments are started at 50, 100, and 150 years, and the 25 % basal melt rate experiments are started at 0, 50, 100, 150, 200, 400, 500, and 600 years after the start of a transient continuation run.

The accelerated collapse resembles the MISI mechanism (Schoof, 2012, 2007). However, the analytical solution of Schoof (2012, 2007) does not account for an initial buttressing shelf and is valid only when there is no change in buttressing. Recent studies have indicated that the current TG ice shelves provide little or no buttressing (Gudmundsson et al., 2023, 2019; Fürst et al., 2016). During and after the modeled collapse, a large buttressing ice shelf forms. This is at least partly because of our conservative stationary calving front; in reality the front might move inland from its present-day position. As a result, ongoing basal melting is needed to overcome the increased buttressing and sustain the retreat. Removing the basal melt entirely would halt the collapse (see Fig. 6). This makes our retreat at least partly driven by forcing and not strictly an unstable ice dynamical system, as would be the case in the theoretical non-buttressed MISI description.

It is still debated whether the presence of warm CDW, about 3–4 °K above the melting point (Jacobs et al., 2011b), under the PIG and TG shelves is a result of natural variability (Silvano et al., 2022) or forced by climate change (Holland et al., 2019). One study (Holland et al., 2019) found that present-day WAIS basal melting is partly caused by internal variability and partly anthropogenically forced. In the case of natural variability, CDW might be (partly) absent in the future, lowering the melt rates. When we first run 50 years into the future with the default initialization and then instantly halve the melt rates, the collapse is delayed but not halted (Fig. 6b).

Figure 7 Future state of the Antarctic Ice Sheet using UFEMISM. Ice thickness change after a 1000-year continuation experiment initialized with present-day mass change rates. The initial modeled GL resulting from the spin-up is shown in green, and the GL position at t=1000 years is shown in black, only where it deviates from the spin-up GL. For this figure, the UFEMISM data on an irregular triangular grid were interpolated to the same regular rectangular grid as is used in the CISM simulations and smoothed with a Gaussian kernel for visibility.

5 Results: UFEMISM simulation

The default experiment shown in Fig. 3 is repeated with the ice sheet model UFEMISM. The initialized state of the UFEMISM simulation is shown in Fig. S4, in terms of ice thickness and ice surface velocity differences with present-day observations. The inverted quantities are shown in Fig. S5. Figure 7 shows the result of an unforced continuation simulation with UFEMISM, initialized with the present-day mass loss rates and using PMP at the grounding line. The run setup is very similar to the default experiment using CISM but differs, as mentioned in Sect. 2.4, in the stress balance, discretization, sliding-law, and grounding-line parameterizations.

The ice thickness change shown in Fig. 7, after 1000 years of an unforced forward run, is again dominated by the collapse of Thwaites Glacier and Pine Island Glacier, just as in the default CISM simulation shown in Fig. 3. The same floating corridor has formed at the Rutford Ice Stream, connecting the Amundsen Sea and the Filchner–Ronne Ice Shelf. The Antarctic Peninsula is still connected to the Bellingshausen sector via grounded ice, a feature that is absent in the CISM simulation shown in Fig. 2. The UFEMISM simulation predicts slightly less grounding-line retreat. In the remaining areas of the AIS, the ice thickness change is largely negative and close to 0 (on the order of 1–5 m of total ice thickness change over 1000 years), with the exception of some small outlet glaciers in Dronning Maud Land. The mass change patterns on the Ross and Filchner–Ronne ice shelves match the patterns shown in Fig. 3. We therefore conclude that the ASE collapse is likely not restricted to CISM but is rather the result of initializing an ice sheet model with the present-day mass change rates.

6 Model choices and parameter value exploration

The results presented in Sect. 4 are based on the best estimate of the current system state and a single set of preferred model physics choices and parameter values. To assess the sensitivity of our results to those choices and parameter values, we varied them within ranges we viewed as physically realistic. The aim of these experiments is to determine which processes and parameter choices can either accelerate or slow (and perhaps prevent) the projected collapse. Where appropriate, we redid the transient initialization with the new choices (runs with an asterisk (*) in Table 3). In other cases, we recomputed free parameters to start from the same initialized state (same ice thickness and velocities and initial $\mathrm{d}h/\mathrm{d}t$) as our default experiment to minimize the effect of a slightly different initialized state. We verified for the simulations starting with a new initialization that there is little model drift and that PIG and TG would not collapse without applying the mass change rates.

Figure 8 Sensitivity analysis. Ice mass above floatation (left y axis) and equivalent sea level rise (SLR; right axis) in the combined Thwaites Glacier and Pine Island Glacier basins for different experiments (22 in total). The mass change is calculated as the percentage difference relative to the initial value in the TG and PIG basins. The default run shown in Fig. 3 is plotted with black diamonds. The same run without applying the present-day mass change rates, i.e., the model drift experiment, is shown as a dashed line. The UFEMISM runs are shown in green, and the extrapolated observations as in Fig. 2b are shown in marker-free black. Most runs start with the observed mass change rates but deviate early in the run due to dynamic feedbacks. Some runs, e.g., SMB increases by 20 % instantly, start with different mass change rates. This explains the spreading trajectories in the first 100 years. The colored bars at the top indicate the range of times that the simulations need to reach 50 % of their present-day ice volume to highlight the uncertainty related to different categories. See Table 3 for a description of each run.

Table 3 Summary of the sensitivity experiments. Descriptions of the runs shown in Fig. 4. Except as noted otherwise, all settings and resolutions are the same as for the default initialization and continuation run. For some runs, denoted by an asterisk (*), a new initialized state was needed. This new initial state is tested on model drift and verified that a PIG and TG collapse was not found when not applying the present-day mass change rates. An extended description of each run can be found in Table S1. The year of collapse initiation is defined as the year when the rate of basin-wide mass loss exceeds a 10 % change in initial mass per 75 years (0.13 % yr⁻¹), which is equivalent to 1.72 mm GMSL yr⁻¹. The end of the collapse is defined as the first time the mass loss drops below this threshold again.

Figure 8 and Table 3 show the relative mass loss and resulting GMSL rise relative to the initial state for this set of model choices and parameter values. For some simulations, marked with an asterisk (*) in Table 3, a new initialized state was necessary. We made sure that this new initial state was close to the default initial state in terms of thickness, velocity, and grounding-line position RMSE. We ran model drift experiments to verify that, for the new initialized state, the WAIS would not collapse without adding the present-day mass change rates. Notably, the continuation of present-day climate forcing leads to TG and PIG collapse in all experiments, except the default simulation without applying the present-day mass change rates (the model drift experiment, dashed line in Fig. 8). In this experiment, PIG starts to readvance to a downstream seabed ridge, creating some positive model drift.

Once a collapse is initiated (i.e., when the mass loss accelerates), mass loss rates are similar and all experiments show a similar collapse duration (150–250 years; see Table 3). During the fast retreat phase, the mass loss rate is about 3 mm GMSL yr⁻¹ for almost every experiment. This is irrespective of the parameter setting, suggesting that the collapse rate during the fast phase is controlled by the bed topography and/or internal system dynamics rather than detailed physics or parameter choices.

The collapse timing is insensitive to the basal friction law (in line with Barnes and Gudmundsson, 2022), model resolution (8 or 4 km), and a 20 % increase in the integrated SMB. Also, glacial isostatic adjustment (GIA), estimated using an elastic lithosphere–relaxing asthenosphere (ELRA) model (Lambeck and Nakiboglu, 1980) with a relaxation time of 3000 years, has little impact on glacier stability. GIA with a very short relaxation time of 100 years delays the collapse by several centuries but does not prevent it. These findings are similar to those of Berdahl et al. (2023).

The collapse timing is sensitive to model choices related to basal melting. A weaker (stronger) dependence of the basal melt on the thermal forcing at the ice draft depth (Jourdain et al., 2020), i.e., a higher (lower) value for γ₀ in Eq. (5), delays (advances) the start of the collapse (as in Berdahl et al., 2023). Allowing no melt in a grid cell containing the GL (NMP) instead of partial melt (PMP) delays the collapse for several centuries but does not prevent a collapse.

The UFEMISM simulations show a similar pattern of mass loss as the CISM simulations. Again, there is a gradual linear phase of ice mass loss. When approximately 30 % of the original ice mass above floatation in the ASE has disappeared, the accelerated collapse begins; this phase lasts approximately 200 years with FCMP and 250 years with NMP. The similar behavior of the UFEMISM simulations suggests it is unlikely that the WAIS collapse modeled in CISM is a model artifact. The timing at the start of the collapse differs in our sensitivity experiment by about 1200 years. The large difference in collapse timing highlights that the choice of an ice sheet model is a major source of uncertainty.

7 Discussion

We have developed a modified ice sheet initialization method and applied it in Antarctic Ice Sheet simulations using two models, CISM and UFEMISM. This method enables an ice sheet model to be spun up to a state consistent with the observed ice sheet geometry while also matching the observed mass change rates. Previous initialization methods, both spin-up and data assimilation, were unable to accurately represent the observed mass change rates.

In model simulations spanning a wide range of physics and parameter choices, Thwaites Glacier and Pine Island Glacier collapse on a timescale of several centuries without additional climate forcing beyond the ocean forcing required to match present-day melt rates. In other words, the present ocean-driven mass imbalance is a precursor for large-scale deglaciation of the WAIS. A period of gradual retreat over several centuries is followed by a phase of rapid collapse over about 2 centuries, with GMSL rising by about 3 mm yr⁻¹ during the collapse phase. This behavior, including the collapse rate, is consistent with several previous studies (Joughin et al., 2014; Golledge et al., 2021; Coulon et al., 2024; Feldmann and Levermann, 2015). The timing of collapse initiation differs widely across our experiments, as well as in previous studies, and remains a source of large uncertainty (Fig. 8).

Each sensitivity experiment (Table 3) varies a single parameter or physics choice relative to the default experiment. It would be interesting to study simultaneous changes in multiple parameters and physical choices. However, when combining two conservative model choices, e.g., NMP and low melt sensitivity to thermal forcing, the inversion procedure often fails to reproduce the observed GL location, suggesting that the combination is unrealistic. For this reason and because most combinations require a separate, computationally expensive initialization, we did not consider these combined effects. To generate a single conservative end-member, we combined the two most conservative runs (no-melt parameterization at the GL and ELRA isostasy with a relaxation time of 100 years in Table S1) into a “most conservative estimate.” In this simulation, TG and PIG collapse after 2500 years (not shown in Fig. 8). Applying no melt in a cell containing the grounding line is conservative, since the transition between grounded and floating ice occurs over an extended grounding zone possibly influenced by tides (Graham et al., 2022).

We note three missing processes that could affect the ice evolution. First, the mass change rates on grounded ice are controlled by parameters in the basal sliding law, which are not necessarily constant over time. After the initialization phase, C_c is kept constant during forward runs, not accounting for future evolution of the ice–bed interface, such as increased basal meltwater that could enhance sliding (Kazmierczak et al., 2024). In addition, possible feedbacks between the ice flow and the bed (e.g., Rémy and Legresy, 2004; van der Wel et al., 2013) are only partly represented by the basal friction law and the calculation of the effective pressure on bedrock below sea level, representing the hydrological connection to the ocean.

Second, as ice shelves melt and retreat, they add freshwater to the ocean, but the current model setup ignores the potential feedback of the additional freshwater on ocean temperature or salinity. There are several ways for the freshwater flux from a collapsing WAIS to influence the ocean properties of the Amundsen Sea and the Southern Ocean. For example, increased freshwater can cool the Southern Ocean sea surface (Bintanja et al., 2015) and reduce Antarctic Bottom Water formation (Williams et al., 2016). Alternatively, the freshwater input can stratify the ocean just beyond the shelves, thus trapping heat in the subsurface ocean and ultimately increasing heat transport into cavities, resulting in higher basal melt rates (Flexas et al., 2022). It has also been suggested that interactions between freshwater flux and sea ice formation and the presence of warm, dense water in large AIS cavities can stabilize and stratify the ocean flow in cavities, preventing further warm water from flowing into the cavity (Hellmer, 2004). Hence, the exact effect of a large freshwater pulse into the ocean surrounding the AIS remains unknown based on model studies (Swart et al., 2023). Since we use a standalone ice sheet model in this study and no consensus has been reached on the effect of the freshwater flux, we did not parameterize this effect. In future studies, freshwater feedbacks could be added by coupling the ice sheet model to a cavity-resolving ocean model.

Third, bedrock uplifting could slow down glacier retreat in MISI-prone areas (Book et al., 2022; Kachuck et al., 2020; van Calcar et al., 2023). In our default experiment, we do not simulate bedrock height changes. We perform sensitivity tests with CISM's ELRA model of glacial isostasy in Sect. 6. However, these experiments assume, perhaps unrealistically, that the bedrock is in equilibrium at the start of the forward simulation.

8 Conclusions and outlook

We have developed a method to incorporate present-day mass change rates in an Antarctic Ice Sheet model initialization using a spin-up. Given these observed rates, we show that for a wide range of model choices, Thwaites Glacier and Pine Island Glacier will ultimately collapse, likely within 1000 years under present-day conditions. These results are consistent with several previous studies but carried out with an improved initialization scheme, with a much wider range of model configurations, and without additional climate forcing. The timing of collapse onset is highly dependent on model details and varies in our study from 300 to 2500 years. The length of the fast collapse phase (about 200 years) and the maximum mass loss during the collapse (∼3 mm SLE yr⁻¹, sea level equivalent) are consistent across experiments. These experiments suggest that eventual collapse is inevitable for the current climate and can probably be averted only by a large decrease in sub-shelf ocean temperatures. The collapse will likely begin sooner if Amundsen sub-shelf cavities were to warm further in the next few decades, as projected by a recent ocean modeling study (Naughten et al., 2023).

The present-day observed imbalance in the ASE is likely caused by the presence of warm ocean conditions, and the current ice sheet geometry reacts to higher basal melt rates and the loss of buttressing. While some of the warming may be the result of natural variability, it is plausible that the cavities warmed during the past century when a different geometry was present. In that case, GCMs and high-resolution regional models could be used to simulate the coupled ice–ocean–atmosphere processes that might have triggered the recent warming. These models could also explore whether future anthropogenic action could cool these cavities enough to prevent collapse. Replication of our experiments with other ice sheet models could substantiate whether the collapse of Thwaites Glacier and Pine Island Glacier is indeed unavoidable without a considerable reduction in basal melt, as indicated in our study.