Introduction

Human-induced global warming, the temperature increase in combined surface air and sea surface temperatures due to fossil fuel emissions, has reached about 1.0 °C since the preindustrial era, and is currently increasing at 0.2 °C per decade1. As a consequence of continued warming, global mean sea level is projected to rise by 0.25–0.59 m under a low-emissions Representative Concentration Pathways (RCP) 2.6 scenario and by 0.61–1.10 m under a high-emissions RCP8.5 scenario over the course of the twenty-first century2. Just 0.20 m of sea-level rise is predicted to more than double the frequency of extreme coastal flooding events and impair the habitability of low-lying Pacific Islands3. Of all the components of the global sea-level budget, the contribution of the Antarctic Ice Sheet (AIS) is the most uncertain4,5,6. This uncertainty is clear in the recent community effort of the Ice Sheet Model Intercomparison Project 6 (ISMIP6), in which AIS contributions projected over the next century ranged from −0.08 to +0.30 m of sea-level equivalent under an RCP8.5 scenario7.

While some of this range is related to uncertainties in the climate response to a given emissions scenario and differences in model initialization8, ice sheet model projections are also subject to uncertainties in the parameterization of physical processes9,10,11. In particular, many ice sheet projection studies have emphasized uncertainties related to the parameterization of sub-ice shelf basal melting and surface mass balance7,8,12,13,14. However, unconstrained model parameters related to the underlying bedrock and substrate conditions, basal friction, and ice rheology, which are difficult to observe, influence ice sheet sensitivity to climate forcing on various timescales15,16,17,18. As a result of all of these uncertainties, AIS does not show clear emissions scenario dependence in twenty-first century sea-level projections6, and the timescale at which emissions scenario dependence will emerge is unclear. To address this, we first apply a combination of ice sheet model simulations and a statistical emulator to quantify the effect of ice flow-related uncertainty and determine the time of emergence of AIS sea-level contribution under different RCP scenarios. We then relate our findings to global warming levels (GWLs), defined as 20-year moving windows with respect to the preindustrial 1850–1900 value, to provide AIS sea-level projections per degree of warming to 2300.

Results

Historically constrained sea-level projections

We perform ensembles of ice sheet simulations using the Parallel Ice Sheet Model (PISM) for RCP2.6 and RCP8.5 scenarios following ISMIP6 protocol19, but extending the simulations to 2300 (see Methods). Our two standard RCP2.6 and RCP8.5 ensembles use climate forcing (i.e., surface air temperature, precipitation, ocean temperature, and salinity) derived from the NorESM1-M of the Coupled Model Intermodel Comparison Project 5 (CMIP5)20. This atmosphere-ocean general circulation model (AOGCM) was selected as a Tier 1 model of ISMIP6 due to its low bias with respect to historical metrics of Antarctic climate and Southern Ocean conditions21. In the simulations, we investigate a full-factorial parameter sampling of 4 unconstrained model parameters related to ice flow, rheology, and basal shear stress (Supplementary Table S1). The parameters include enhancement factors of the shallow ice and shallow shelf approximations (ESIA and ESSA), a sliding law exponent (q) and the minimum till friction angle (ϕmin). The enhancement factors control the creep and sliding velocities, whereas the latter parameters influence the basal shear stress.

The process-based ice sheet model simulations using the standard RCP2.6 and RCP8.5 forcings are then used to train a statistical emulator that employs Gaussian Process (GP) regression (GPR) (Fig. 1). While previous studies have employed similar statistical techniques for determining the probabilistic sea-level contribution of the ice sheet or ice sheet catchment for a given time period10,12,22, here we also consider the temporal evolution of the sea-level contribution. In particular, the GPR uses the four ice sheet model parameters, and a combination of (i) the direct effect, (ii) the cumulative effect, and (iii) the committed effect of global warming as independent variables. This is realized in terms of (i) the global mean temperature (from the RCP scenarios), (ii) the cumulative global mean temperature, and (iii) the time since the last global mean temperature change. The latter effect is important to capture the committed sea-level rise that continues even after global mean temperatures are held constant after year 2100. Each point in time creates a training data point, for a total of 45,846 data points from which we use 5% for training. Once fitted, the emulator is in excellent agreement with the individual ice sheet model simulations (Supplementary Fig. S7 and Supplementary Table S2) and able to explore the sensitivity of Antarctica’s sea-level contribution to model parameters much more efficiently and, therefore, in much greater detail than could be achieved with process-based ice sheet modeling alone. Furthermore, we can explore a whole range of different scenarios, for example, additional RCPs, such as RCP4.5 or RCP6.0 (Fig. 1).

Fig. 1: Historically constrained AIS sea-level contribution projections.
figure 1

a AIS sea-level contribution (m) projected by the statistical emulator relative to 2000. Gray shading indicates the full ice sheet model range, with darker color marking overlap between the two RCP scenarios. Color shading shows the historically constrained projections for the RCP scenarios after discounting parameter ranges that fall outside historical AIS sea-level contribution estimates25. The dark blue and red lines show the median sea-level contribution for the RCP2.6 and RCP8.5 scenario, respectively. The two dotted lines indicate the likely and the very likely time of emergence, defined as the time when the difference between the RCP8.5 and RCP2.6 medians is larger than the 68% and 95% confidence interval of the RCP2.6 ensemble. The right panel shows projected AIS contributions at 2300 relative to 2000 of other studies compared to this study (LOW21). GOL15 refers to Golledge et al. (2015)56, BUI19 refers to Bulthuis et al. (2019)73, and BAM19 refers to Bamber et al. (2019)74. EDW19 refers to the emulated projections of Edwards et al. (2019)10 using the Deconto and Pollard (2016)45 ice sheet simulations without marine ice cliff instability (MICI). Inclusion of MICI in the ice sheet simulations increases the upper bounds of the RCP4.5 and RCP8.5 scenarios by >3.5 m10. b Probability distribution functions (PDFs) of AIS sea-level contribution (m) for the RCP scenarios through time. The vertical axis is the density of the probability. The PDFs are for a normal distribution fitted to the emulator runs, and agree well with their raw histograms.

