Global warming levels control time of emergence of Antarctic ice loss

Of all the components of the global sea level budget, the future contribution of the Antarctic Ice Sheet (AIS) is the most uncertain in sea level rise projections. Dynamic ice sheet model simulations show considerable overlap in the projected AIS sea level contribution under various greenhouse gas emissions scenarios and the timescale at which scenario-dependence will emerge is unclear. With historically-constrained ice-sheet simulations and a statistical emulator, we demonstrate 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 ﬂow and sliding. However, the results also indicate that the total global warming that occurs over the 21st century controls the resulting long-term AIS sea level commitment, with multi-meter differences between the highest and lowest emissions scenarios in subsequent centuries.

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 (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 period 9,18,19 , here we also consider the temporal evolution of the sea level contribution. In particular, the GP regression uses the 4 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 realised 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 ( Fig S1) 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 modelling alone. Furthermore, we can explore a whole range of different scenarios, for example, additional RCPs, such as RCP4.5 or RCP6.0 (Fig. 1).
We discount parameter combinations that produce sea level contribution trajectories that exceed or underestimate historical estimates of the AIS sea level contribution [20][21][22] (grey shading in Fig 1; Fig S2). 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, E SIA and E SSA ) associated with higher sea level contributions and vice versa. While the full range of projected sea level contributions at 2100 exceeds those of ISMIP6 for the two RCP scenarios 7 , considering only the historically constrained parameter combinations eliminates the most extreme scenarios, thereby narrowing the projected ranges to 0.12 to 0.44 m for RCP2.6 and 0.21 to 0.56 m for RCP8.5 at 2100 relative to 2000 (95% confidence interval, blue and orange shading in Fig 1, respectively). The longer-term response to 2300 results in ranges of 0.45 to 1.57 m and 1.96 to 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% confidence interval (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 (Table S2). 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 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 (standard deviation > 100 m) of the ASE and Wilkes Land outlet glaciers.
Accounting for the high variance in regional Southern Ocean temperatures among RCP8.5 climate projections of different AOGCMs of CMIP5 17 , 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 standard deviations in ice thickness for the full ensemble (Fig 2g,h).
The difference in ice shelf thickness between the scenarios exceeds the standard deviation of a control ensemble forced with a constant modern climate (defined here as threshold exceedance) 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  3i-l). This 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. 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. In the decades prior to the time of emergence (95% CI of 2189), the difference in sea level contribution between the RCP scenario ensemble averages increases most in the Ross and Weddell Sectors (Fig 3). In particular, the Ross Ice Shelf shows marine ice formation and ice sheet advance in RCP2.6, but high melt rates and ice shelf disintegration by 2100 in RCP8.5. The difference between melt rates of the Ronne-Filchner Ice Shelf in the Weddell Sector is less extreme, but the RCP8.5 does show partial ice shelf collapse by 2150. Although both scenarios show relatively high basal melt rates of ASE, the difference in sea level contribution increases over this time period. 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 ( Fig S5).