Hindcasting can be a useful tool for ice sheet model validation by ensuring a model can accurately capture observed rates of change23. Here, we use the emulator to hindcast ice mass changes for the full parameter space through the past decades and compare to historical estimates of the AIS mass loss24,25,26 (Supplementary Fig. S8). We discount parameter combinations that produce sea-level contribution trajectories that exceed or underestimate these estimates, an assumption that implies the current ice sheet is not in steady state. While the full range of projected sea-level contributions at 2100 exceeds those of ISMIP6 for the two RCP scenarios7, considering only the historically constrained parameter combinations eliminates the most extreme scenarios. For these eliminated parameter combinations, the ice sheet is either especially sensitive or resistant to climate forcing, with reduced basal shear stress and high enhancement (i.e., low ϕmin, high q, ESIA and ESSA) associated with higher sea-level contributions and vice versa. The historically constrained parameter space narrows the projected ranges to 0.12–0.44 m for RCP2.6 and 0.21–0.56 m for RCP8.5 at 2100 relative to 2000 (95% confidence interval (CI), blue and red shading in Fig. 1, respectively). These historically constrained ranges are higher than the ISMIP6 projections for both scenarios, and consistent with observations that ice mass loss of AIS and the Greenland Ice Sheet has generally exceeded ice sheet model predictions in recent years27,28. The longer-term response to 2300 results in ranges of 0.45–1.57 m and 1.96–3.79 m for RCP2.6 and RCP8.5, respectively.

The median sea-level contributions of the two RCP scenarios diverge after 2050, but the model parameter-related uncertainty remains high even with the parameter constraints, with substantial overlap between the ranges of the trajectories of the RCP2.6 and RCP8.5 scenarios for the following century (Fig. 1). The time of emergence of RCP8.5 relative to RCP2.6, defined here as the point in time when the difference between the medians of the scenarios is greater than the 95% CI of the RCP2.6 scenario, occurs at 2189. At this point in time, the distributions of the sea-level contribution are distinct between the emissions scenarios, with the difference increasing to the end of the simulation period.

Spatial patterns of ice sheet change

The ice sheet simulations that use parameter combinations within the range acceptable for historical sea-level contributions can provide insight into the controlling physical processes (Supplementary Table S2 and Supplementary Figs. S8S10). The RCP2.6 and RCP8.5 ensemble averages show similar spatial coverage of ice thinning and grounding-line retreat of the Amundsen Sea Embayment (ASE) and East Antarctic outlet glaciers at 2100 relative to 2015 (Fig. 2a, b). The clearest contrasts between the scenarios occur with the larger Ronne–Filchner (Weddell sector) and Ross ice shelves, with grounding-line advance of the Siple Coast of the Ross ice shelf in RCP2.6 versus thinning and collapse of the shelf in RCP8.5. The respective ensembles display low standard deviation (SD) in ice shelf thickness because floating ice is unaffected by model parameters that influence the basal shear stress (Fig. 2e, f). However, both scenarios indicate high uncertainty (SD > 100 m) of the ASE and Wilkes Land outlet glaciers.

Fig. 2: Modeled ice thickness changes and threshold exceedance.
figure 2

ad Ice thickness change (m) for the period of 2015–2100 for the RCP2.6, RCP8.5, RCP8.5-CO, and RCP.85-WO ensemble averages. eh Ice thickness standard deviation at 2100. il Time when the difference between the scenario-based ensemble average and a control (constant climate) ensemble average exceeds the standard deviation of the control ensemble. ASE and WL refer to the Amundsen Sea Embayment and Wilkes Land, respectively.

Accounting for the high variance in regional Southern Ocean temperatures among RCP8.5 climate projections of different AOGCMs of CMIP521, we also consider warm and cool ocean RCP8.5 cases from other models (RCP8.5-WO and RCP8.5-CO, respectively; see Methods). In comparison to the standard RCP8.5 case, the ensemble averages for the RCP8.5-CO and RCP8.5-WO scenarios display differing amounts of ice shelf thinning in accordance with differences in sub-ice shelf melt rates. The Ross ice shelf thins but remains intact in the RCP8.5-CO ensemble average, with little change in grounding-line position (Fig. 2c). The RCP8.5-WO ensemble average shows substantially more thinning in the Ross sector and partial collapse of the Ronne–Filchner ice shelf (Fig. 2d). Both ensemble averages display similar ice sheet thinning in ASE and Wilkes Land in the ensemble average, but with high SDs in ice thickness for the full ensemble (Fig. 2g, h).

In terms of emergence between scenarios, we consider the time when the ice thickness difference for a given climate scenario to a control ensemble (constant climate) exceeds the SD of the control ensemble, referred to here as threshold exceedance. The change in ice shelf thickness exceeds this threshold in the first decade of the simulation in the Ross Sea sector for each scenario, and in the Weddell sector and Antarctic Peninsula with the RCP8.5-CO and RCP8.5-WO forcings (Fig. 2i–l). Threshold exceedance occurs across all of the remaining ice shelves within the following decades. In terms of the grounded ice, the majority of the ice sheet does not exceed the threshold with the RCP2.6 climate forcing (Fig. 2i), with the exception of parts of the Ross Sea sector (before 2025). In contrast, threshold exceedance progresses into the interior of WAIS from the Ross, Amundsen, and Weddell Sea sectors under the RCP8.5 climate forcings (Fig. 2j–l). The timescale is on the order of decades in RCP8.5-WO, but generally >100 years in RCP8.5 and RCP8.5-CO for the majority of WAIS.

Regions sensitive to ice flow parameters

The region of highest sensitivity to the ice flow parameters is the ASE (Fig. 2e–h), which is currently experiencing the greatest mass loss of the entire ice sheet due to high rates of sub-ice shelf melting25,29. Importantly, the Thwaites and Pine Island Glaciers in the ASE overlie a deep marine basin with downward-sloping basal topography that lies substantially below the sea level (Fig. 3a), which creates an ice sheet configuration prone for marine ice sheet instability (MISI)30. While ice sheet models suggest that continued ocean warming will lead to ice shelf thinning and grounding-line retreat in this sector13,31, basal friction and ice rheology influence the rate of retreat16,17, which can have a large impact on the sector’s sea-level contribution32.

Fig. 3: Surface velocity changes in the ASE.
figure 3

a Initial basal topography and b surface ice velocity of the Amundsen Sea Embayment (ASE). PIG and TG refer to Pine Island Glacier and Thwaites Glacier, respectively. ch Change in ice surface velocity at 2100 relative to the start of the projection experiment (2015) for the control (CTRL), RCP2.6 and RCP8.5 simulations with the c, e, g minimum and d, f, h maximum AIS sea-level contribution, respectively. The control experiments use present-day climate for the duration of the model run. i, j Rate of sea-level contribution for the simulations shown in panels c through h.

In our simulations, the Thwaites and Pine Island Glaciers show increased ice surface velocity relative to the initial condition, regardless of model parameter combination or climate forcing. However, increasing the enhancement factors and decreasing q, which allows the ice to flow and slide more easily, leads to a larger area of increased ice surface velocity regardless of the RCP scenario (i.e., comparing Fig. 3c, e, g to Fig. 3d, f, h). In terms of the rate of sea-level contribution, the simulations show significant overlap between the RCP2.6 and RCP8.5 scenarios in the ASE sector until 2070 (Fig. 3j), when RCP8.5 ocean temperatures increase substantially (Supplementary Fig. S13). This results in a faster flowing ice shelf by 2100 and a more rapid increase in the ASE sea-level contribution in the following decades in the RCP8.5 scenario.

Similar to the ASE, the rates of ice loss under our two emissions scenarios for a given parameter combination start to diverge around 2070 for the entire ice sheet (Fig. 3j), consistent with previous studies10,33, though scenario overlap persists until approximately 2100. The increase in the rate for the RCP8.5 scenario is more significant for the model parameter combination with the highest historically constrained sea-level contribution, as the ice is more sensitive to ocean warming in this case. In contrast, the RCP2.6 rate of sea-level contribution remains relatively constant for AIS, but depends on the parameter combination.

The longer-term 285-year response of AIS under a modern climate forcing shows continued mass loss in both the ASE and WL for the full range of historically constrained parameters (Fig. 4). The main difference in ice thickness change between the simulations with the lowest and highest sea-level contributions is the magnitude of ice thickness change in WL and the wider EAIS. On top of this committed response, the RCP forcings result in additional decreases in ice thickness in the ASE and WL, though the decrease is much more substantial for RCP8.5. In contrast to these parameter-dependent regions, the Ross Ice Shelf and Ronne–Filchner Ice Shelf are more strongly influenced by the climate forcing. As a result, the key difference between the long-term forced responses between the RCP scenarios is the disintegration of WAIS in RCP8.5 from the collapse of the ice shelves and more significant retreat in the ASE.

Fig. 4: Committed and forced ice sheet responses.
figure 4

ac Committed ice thickness changes (m) at 2300 relative to 2015 for the parameter combinations with the minimum, ensemble average, and maximum sea-level contributions under present-day (PD) climate forcing (Ctrl). df Additional ice thickness change (m) from the RCP2.6 forcing (i.e., RCP2.6 - Ctrl). gi Additional ice thickness change (m) from the RCP8.5 forcing (i.e., RCP8.5, Ctrl). Black texts in the top left corners indicate the total sea-level contribution relative to 2015. Navy texts indicate the sea-level contribution of the forced response only (i.e., RCP experiment, Ctrl). ASE and WL refer to the Amundsen Sea Embayment and Wilkes Land, respectively.

Influence of basal and surface mass balance changes

The primary process driving the divergence in sea-level contribution between RCP2.6 and RCP8.5 is ocean thermal forcing, which controls ice shelf basal melt rates. Surface air temperature and precipitation forcing play a relatively minor role in terms of the difference in sea-level contribution, with relatively higher precipitation having a slight compensating effect in RCP8.5 (Supplementary Fig. S12). However, other AOGCMs used as forcing in ISMIP6 display a greater increase in snow accumulation in high-emissions scenarios19; as a result, some ISMIP6 models project negative end-of-century sea-level contributions with forcing from a high-emissions CMIP6 Shared Socioeconomic Pathways (SSP) 585 scenario (roughly equivalent to CMIP5 RCP8.5), despite the relatively higher increase in ocean temperatures and basal melt rates14,34.

In comparing simulations using the CNRM-CM6-1 SSP126 (roughly equivalent to RCP2.6) and SSP585 forcings within the historically constrained ice flow parameter range, we observe little difference in ice thickness changes in grounded ice regions or sea-level contribution at 2100 between the scenarios (Fig. 5). However, in extending the simulations to 2300, the thinning and collapse of the Ross and Ronne–Filchner Ice Shelves and retreat in the ASE lead to significant ice loss in WAIS. The compensating effect of surface mass balance to ocean-driven ice mass loss leads to an even longer period of overlap between the SSP126 and SSP585 sea-level trajectories than the RCP2.6 and RCP8.5 cases (Fig. 5e).

Fig. 5: Projected ice sheet changes for SSP scenarios.
figure 5