Uncertainty in the rate of marine ice sheet instability
Although the ISMIP6 projections extend to 2100, this period is insufficient for identifying scenario dependence of the AIS sea level contribution considering the influence of unconstrained model parameters. The region of highest uncertainty during this period is the ASE (Fig 2e-h), where AIS mass loss is currently highest due to high rates of sub-ice shelf melting 21,28 . Importantly, the Thwaites and Pine Island Glaciers in the ASE overlie a deep marine basin with downward-sloping basal topography that lies substantially below sea level (Fig 4a), which creates an ice sheet configuration prone for marine ice sheet instability (MISI) 29 . While ice sheet models suggest that continued ocean warming will lead to ice shelf thinning and grounding-line retreat in this sector 30,31 , basal friction and ice rheology influence the rate of MISI 12,13 , which can have a large impact on the sector's sea level contribution 32 . The historically constrained ice sheet model simulations that lead to the lowest and highest AIS sea level contribution (referred to as min and max in Fig 4, respectively) indicate that by the end of the century (2090-2100), the parameter combination is more influential than the emissions scenario in terms of ice sheet dynamics and sea level contribution (Fig 4c-f). In all cases, the Thwaites and Pine Island Glaciers show increased ice surface velocity relative to the initial condition. 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 4e,f to Fig 4c,d). Notably, less grounding-line retreat of the Thwaites Glacier occurs with the model parameter combination that leads to the highest sea level contribution.
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 4g), when RCP8.5 ocean temperatures increase substantially (Fig S6). The difference in the ASE sea level contribution increases rapidly in the following decades between emissions scenarios, particularly after 2100 (Fig 3a). For the entire ice sheet, the RCP8.5 forcing also produces an increase in the rate of sea level contribution in the last decades of the century. This increase is more significant for the model parameter combination with the highest sea level contribution (i.e., high enhancement factors, low q) 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 is sensitive to the parameter combination. The overlap between RCP scenarios with parameter combinations with the lowest and highest sea level contributions persists until after 2150, just prior to the time of emergence indicated by the emulator (95% CI of 2189).

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 • C 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 5a). Despite these substantial differences, the range of the projected AIS sea level contribution is nearly indistinguishable between these scenarios at 2100 (Fig 5b). 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.
In considering that the Paris Agreement sets a target of 2 • C of global warming, we explore how 2 • C of warming by the end of the century compares to different warming magnitudes (Fig 5c,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 5c,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 5e,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 temperature 33,34 , precipitation 35 , sea surface temperature and ocean carbon uptake 36 , are expected to permanently exceed recent natural variability within decades due to anthropogenic influence. The projections analysed 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 decades 37,38 .
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 CMIP5 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 Seas 17 . It has been estimated that the Southern Ocean accounts for about 75% of global ocean heat uptake 39 , and the combined effects of wind and meltwater may be causing local warming to accelerate by increasing upwelling of warm, saline circumpolar deepwater 40 . Meltwater feedback, not considered in the CMIP5 modesl used here or CMIP6 models 41 , has also been suggested to act as a positive feedback that increases AIS ice loss 42 . 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 21st century controls the resulting AIS sea level contribution, independent of the rate of warming. Multi-meter differences in the total contribution over centennial time scales are projected between 2 and 4 • C of global warming. However, despite recent improvements in mapping bed topography 43 , the limited constraints on substrate conditions and ice rheology, which impact ice flow and sliding, are especially important over decadal timescales in regions prone to MISI, such as the ASE and Wilkes Land. 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 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 time scales under GWLs >2 • C. The disintegration of these largest AIS ice shelves under these high

7/13
GWLs ultimately results in widespread retreat of WAIS, which is likely not reversible on human timescales 44 , and a sea level commitment that is multiple meters higher than the potential range under a GWL limited to 2 • C.

Ice sheet model simulations
We perform ice sheet simulations using the Parallel Ice Sheet Model (PISM) version 1.1, a three-dimensional thermo-mechanical ice sheet-ice shelf model 45 . 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 ice 46 . 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 field 47 , and can migrate freely. We do not include a sub-grid basal melt scheme for partially floating cells 23 , however, as this may overestimate the rate of grounding line retreat 48 . We employ an ocean model based on Holland and Jenkins (1999) 49 to compute basal ice shelf melt rates and temperature from thermodynamics in an ocean boundary layer.
To initialize the ice sheet model, we use an identical approach to the Victoria University of Wellington (VUW) contribution to ISMIP6 7 . Starting from initial bedrock and ice thickness conditions from Morlighem et al. (2020) 43 and reference climatology from Van Wessem et al. (2014) 50 , 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) 42 , 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 shallow-ice approximation 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 .
The projection experiments are initiated in the year 2015 from the thermally-equilibrated state after again resetting the ice thicknesses and bed topography. Simulations are run for 285 years at 16 km horizontal resolution. For the main RCP2.6 and RCP8.5 model ensembles, we use climate forcing derived from the NorESM1-M CMIP5 model 16 , the only Tier 1 model of ISMIP6 that included both RCP scenarios 7,17 . 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 cavities 51 . 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. Considering the substantial differences in regional ocean temperature changes in the Southern Ocean of CMIP5 models 17 , we also perform separate parameter-constrained ensembles using climate forcing from the other ISMIP6 Tier 1 models with the relatively lowest and highest ocean temperature anomalies, CCSM4 RCP8.5 (RCP8.5-CO) and HadGEM2-ES RCP8.5 (RCP8.5-WO), 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 stress 52 . These specific parameters were selected given their known influence in ice sheet model simulations, yet poor constraints 10 , and are described in more detail below. We use three values per four model parameters for a total of 81 combinations (see Table S1).
The parameter values cover a wide range of values previously used in PISM simulations of the Antarctic Ice Sheet 10,14,42,45,53 . Comparison of the simulations to observed Antarctic mass loss of the last decade 21 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.

Ice sheet model parameters
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 uniformly 45 . 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 E SIA > 1 and E SSA < 1 54,55 . Increasing E SIA allows the ice to deform more easily, and increasing E SSA produces faster ice streams and thinner ice shelves. Based on enhancement factors previously

8/13
used for AIS simulations with PISM 10,14,42,45,53  PISM parameterizes ice sheet sliding in the form of a power law that ranges from plastic Coloumb sliding to a linear sliding law 46,52 : In the power law, τ b is the basal shear stress, τ c is the yield stress, and u is velocity, where u threshold is a threshold velocity. 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 φ : where c 0 is the till cohesion, N till is the effective pressure of the till, and φ is the till friction angle. As in previous studies 10,45,46,56 , φ 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 weaker 57 : In the above, b refers to the bed elevation, M is defined as (φ max − φ min )/(b max − b min ), 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 • , b min of -300 m, and b max of 200m.

Statistical Emulator
Statistical emulation has been used in a number of studies investigating ice sheet sea level contribution and regional sea level change 6,9,18,19,58 . In this case, the statistical emulator is a regression model that is based on Gaussian Processes (GP). GP regression (GPR) is a non-parametric approach to fitting a function to a set of training data x.
In this notation, the function f , which is the 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(x, x ). A popular choice of the covariance function k is the radial basis function (RBF) kernel with different length scales for each feature, 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 7 feature, four PISM parameters, E SSA , E SIA , 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 (t). Those forcings are time series that have been split into individual data points, for example GMT (t = 2018), GMT (t = 2019), etc. The data set, i.e., the 9/13 x in Eq. (4) and (5), has the following structure: SSA , E 1 SIA , q 1 , φ 1 , GMT RCP2.6 (t 1 ), ∑ GMT RCP2.6 (t 1 ),t RCP2.6 1 } {E 1 SSA , E 1 SIA , q 1 , φ 1 , GMT RCP2.6 (t 2 ), ∑ GMT RCP2.6 (t 2 ),t RCP2.6 2 } . . .
The length scales of the RBF kernel after fitting are l = 2.14, 0.263, 0.231, 5.47, 3.01, 3.41 × 10 4 , 2.17 × 10 3 T and the signal variance is The resulting covariance matrix, as computed using Eq. (5) for the training data is shown in Fig S3. 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 Fig S1. 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-M 16 . 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 highest 17 . If we consider a more conservative oceanic warming, as is the case with the RCP8.5-CO ensemble (Table S3), then emergence between RCP and GWL scenarios may require an even longer timespan than in the standard projections. 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.