Ice thickness change (m) at a, b 2100 and c, d 2300 relative to 2015 for the historically constrained SSP126 and SSP585 ensemble members with the median sea-level contribution, respectively. e Sea-level contribution timeseries for the historically constrained SSP126 and SSP585 ensemble members with the minimum, median, and maximum sea-level contribution. Light red and blue shading fills the area between the SSP585 minimum and maximum and SSP126 minimum and maximum, respectively. Purple shading indicates SSP scenario overlap.

The main driver of scenario dependence in these simulations is therefore the breakup of the large ice shelves. In the decades prior to the time of emergence of RCP2.6 and RCP8.5 (2189), the difference in sea-level contribution between the ensemble averages increases most in the Ross and Weddell Sea Sectors (Fig. 6a). In particular, the Ross Ice Shelf shows refreezing and ice sheet advance in RCP2.6, but high melt rates and ice shelf disintegration by 2100 in RCP8.5 (Fig. 6b–g). The refreezing, which is also observed with SSP126 forcing, may be an artifact of the basal melt parameterization resulting from cooling and freshening of the Ross Sea in the climate models. Modern observations indicate that refreezing of this magnitude occurs as the result of supercooling of buoyant meltwater plumes29. The scenario difference between the modeled melt rates of the Ronne–Filchner Ice Shelf is less extreme than in the Ross, but the RCP8.5 does show partial ice shelf collapse by 2150. Although both scenarios show relatively high basal melt rates in the ASE, the difference in sea-level contribution increases over this time period.

Fig. 6: Scenario differences by AIS sector and sensitivity to basal melt rates.
figure 6

a Difference in ensemble average sea-level contribution (RCP8.5–RCP2.6) by AIS sector. We use the same sector boundaries as defined in Levermann et al. (2020)75. Decadally averaged basal melt rate in 50-year increments for bd the RCP2.6 ensemble average and eg the RCP8.5 ensemble average. h Sea-level contribution (m) relative to 2020 versus cumulative basal melt (Gt yr−1) in the historically constrained ice sheet simulations for the same time period as panel a, with color indicating the climate forcing scenario.

Considering the parameterization of ice shelf basal melting is highly uncertain, but consequential for ice sheet projections7, we perform ice sheet simulations using a wide range of climate forcings, which together cover a large range of potential basal melt rates. For the same climate forcing, the ranges in both the sea-level contribution and cumulative melt rate are due to the uncertainty of ice flow parameterization (Fig. 6h). Particularly at lower levels of cumulative basal ice shelf melt, the range in sea-level contribution is greater than 1 m. In general, the sea-level equivalent sensitivity to cumulative basal melt increases, with some deviation from the SSP585 scenario with relatively higher snow accumulation compared to the other climate forcings.

Antarctic changes with global warming levels

Given that the emulator is trained with end-members of the RCP scenarios (i.e., RCP2.6 and RCP8.5), it can be applied to the intermediate scenarios of RCP4.5 and RCP6.0 without running additional process-based ice sheet model simulations. These intermediate scenarios achieve GWLs of 2.3 and 2.6 °C relative to the Preindustrial Era, respectively, as compared to 4.0 °C of warming in the high emission RCP8.5 scenario (Fig. 7a). Despite these substantial differences, the range of the projected AIS sea-level contribution is nearly indistinguishable between these scenarios at 2100 (Fig. 7b). By 2300, although the AIS sea-level contribution has a multi-meter gap between RCP2.6 and RCP8.5, the RCP4.5 and RCP6.0 scenarios overlap with all of the emissions scenarios, indicating that a much longer timescale is required to achieve a scenario-based divergence. However, it should be noted that this long timescale is partly due to the constant climate conditions of the ice sheet simulations after 2100. If mean global temperature continued to increase after 2100, the divergence between scenarios would occur earlier and the sea-level commitment would be higher; hence, the emulator is projecting the minimum sea-level commitment for these RCP scenarios.

Fig. 7: Projected AIS sea-level contribution per GWL.
figure 7

a Smoothed global mean surface air temperature of the RCP scenarios, which is the forcing for the statistical emulator (details about the smoothing can be found in the Methods). b Left: sea-level rise per °C global warming relative to 2000 for each RCP scenario up to 2100. Right: range of emulated AIS sea-level rise contribution at 2100, 2200, and 2300 relative to 2000 for each RCP scenario. The different shadings indicate the 50, 90, and 95% percentile ranges of the historically constrained emulator runs. The circles indicate the median contribution of a given scenario and the gray bar shows the difference in y-axis scales. c, d Same as a and b, but for idealized experiments of 1, 2, 3, 4, and 5 °C global warming by 2100. e Idealized experiments altering the rate at which 2 °C warming is reached. f Left: sea-level rise per °C global warming relative to 2000 for each warming rate scenario up to the time at which 2 °C is reached. Right: range of emulated AIS sea-level rise contribution at 2100, 2200, and 2300 relative to 2000 for each warming rate scenario.

In considering that the Paris Agreement sets a target of well below 2 °C of global warming, we explore how 2 °C of warming by the end of the century compares to different warming magnitudes (Fig. 7c, d). In these GWL experiments, the global mean temperature is linearly increased to a given GWL and then held constant. Similar to the RCP experiments, the projected sea-level contribution at 2100 shows overlap among a GWL range between 1 and 5 °C (Fig. 7c, d). The projected range of sea-level contribution under a 2 °C GWL falls entirely below the 4 °C GWL range at 2300, but still shows overlap with the 3 °C GWL range. If we vary the rate of warming to the 2 °C global warming target over the next century (Fig. 7e, f), the sea-level contribution at the time of 2 °C warming is lower if this occurs earlier. However, the long-term committed response is virtually identical; in comparing 2 °C warming by 2020 with 2 °C warming by 2100, the difference in median sea-level contribution in the year 2300 is only 2 cm. This result implies that as long as global warming stays below 2 °C within this century, the committed AIS response will remain similar, independent of the pathway to 2 °C of warming. However, exceeding the 2 °C global warming target by just 1 °C increases Antarctica’s median sea-level contribution by another 69 cm, a 50% increase, at 2300.

Discussion

A number of climate parameters, including surface air temperature35,36, precipitation37, sea surface temperature, and ocean carbon uptake38, are expected to permanently exceed recent natural variability within decades due to anthropogenic influence. The projections analyzed here, however, suggest that a high-emissions signature of the AIS sea-level contribution will not unambiguously emerge from the wide potential range of low-emission sea-level projections for over 100 years due to current limitations in our understanding in ice flow and sliding and despite the significant climatic differences between scenarios. Our process-based ice sheet model simulations suggest that widespread thinning and partial collapse of the Ross and Ronne–Filchner ice shelves arising from decades of accumulated oceanic warming will be the earliest warning sign of the multi-meter contributions of AIS that occur under high-emissions scenarios. Observations indicate that such warming trends have already started in the vast majority of the global upper ocean to the depths of modern ice sheet grounding lines over the last decades39,40, though observations are especially limited in the Southern Ocean41.

In comparison to other ice sheet models, our results indicate a markedly higher future AIS sea-level contribution than other ISMIP6 models by 21007. Two main factors contribute to this. First, ISMIP6 protocol required subtracting the committed ice sheet response from the total sea-level contribution19, in order to accommodate models that rely heavily on data assimilation in the initialization, and are therefore subject to model drift due to inconsistencies in input datasets42,43. Since we use a fully thermally evolved model tuned to the dynamic equilibrium of present day, we include the committed AIS response, which is non-trivial but highly uncertain, especially over long timescales (Fig. 4). Secondly, and more importantly, the historical constraints limit the plausible range of our projections. Notably, ISMIP6 models have also been shown to underestimate ice loss of the Greenland Ice Sheet as compared to observations28.

In addition to ISMIP6, our results are consistent with other recent ice sheet projection studies that show divergence in the rate of ice loss between different emissions scenarios between 2060 and 2080 due to enhanced thinning and retreat of Thwaites Glacier in the ASE10,33,44. This time period marks the onset of accelerated ice sheet retreat under the RCP8.5 scenario that in our experiments leads to partial collapse of WAIS by 2300 (Fig. 4h). In other models that include a parameterization of continual ice cliff failure triggered by the removal of ice shelves, the projected sea-level contribution is considerably higher at 2300 for RCP8.544,45; however, dynamic thinning and the resistive force of calved debris may reduce the potential of such catastrophic collapse from occurring in vulnerable locations, such as Thwaites Glacier46.

As ocean thermal forcing is the primary driver of AIS mass loss, the main uncertainty affecting the prediction of the AIS sea-level contribution in addition to ice sheet flow and sliding is therefore the rate and magnitude of Southern Ocean heat uptake and redistribution. In particular, the CMIP models used as forcing for the ice sheet model simulations differ in the overall magnitude of Southern Ocean warming they predict, particularly in the Ross and Weddell Seas21. It has been estimated that the Southern Ocean accounts for about 75% of global ocean heat uptake47, and the combined effects of wind and meltwater may be causing local warming to accelerate by increasing upwelling of warm, saline circumpolar deepwater41. Meltwater feedback, not considered in the CMIP5 and CMIP6 models used here48, has also been suggested to act as a positive feedback that increases AIS ice loss33. Since no mechanism currently exists for removing heat from the ocean, future observations of substantial ice shelf thinning and collapse may imply that a higher-end AIS sea-level contribution is already committed.

Here, we have addressed the significant level of uncertainty in AIS projections under given emissions scenarios and GWLs. The results demonstrate that the total warming that occurs over the twenty-first century controls the resulting AIS sea-level contribution, independent of the rate of warming. However, despite recent improvements in mapping bed topography49, the limited constraints on substrate conditions and ice rheology, which impact ice flow and sliding, are especially important over decadal timescales in regions with retrograde bed slopes prone to MISI, such as the ASE and WL16,49,50,51. Optimizing ice flow simulation in these regions with additional observations may therefore considerably reduce the range of the potential AIS sea-level contribution projected under a given RCP scenario or GWL. In contrast to MISI-prone regions, ice thinning of the Ross and Ronne–Filchner ice shelves is highly climate scenario-dependent on a scale of years to decades, driven by their sensitivity to ice shelf basal melt rates. However, the effect on the sea-level contribution from these sectors is not observed without at least partial ice shelf collapse, which occurs on centennial timescales under GWLs > 2 °C. The disintegration of these largest AIS ice shelves under these high GWLs ultimately results in widespread retreat of WAIS, which is likely not reversible on human timescales52, and a sea-level commitment that is multiple meters higher than the potential range under a GWL limited to 2 °C.

Methods

Ice sheet model setup

We perform ice sheet simulations using the PISM version 1.1, a three-dimensional thermo-mechanical ice sheet-ice shelf model53. The stress balance model of PISM is a hybrid approximation of the Stokes stress balance, in which velocities are calculated by the superposition of the shallow ice approximation (SIA), which dominates in grounded ice regions, and the shallow shelf approximation (SSA), which dominates in ice shelves and ice streams and determines the basal sliding velocity of grounded ice54. The ice sheet grounding line, the transition between grounded and floating ice, is refined by a sub-grid-scale parameterization that smooths the basal shear stress field55, and can migrate freely. We do not include a sub-grid basal melt scheme for partially floating cells56, however, as this may overestimate the rate of grounding-line retreat57.

To initialize the ice sheet model, we use an identical approach to the Victoria University of Wellington (VUW) contribution to ISMIP67. Starting from initial bedrock and ice thickness conditions from Morlighem et al. (2020)49 and reference climatology from Van Wessem et al. (2014)58, we perform a multi-stage spin-up that guarantees well-evolved thermal and dynamic conditions without loss of accuracy in terms of geometry. This is achieved through an iterative nudging procedure described in Golledge et al. (2019)33, in which incremental grid refinement steps are employed that also include resetting of ice thicknesses to initial values. Drift is thereby eliminated, but thermal evolution is preserved by remapping of temperature fields at each stage. We start with an initial 32-km resolution 20-year smoothing run in which only the SIA is used. Then, holding the ice geometry fixed, we run a 250,000 year, 32-km resolution, thermal evolution simulation in which temperatures are allowed to equilibrate. Refining the grid to 16 km and resetting bed elevations and ice thicknesses, we run a further 1500 years using full model physics and a present-day climate. After resetting bed elevations and ice thicknesses, we perform a 65-year historical simulation (1950–2015). This is a longer historical run than the VUW-PISM contribution to ISMIP6, in order to account for the historical effects of ocean forcing in projections using PISM59.

Ice sheet model projections are sensitive to the initialization method and resulting state8,60. The initialized ice sheet configuration resulting from our spin-up is similar to the VUW contribution to ISMIP67, and well matches observations in terms of ice geometry, grounding-line position, and ice dynamics (Supplementary Fig. S1). In comparison to other ISMIP6 models, the root mean square error (RMSE) is within the range for both ice thickness and ice surface velocity (Supplementary Fig. S2). To account for uncertainties in terms of the model parameters used in this spin-up, we performed sensitivity experiments using parameter combinations with relatively higher and lower sensitivity to climate forcing. These alternative initialized states also show RMSE within the range of the ISMIP6 models. Projections from these initial states indicate a consistent time of emergence within the historically constrained parameter range (Supplementary Fig. S3).

We employ an ocean model based on Holland and Jenkins (1999)61 to compute basal ice shelf melt rates and temperature from thermodynamics in an ocean boundary layer. Our approach to optimizing basal melt rates is similar to that of Golledge et al. (2019)33 in that we apply empirical tuning using a melt factor that reproduces present-day melt rates (Supplementary Fig. S4). The targets of our model are the melt rate observations of Depoorter et al. (2013)62, which indicate a total basal mass loss of 1454 ± 174 Gt yr−1, with annual melt rates of up to 16 m yr−1 of Pine Island Glacier and especially low melt rates in the cold ocean cavities of the Ross and Ronne–Filchner ice shelves. Our model is initialized with 1443 Gt yr−1, with a maximum Pine Island Glacier melt rate of 17 m yr−1. A comparison of the model results presented here to other ISMIP6 models for a control run with present-day climate is shown in Supplementary Fig. S4. In the projection experiments, the ocean temperature and salinity fields are applied as anomalies relative to the present-day field, though the optimization using the dimensionless coefficient as a melt rate factor introduces uncertainty, as it assumes that the relationship between temperature, salinity, and the melt rate factor is unchanged through time. We address this by performing experiments with multiple CMIP5 models and one CMIP6 model, which together produce a wide range of future melt rates, as explained in the subsequent section.

Ice sheet model projections and parameters

The projection experiments are initiated in the year 2015 from the same dynamically equilibrated state following the ice sheet model spin-up and historical simulation. Simulations are run for 285 years at 16-km horizontal resolution. While this is a relatively coarse resolution, this enables us to perform a greater number of simulations for a wider range of parameter combinations and climate forcings, though uncertainty related to model resolution exists. For the main RCP2.6 and RCP8.5 model ensembles, we use climate forcing derived from the NorESM1-M CMIP5 model20, the only Tier 1 model of ISMIP6 that included both RCP scenarios7,21.

The climate forcing includes surface mass balance and temperature anomalies over the ice sheet and ocean temperature and salinity fields extrapolated to sub-ice shelf cavities63. The first 15 years of climate forcing are identical, at which point the RCP2.6 and RCP8.5 climate forcings diverge. At the year 2101, the climate anomalies remain constant for the respective RCP scenario for remainder of the ice sheet simulation to the year 2300. Ideally, these simulations would be performed with extended CMIP runs to 2300, though for consistency we use models selected for ISMIP6. Given that temperatures do not continue warming after 2100, these projections represent a minimum contribution in the subsequent centuries. Also, considering the substantial differences in regional ocean temperature changes in the Southern Ocean of CMIP5 models21, we also perform separate parameter-constrained ensembles using climate forcing from the other ISMIP6 Tier 1 and 2 models, including CCSM4 RCP8.5 (RCP8.5-CO), HadGEM2-ES RCP8.5 (RCP8.5-WO), CNRM-CM6-1 (SSP126 and SSP585), and IPSL-CM5A-MR (RCP2.6 and RCP8.5), as well as with constant modern climate conditions (control).

Model ensemble members are performed using different combinations of model input parameters related to ice flow and sliding in order to assess their individual and combined impacts. We consider four model parameters: SIA and SSA enhancement factors, a sliding law exponent (q), and the minimum till friction angle (ϕmin). The enhancement factors control the creep and sliding velocities, whereas the basal substrate terms q and ϕmin are related to the basal shear stress64. These specific parameters were selected given their known influence in ice sheet model simulations, yet poor constraints11, and are described in more detail below. We use three values per four model parameters for a total of 81 combinations (see Supplementary Table S1). The parameter values cover a wide range of values previously used in PISM simulations of the AIS11,18,33,53,65. Comparison of the simulations to observed Antarctic mass loss of the last decade25 narrows the range of parameter space, allowing for a parameter-constrained ensemble. Lastly, to test and enhance the statistical emulator (described in detail below), we perform additional experiments with random parameter combinations within ranges of high uncertainty using the main RCP2.6 and 8.5 climate forcings.

Enhancement factors for each part of the stress regime (i.e., SIA and SSA) are used in PISM, as in most continental-scale ice sheet models, to account for the anisotropy of ice and are applied uniformly53. This allows the creep and sliding velocities to be optimized using simple coefficients applied to the respective flow laws. For grounded ice, the horizontal shear stress is the dominant component of the stress regime, whereas the ice shelf is dominated by horizontal tension; as such, realistic values for enhancement factors differ for the SIA and SSA, with ESIA > 1 and ESSA < 166,67. Increasing ESIA allows the ice to deform more easily, and increasing ESSA produces faster ice streams and thinner ice shelves. Based on enhancement factors previously used for AIS simulations with PISM11,18,33,53,65, we use values of 1.2, 2.4, and 4.8 for ESIA and values of 0.4, 0.6, and 0.8 for ESSA.

PISM parameterizes ice sheet sliding in the form of a power law that ranges from plastic Coulomb sliding to a linear sliding law54,64:

$${\tau }_{b}=-{\tau }_{c}\left(\frac{u}{{u}_{{{{{{\rm{threshold}}}}}}}^{q}| u{| }^{1-q}}\right)$$
(1)

In the power law, τb is the basal shear stress, τc is the yield stress, and u is velocity, where uthreshold is a threshold velocity of 100 m yr−1. The q term is the sliding exponent parameter, where q = 0 is purely plastic sliding, and q = 1 is sliding linearly related to the applied stress. The default value for this parameter in PISM is q = 0.25. We use values of 0.25, 0.50, and 0.75 in the ice sheet model ensemble.

The material properties of the till are related to ϕ:

$${\tau }_{c}={c}_{0}+(\tan \phi ){N}_{{{{{{\rm{till}}}}}}}$$
(2)

where c0 is the till cohesion, Ntill is the effective pressure of the till, and ϕ is the till friction angle. As in previous studies11,23,53,54, ϕ is heuristically determined as a piecewise linear function of the bed elevation, based on the assumption that lower-lying till with a marine history is weaker68:

$$\phi(x,y) =\left\{\begin{array}{ccc}\phi_{\min},\hfill & {{{\mbox{if}}}}\; b(x,y) \leq b_{\min}\hfill \\ \phi_{\min} + (b(x,y)-b_{\min})M, & {{{\mbox{if}}}}\; b_{\min}\, < \, b(x,y)\, < \, b_{\max} \\ \phi_{\max},\hfill & {{{\mbox{otherwise}}}}\hfill\end{array}\right.$$
(3)

In the above, b refers to the bed elevation, M is defined as (ϕmax − ϕmin)/(bmax − bmin), and (x, y) refers to a given point in space. We use values of 5, 10, and 15 ° for ϕmin in the ice sheet model ensemble, with ϕmax of 40 °, bmin of –300 m, and bmax of 200 m. While we focus on ϕmin in our ice sheet model ensemble, it should be noted that these other values also have high uncertainties, and in particular, bmin should be further explored in relation to ϕmin.

In PISM, Ntill is based on the amount of water contained in the till. If the effective pressure is lower, more of the weight of the ice is carried by the pressurized water in the till, thereby allowing the ice to slide more easily. We use a non-conserving subglacial hydrology scheme, based on the undrained plastic bed model of Tulaczyk et al. (2000)69, which balances the basal melt rate and constant drainage rate to determine the effective pressure on the saturated till; the scheme adapted for PISM is described in detail in Bueler and van Pelt (2015)70, referred to as the "null" hydrology model. We use the default PISM values for the fraction of effective overburden pressure and the decay rate of the effective water content in the till. Given that these parameters influence the effective water content in the till, this should be considered an additional parameter-related uncertainty.

Statistical emulator and historical constraints

Statistical emulation has been used in a number of studies investigating ice sheet sea-level contribution and regional sea-level change6,10,12,22,71. In this case, the statistical emulator is a regression model that is based on GPs. GPR is a non-parametric approach to fitting a function to a set of training data x.

$$f({{{{{{{\bf{x}}}}}}}}) \sim {{{{{{{\mathcal{GP}}}}}}}}\left(m({{{{{{{\bf{x}}}}}}}}),k({{{{{{{\bf{x}}}}}}}},{{{{{{{\bf{x}}}}}}}}^{\prime} )\right)$$
(4)

In this notation, the function f, which is the Antarctic contribution to eustatic sea level at a specific point in time, is a random variable that is fully described by the independent variable x, a mean function, m(x), which we can assume to be zero if we normalize the data (m(x) = 0), and a covariance function \(k({{{{{{{\bf{x}}}}}}}},{{{{{{{\bf{x}}}}}}}}^{\prime} )\). A popular choice of the covariance function k is the radial basis function (RBF) kernel with different length scales for each feature,

$$k({{{{{{{\bf{x}}}}}}}},{{{{{{{\bf{x}}}}}}}}^{\prime} )={\sigma }^{2}\exp \left(-\frac{1}{2}{({{{{{{{\bf{x}}}}}}}}-{{{{{{{\bf{x}}}}}}}}^{\prime} )}^{T}M({{{{{{{\bf{x}}}}}}}}-{{{{{{{\bf{x}}}}}}}}^{\prime} )\right)$$
(5)

with signal variance σ2 and M = diag(l)−2. Here, l is a vector of length-scale parameters that has to be learnt from the data during training. The data itself is a list of set, each set consisting of the seven feature, four PISM parameters, ESSA, ESIA, q, ϕ, and three forcings that depend on the respective RCP scenarios (either RCP2.6 or RCP8.5), and which are global mean temperature (GMT), the cumulative sum of GMT (∑GMT), and the time since last GMT change (\(\hat{t}\)). Those forcings are timeseries that have been split into individual data points, for example GMT(t = 2018), GMT(t = 2019), etc. The data set, i.e., the x in Eqs. (4) and (5), has the following structure:

$$\left\{\begin{array}{l}\{{E}_{SSA}^{1},{E}_{SIA}^{1},{q}^{1},{\phi }^{1},GM{T}^{RCP2.6}({t}_{1}),\sum GM{T}^{RCP2.6}({t}_{1}),{\hat{t}}_{1}^{RCP2.6}\}\\ \{{E}_{SSA}^{1},{E}_{SIA}^{1},{q}^{1},{\phi }^{1},GM{T}^{RCP2.6}({t}_{2}),\sum GM{T}^{RCP2.6}({t}_{2}),{\hat{t}}_{2}^{RCP2.6}\}\\ \ldots \\ \{{E}_{SSA}^{1},{E}_{SIA}^{1},{q}^{1},{\phi }^{1},GM{T}^{RCP2.6}({t}_{k}),\sum GM{T}^{RCP2.6}({t}_{k}),{\hat{t}}_{k}^{RCP2.6}\}\\ \{{E}_{SSA}^{2},{E}_{SIA}^{2},{q}^{2},{\phi }^{2},GM{T}^{RCP2.6}({t}_{1}),\sum GM{T}^{RCP2.6}({t}_{1}),{\hat{t}}_{1}^{RCP2.6}\}\\ \ldots \\ \{{E}_{SSA}^{N},{E}_{SIA}^{N},{q}^{N},{\phi }^{N},GM{T}^{RCP2.6}({t}_{k}),\sum GM{T}^{RCP2.6}({t}_{k}),{\hat{t}}_{k}^{RCP2.6}\}\\ \{{E}_{SSA}^{1},{E}_{SIA}^{1},{q}^{1},{\phi }^{1},GM{T}^{RCP8.5}({t}_{1}),\sum GM{T}^{RCP8.5}({t}_{1}),{\hat{t}}_{1}^{RCP8.5}\}\\ \ldots \\ \{{E}_{SSA}^{N},{E}_{SIA}^{N},{q}^{N},{\phi }^{N},GM{T}^{RCP8.5}({t}_{k}),\sum GM{T}^{RCP8.5}({t}_{k}),{\hat{t}}_{k}^{RCP8.5}\}\end{array}\right\}$$
(6)

with the ID of each PISM parameter ensemble as superscript  {1, …, N = 81} for the first four, and time index tk {2017, 2018, …, 2300} for the remaining three independent variables. Thus, the full data for the two RCP scenarios and the 81 parameter combinations consist of 2 × 81 × (2300 − 2017) = 45846 rows with seven features each. Solving a GPR analytically requires inverting large matrices, which typically scales as \({{{{{{{\mathcal{O}}}}}}}}({n}^{3})\) operations, where n is the number of training data. Fitting a GP, therefore, gets slow for large n. In our case, however, we required only 5% of the data for the training step, i.e., 2292 randomly picked sets from Eq. (6), to fit the GP.

The length scales of the RBF kernel after fitting are

$${{{{{{{\bf{l}}}}}}}}={\left(2.14,0.263,0.231,5.47,3.01,3.41\times 1{0}^{4},2.17\times 1{0}^{3}\right)}^{T}$$

and the signal variance is

$${\sigma }^{2}=0.55{8}^{2}$$

The resulting covariance matrix, as computed using Eq. (5) for the training data, is shown in Supplementary Fig. S6.

In a second step, we predicted the mean of f in Eq. (4) for the whole data set. The prediction error of the GPR model for each of the PISM ensemble members is shown in Supplementary Fig. S7.

For the historical constraints on the model parameters, we test three sets of estimated Antarctic sea-level contributions24,25,26. The emulator, which is trained using the RCP8.5 and RCP2.6 datasets, is hindcast for the past decades (Supplementary Fig. S8). We discount parameter combinations that produce sea-level contribution trajectories that over- or underestimate these historical estimates. The upper and lower bound of the historical sea-level contributions is set to 3 SD, which corresponds approximately to the 99% CI. While each estimation covers nearly the same parameter space for the emulator, the most conservative case (i.e., more parameter space is included) is the IMBIE (2018) estimation (Supplementary Fig. S9). As a result, we use this parameter space in the emulator and ice sheet model analyses. In the ice sheet model ensemble, 7 of the original 81 parameter combinations meet this criterion (Supplementary Fig. S10).

A key caveat of these sea-level projections is that the emulator is trained using ice sheet simulations that use climate forcing from a single AOGCM, the NorESM1-M20. While NorESM1-M exhibits high skill in reproducing historical metrics of the Southern Ocean and Antarctic climate, it does show relatively high oceanic warming in the RCP8.5 scenario compared to other models, though not the highest21. If we consider a more conservative oceanic warming, as is the case with the RCP8.5-CO ensemble (Supplementary Table S3), then emergence between RCP and GWL scenarios may require an even longer timespan than in the standard projections. This is also true for high-emissions scenarios with relatively higher surface mass balance (Supplementary Figs. S14 and S15). The reverse is true considering a less conservative case of oceanic warming (i.e., RCP8.5-WO), in which emergence would occur earlier and the AIS sea-level commitment would be greater.

The global-mean-temperature timeseries for the different scenarios, which are inputs to the emulator, were smoothed using a Savitzky–Golay filter72 to remove the more rapid short-term climate variations (see Supplementary Fig. S5). The Savitzky–Golay filter performs a polynomial least-squares fit within the filter window for the original, unfiltered RCP timeseries. The smoothing can be understood as a weighted moving average filter. The window size is set to 51, and the polynomial degree is 3.