ISMIP6 projections of ocean-forced Antarctic Ice Sheet evolution using the Community Ice Sheet Model

The future retreat rate for marine-based regions of the Antarctic Ice Sheet is one of the largest uncertainties in sealevel projections. The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) aims to improve projections and quantify uncertainties by running an ensemble of ice sheet models with atmosphere and ocean forcing derived from global climate models. Here, ISMIP6 projections of ocean-forced Antarctic Ice Sheet evolution are illustrated using the Community Ice Sheet Model (CISM). Using multiple combinations of sub-ice-shelf melt parameterizations and calibrations, CISM is spun up to 5 steady state over many millennia. During the spin-up, basal friction parameters and basin-scale thermal forcing corrections are adjusted to nudge the ice thickness toward observed values. The model is then run forward for 500 years, applying ocean thermal forcing anomalies from six climate models. In all simulations, the ocean forcing triggers long-term retreat of the West Antarctic Ice Sheet, including the Amundsen, Filchner-Ronne, and Ross Basins. Mass loss accelerates late in the 21 century and rises steadily over the next several centuries without leveling off. The resulting ocean-forced SLR at year 2500 varies from 10 about 10 cm to nearly 2 m, depending on the melt scheme and model forcing. Relatively little ice loss is simulated in East Antarctica. Large uncertainties remain, as a result of parameterized basal melt rates, missing ocean and ice sheet physics, and the lack of ice–ocean coupling.

Rather than make quantitative predictions (which remain elusive because of model and forcing uncertainties), our main intent is to explore parameter space, taking advantage of the ISMIP6 framework to build on previous multi-century Antarctic simulations (e.g., Pollard and Deconto, 2009;Cornford et al., 2015;Pollard and DeConto, 2016;Larour et al., 2019). The ice-sheet physics is conventional in the sense that it includes well-understood retreat mechanisms such as MISI, but not hydrofracture or cliff collapse (Pollard et al., 2015;Pollard and DeConto, 2016). Some important long-term feedbacks such as solid-Earth 5 rebound and relative-sea-level changes (Gomez et al., 2010;Larour et al., 2019) are omitted, and there is no ice-ocean coupling (e.g., De Rydt and Gudmundsson, 2016;Seroussi et al., 2017;Favier et al., 2019). Thus, the timing and magnitude of simulated ice sheet retreat are imprecise, but we can identify responses that are robust across simulations, and also draw attention to the largest sources of uncertainty.
Section 2 gives an overview of CISM and summarizes the protocol and ocean data sets for ISMIP6 Antarctic projections. We The Community Ice Sheet Model (CISM) is a parallel, higher-order ice sheet model designed to perform continental-scale 15 simulations on timescales of decades to millennia. CISM is a descendant of the Glimmer model (Rutt et al., 2009) and is now an open-source code developed mainly at the National Center for Atmospheric Research, where it serves as the dynamic ice sheet component of the Community Earth System Model (CESM). The most recent documented release was CISM v2.1 , coinciding with the 2018 release of CESM2 (Danabasoglu et al., in review). The model performs well for community benchmark experiments, including the ISMIP-HOM experiments for higher-order models (Pattyn et al.,20 2008) and several stages of the Marine Ice Sheet Model Intercomparison Project: the original MISMIP (Pattyn et al., 2012), MISMIP3d (Pattyn et al., 2013), and MISMIP+ (Asay- Davis et al., 2016;Cornford et al., in review). CISM participated in two earlier ISMIP6 projects focused on ice sheet model initialization: initMIP-Greenland (Goelzer et al., 2018) and initMIP-Antarctica . More recently, CISM results were submitted for the ISMIP6 projections (Goelzer et al., in review; Seroussi et al., in review), the LARMIP-2 experiments (Levermann et al., 2019), and the Antarctic BUttressing Model 25 Intercomparison Project (ABUMIP; Sun et al., in review).
CISM runs on a structured rectangular grid with a terrain-following vertical coordinate. The simulations in this paper were run on a 4-km grid with 5 vertical levels, as for the CISM contributions to ISMIP6 projections (Seroussi et al., in review).
Scalars (e.g., ice thickness H and temperature T ) are located at grid cell centers, with horizontal velocity u = (u, v) computed at vertices. The dynamical core has parallel solvers for a hierarchy of approximations of the Stokes ice-flow equations, including 30 the shallow-shelf approximation (MacAyeal, 1989), a depth-integrated higher-order approximation (Goldberg, 2011), and the 3D Blatter-Pattyn higher-order approximation (Blatter, 1995;Pattyn, 2003). (The latter two approximations are classified as L1L2 and LMLa, respectively, in the terminology of Hindmarsh, 2004.) The simulations for this paper use the depth-integrated solver, which gives a good balance between accuracy and efficiency for continental-scale simulations .
CISM supports several basal sliding laws. For ISMIP6 projections we use a sliding law based on Schoof (2005), with a functional form suggested by Asay-Davis et al. (2016): 5 where τ b is the basal shear stress, u is the basal ice velocity, N is the effective pressure, m = 3 is a power-law exponent, and C p and C c are empirical coefficients. In the ice sheet interior, where the ice is relatively slow-moving with large effective pressure, this law gives power-law behavior: Where the ice is fast-moving with low effective pressure, we have Coulomb behavior: Since the power-law coefficient C p is spatially variable and poorly constrained, we adjust C p (x, y) to nudge the ice thickness toward observed values, as described in section 3.1. The dimensionless Coulomb coefficient is set to C c = 0.5, close to the value of 0.42 used by Pimentel et al. (2010). The effective pressure is set to the ice overburden pressure, N = ρ i gH, where ρ i is ice density and g is gravitational acceleration. We initially made N proportional to (1 − H f /H), where H f is the flotation 15 thickness, to represent the connection of subglacial water to the ocean (Leguy et al., 2014), but this choice led to excessive grounding-line retreat in some regions. Setting N to overburden pressure implies power-law behavior in nearly all of the ice sheet.
CISM also supports several calving laws, but none was found to give calving fronts in good agreement with observations for both large and small Antarctic ice shelves. Instead, we use a no-advance calving mask, removing all ice that flows beyond the 20 observed calving front. The calving front can retreat where there is more surface and basal melting than advective inflow, but more often the calving front remains in place.
Many studies (e.g., Pattyn, 2006;Schoof, 2007) have emphasized the challenges of simulating ice dynamics in the transition zone between grounded and floating ice, especially when the grid resolution is ∼ 1 km or coarser, as is the case for most simulations of whole ice sheets. Grounding line parameterizations (GLPs), which give a smooth transition in basal shear 25 stress across the transition zone, have been shown to improve numerical accuracy in models with relatively coarse grids (e.g., Gladstone et al., 2010;Leguy et al., 2014;Seroussi et al., 2014). CISM has a GLP that was not used for the Greenland simulations in Lipscomb et al. (2019), but is used for Antarctic simulations and is described in Appendix A.

Protocols based on ISMIP6 for Antarctic projections
Protocols for the ISMIP6 Antarctic projections are described in detail by Nowicki et al. (in review) and Jourdain et al. (in review). Here, we give a brief summary, noting where the simulations in this paper differ from the protocols. The ISMIP6 projection experiments are run with standalone ice sheet models, forced by time-varying atmosphere and ocean fields derived from the output of CMIP5 and CMIP6 AOGCMs. (The simulations described here use only the ocean forcing.) The CMIP5 5 models were selected by a procedure described by Barthel et al. (2019) to sample the diversity in climate evolution around Antarctica. For CMIP6 models, there was not time for a formal selection, but output from selected models was processed as it became available. Each experiment runs for 86 years, from the start of 2015 to the end of 2100. Model initialization methods are left to the discretion of each group, and are detailed for specific models in Seroussi et al. (2019, in review). If the initialization date is before the start of the projections, a short historical run is needed to advance the ISM to the end of 2014.

10
For projections, ISMIP6 provides data sets of yearly atmosphere and ocean forcing on a standard 8-km grid. The ocean forcing consists of 3D fields of thermal forcing (i.e., the difference between the in situ ocean temperature and the in situ freezing temperature), obtained by a complex process described by Jourdain et al. (in review). It is not possible to use ocean temperature and salinity directly from AOGCMs, because current CMIP models do not simulate ocean properties in ice shelf cavities. Moreover, global ocean models are usually run at a resolution of ∼ 1 • , too coarse to give an accurate mean state for 15 ocean properties near Antarctic ice shelves. Instead, Jourdain et al. (in review) combined recent data sets (Locarnini et al., 2019;Zweng et al., 2019;Good et al., 2013;Treasure et al., 2017) to construct a 3D gridded climatology of ocean temperature and salinity. This climatology was interpolated to fill gaps and extrapolated into ice-shelf cavities. To obtain forcing fields for projections, temperature and salinity from the various AOGCM ocean models were extrapolated to cavities, and their anomalies were added to the observationally derived climatology. The result is a time-varying 3D product that can be vertically 20 interpolated to give the thermal forcing at the base of a dynamic ice shelf at any time and location in the Antarctic domain.
To compute basal melt rates beneath ice shelves, ISMIP6 models can use a standard approach, an open approach, or both.
The open approach is chosen independently by each modeling group; the only requirement is to use the ocean data provided by Jourdain et al. (in review). For the standard approach, basal melt rates beneath ice shelves are computed as a quadratic function of thermal forcing as suggested by Favier et al. (2019), with a thermal forcing correction suggested by Jourdain et al.

25
(in review): where m is the melt rate, T F is the thermal forcing, γ 0 is an empirical coefficient, ρ w and c pw are the density and specific heat of seawater, and ρ i and L f are the density and latent heat of melting of ice. The brackets in < T F > denote the average over a drainage basin or sector, and δT sector is a thermal forcing correction with units of temperature, with one value per An alternative local parameterization is obtained by replacing < T F > in the last term: giving a quadratic dependence on the local thermal forcing. A third parameterization, which we call nonlocal-slope or simply slope, is obtained by multiplying Eq. 4 by sin(θ), where θ is the local angle between the ice-shelf base and the horizontal. At the same time, γ 0 is increased by a factor of ∼ 100, since sin(θ) typically is ∼ 10 −2 . This scheme, suggested by Little et al.

5
(2009) and Jenkins et al. (2018), generally gives larger melt rates near grounding lines and lower melt rates near the calving front, in agreement with observational estimates and ocean simulations (e.g., Rignot et al., 2013;Favier et al., 2019). Total melt over the ice sheet is the same, if γ 0 is tuned based on a total-melt criterion. CISM simulations in Seroussi et al. (in review) used the slope parameterization as an open approach.
Sub-ice-shelf thermal forcing and melt rates are uncertain, as is the functional relationship between thermal forcing and The simulations in our paper differ from the ISMIP6 protocols in the treatment of δT sector . Although we use calibrated 5 values of γ 0 , we tune δT sector to better match the mean ice thickness near the grounding line, as described in Sect. 3.1. This procedure calibrates δT sector to optimize melt rates near the grounding line where they matter most for ice extent, but also gives basin-average melt rates that differ from the observational estimates.
In summary, we have presented three basal melt parameterizations (local, nonlocal, and slope) and two calibration methods  3 Model spin-up

Nudging procedure
When spinning up the model, we aim to to reach a stable state with minimal drift under modern forcing, while also simulating 15 ice sheet properties that agree with observations. It is challenging to achieve both goals at once. Long spin-ups typically yield a steady-state ice sheet with large biases in thickness, velocity, and/or ice extent, while initialization methods that assimilate data to match present-day conditions often have a large initial transient. We compromised by using a hybrid method similar to that of Pollard and DeConto (2012), running a long spin-up while continually nudging the ice sheet toward the observed thickness.
Beneath grounded ice, we adjust C p (x, y), a poorly constrained, spatially varying friction coefficient in the basal sliding 20 law (Eq. 1). This coefficient controls the power-law behavior (Eq. 2) in most of the ice sheet, with higher C p giving greater friction and less sliding. Wherever grounded ice is present, C p is initialized to 50, 000 Pa m −1/3 y −1 . During the spin-up, C p is decreased where H > H obs and increased where H < H obs , based on the idea that lower friction will accelerate the ice and lower the surface, while higher friction will slow the ice and raise the surface. The rate of change of C p is given by where H obs is an observational target, H 0 = 100 m is a thickness scale, and τ c = 500 y is a time scale for adjusting C p . Eq. 6 is based on the equation for a critically damped harmonic oscillator, where the first term in brackets nudges H toward H obs , and the second term damps the nudging to prevent overshoots. (The damping is not exactly critical, however, because dC p /dt 5 is not exactly proportional to d 2 H/dt 2 .) We hold C p within a range between 10 2 and 10 5 Pa m −1/3 y −1 , since smaller values can lead to numerical instability, and larger values do not significantly lower the sliding speed.
This method works well in keeping most of the grounded ice near the observed thickness. Also, since C p is independent of the ice thermal state, we remove low-frequency oscillations associated with slow changes in basal temperature, resulting in a better-defined steady state. In forward runs, however, C p (x, y) is held fixed and cannot evolve in response to changes in basal 10 temperature or hydrology. As a result, we can have unphysical basal velocities when the ice dynamics differs from the spun-up state.
Another nudging method is needed to initialize floating ice shelves, where C p = 0. In CISM simulations for initMIP-Antarctica and ISMIP6 projections, (Seroussi et al., 2019, in review), basal melt rates were obtained by nudging. That is, an equation similar to Eq. 6 was used to nudge H toward H obs by adjusting the melt rate m(x, y) in each floating grid cell. In 15 climate change experiments, the spun-up melt rates were added to prescribed melt rate anomalies in each basin. Although this method yields ice-shelf thicknesses and grounding-line locations that agree well with observations, it also overfits the observations, giving noisy melt rates that compensate for other errors without being tied to ocean temperatures. More complications arise when applying this spin-up method to the ISMIP6 projections, which prescribe thermal forcing anomalies instead of melt rate anomalies. In climate change experiments, a melt rate anomaly computed from the thermal forcing anomaly must be added 20 to the spun-up melt rate, instead of computing the evolving melt rate directly from the evolving thermal forcing.
For the simulations described here (which were done too late to be included in Seroussi et al., in review), we take a different approach. During the spin-up, basal melt rates are computed directly from the thermal forcing, using the climatological data set and melt parameterizations described in Sect. 2.2. When we use the calibrated values of both γ 0 and δT sector , many grounding lines drift far from their observed locations. The drift can be reduced (but not eliminated) by continually adjusting δT sector 25 in each of 16 sectors (see Fig. 1), nudging toward an ice thickness target in a region near the grounding line. Here, "near the grounding line" is defined as having a value of a function f float whose magnitude is less than a prescribed value: where b is the seafloor elevation (negative below sea level) and H thresh is the prescribed threshold thickness. Where f float > 0, it is equal to the thickness of the ocean cavity, and where f float < 0, the ice is grounded. We set H thresh = 400 m, so that the 30 target region includes most of the ice likely to switch between floating and grounded during a spin-up, without extending too far upstream. During the spin-up, δT in each sector is adjusted as follows: whereH is the mean ice thickness over the target region,H obs is the observational target for this region, m T is a scale for the rate of change of basal melt rate m with temperature T , and τ m is a timescale for adjusting δT . As in Eq. 6, the firstderivative term damps oscillations. After some experimentation, we set m T = 10 m y −1 K −1 and τ m = 10 y. In most basins, 5 this adjustment keepsH close toH obs with modest values of δT (∼ 1 K or less). In some basins,H <H obs no matter how much δT is lowered, and in these basins δT is capped at −2 K.
Given the calibrated thermal forcing, the melt rate m is computed for floating grid cells using one of the three parameteriza- Running CISM for several millennia while adjusting for C p in each grounded grid cell and δT in each basin, the ice sheet approaches a steady state. We evaluate the spun-up state below.

Spun-up model state
We carried out six Antarctic spin-ups, each with a different pairing of the three basal melt parameterizations (local, nonlocal, 15 and nonlocal-slope) and the two calibrations (MeanAnt and PIGL). For each spin-up, the ice sheet is initialized to presentday thickness using the BedMachineAntarctica data set  with analytic vertical temperature profiles.
The surface mass balance and surface temperature are provided by a regional climate model, RACMO2.3 (van Wessem et al., 2018), and the geothermal heat flux is from Shapiro and Ritzwoller (2004). During each spin-up, the simulated ice thickness is nudged toward observations by adjusting C p (x, y) and δT sector as described in Sect. 3.1. The spin-ups were run for 20,000 20 model years, allowing the ice sheet to come very close to steady state; at this time the total mass is changing by < 1 Gt y −1 .
Although γ 0 varies widely among the different spin-ups, the spun-up states are similar across parameterizations and calibrations, because of the freedom to adjust C p and δT sector independently for each run. The simulated thickness, velocity, and ice extent are broadly in agreement with observations, but with some persistent biases. Some biases can likely be attributed to errors in ocean thermal forcing (which is treated simply by the basin-scale melt parameterizations) and seafloor topography 25 (e.g., an absence of pinning points, resulting in grounding-line retreat that can be compensated by spurious ocean cooling).    . Overall, the agreement with observations is very good for both grounded and floating ice, even though the model is nudged toward observed velocities only indirectly, via the thickness field. In general, we find that good agreement in thickness implies good agreement in velocity as well. (This is true, at least, when using BedMachineAntarctica thicknesses, which are obtained using the mass conservation method of Morlighem et al., 2011.) One place of disagreement is the Kamb Ice Stream on 5 the Siple Coast, which is present in the model but absent in the observations, having stagnated in the 1800s but left a signature in the thickness field (Ng and Conway, 2004).
For the same spin-up, Fig. 4 shows ice surface speeds in the Amundsen Sea sector, including Pine Island and Thwaites Glaciers. The observations reveal a dual structure in the Thwaites velocity field, with a fast western core where the glacier flows into the Thwaites Ice Tongue, and slower speeds in the east where flow is impeded by an ice rise (cf. Fig. 1 in Rignot 10 et al., 2014). The model, however, lacks a sharp division between east and west, and the Thwaites grounding line is retreated compared to observations, perhaps because the interaction with seafloor topography is not captured correctly. At the same time, the Pine Island grounding line is advanced compared to observations. The grounding lines of both glaciers have been retreating since at least the 1990s (Rignot et al., 2014), and the spin-up method is not well suited to initialize dynamic grounding lines in their observed locations.   δT is capped at −2 K, meaning that ice near the grounding line is biased thin, even when basal melting is reduced to zero.
Also, δT is strongly negative for the Amundsen Sea sector, with values < −1 K for the PIGL calibration. This is the model's attempt to curtail the grounding-line retreat that occurs, especially for Thwaites Glacier, when γ 0 is large and δT is near 0.
This means that melt rates might be underestimated in the spin-up, which could alter the sensitivity to future warming. A goal for future work is to reduce δT for these sectors through more realistic ocean forcing and/or seafloor topography. Finally, we note that although the ice sheet mass budget is nearly in balance at the end of the spin-ups, total basal melt rates are 500-600 Gt y −1 , about half the values estimated by Depoorter et al. (2013) and Rignot et al. (2013). Lower-than-observed 12 https://doi.org/10.5194/tc-2019-334 Preprint. Discussion started: 14 January 2020 c Author(s) 2020. CC BY 4.0 License. melt rates are caused mainly by negative δT corrections, again suggesting that δT could be compensating for other model errors, such as too-fast flow in ice shelves that are too lightly buttressed.

Results of projection experiments
After finishing the six spin-ups described in Sect. 3, we ran 36 projection experiments with ocean thermal forcing from six CMIP AOGCMs. The four CMIP5 models are CCSM4, HadGEM2-ES (hereafter HadGEM2), MIROC-ESM-CHEM (hereafter 5 MIROC), and NorESM1-M (hereafter NorESM1). These were among the top six CMIP5 models chosen by Barthel et al. (2019). CCSM4, MIROC, and NorESM1 were used for the ISMIP6 Tier 1 experiments, and we added HadGEM2 (a Tier 2 selection) to sample a model with relatively high ocean warming. The two CMIP6 models are CESM2 and UKESM, which are successors of CCSM4 and HadGEM2, respectively. All the AOGCMs ran high-end emissions scenarios: RCP8.5 for the CMIP5 models and SSP5-85 for the CMIP6 models. Each spun-up ice sheet state was run forward with ocean forcing from the 10 six AOGCMs. We ran each projection experiment for 500 years starting in 1995, the nominal date of the spin-up. In the ISMIP6 protocols, the first 20 years (to the start of 2015) constitute a historical run and the remainder a projection run, but for our simulations the historical and projection periods were forced in the same way. From 1995-2100, we applied the yearly thermal forcing provided 20 by ISMIP6. After 2100, we cycled repeatedly through the last 20 years of forcing, 2081-2100, to evaluate the committed SLR associated with a late-21 st century climate. For comparison, we ran a subset of forcing experiments using a fixed climatology, computed as the mean of the 2081-2100 thermal forcing. Cycling through the yearly forcing drives greater mass loss (by ∼ 15%) than does the fixed climatology, suggesting that years with high thermal forcing have a disproportional influence on long-term mass loss, given the quadratic relationship between thermal forcing and melt rates. Several patterns emerge. First, SLR starts slowly for all experiments, then accelerates near the end of the 21 st century, 30 suggesting that a threshold has been crossed based on some magnitude or duration of thermal forcing. After 2100, SLR is fairly linear and shows no sign of leveling off after 500 years. This is consistent with retreat driven by MISI in large reverse- sloping basins. Once the retreat is under way, it continues until reaching a stable seafloor configuration, up to several hundreds of km upstream.
Second, the ice sheet is more sensitive to some melt combinations than others. Sensitivity to local versus nonlocal parameterizations is similar. However, the nonlocal-slope parameterization yields more SLR than local and nonlocal, and the PIGL calibration gives more SLR than MeanAnt. The greater sensitivity for the slope scheme can be attributed to larger melt rates 5 at steep slopes near grounding lines, where melting is most effective in driving ice retreat. The greater sensitivity of PIGL than MeanAnt results from the higher values of γ 0 . During the spin-ups, PIGL runs acquire more negative values of δT than MeanAnt, to compensate for greater γ 0 . Then, for a given thermal forcing anomaly during the forward runs, the rate of change of m (as given by Eqs. 4 and 5) is proportional to γ 0 and thus is larger for PIGL.
Third, the CMIP model rankings are consistent across melt schemes. HadGEM2 is the warmest model and drives the most 10 SLR, followed by UKESM and then NorESM1. CCSM4 and CESM2 give similar rates of SLR, and MIROC yields the least.
This ranking follows the magnitude of the thermal forcing in the various WAIS basins (Fig. 6). For a given model, the most sensitive melt scheme (slope-PIGL) yields up to 5 times as much SLR as the least sensitive schemes (local-and nonlocal- MeanAnt), and the warmest model (HadGEM2) drives 2 to 4 times as much SLR as the coolest (MIROC). Thus, the total SLR over 500 years varies by a factor of about 20, from ∼ 10 cm to almost 2 m. Fig. 8 shows spatial plots of thickness changes from a subset of projection runs. These plots are similar to Fig. 2, except that the field shown is not the difference in ice thickness, but rather the difference in thickness above flotation-in other words, the part of the ice column that contributes to SLR. The figure shows results from projection runs for three AOGCMs (HadGEM2, 5 NorESM1, and CESM2) and for nonlocal-MeanAnt and nonlocal-slope-PIGL, the least and most sensitive melt combinations. Most of the SLR contribution comes from the large WAIS basins on reverse-sloping beds: Amundsen Sea, Filchner-Ronne, and Ross. In most cases the Amundsen Sea contribution is smaller than the Filcher-Ronne and/or Ross contributions, which is unexpected given that the Amundsen Sea sector dominates current WAIS retreat. It is possible that the melt parameterizations underestimate the delivery of heat to Amundsen Sea grounding lines, in part because of the large negative δT corrections 10 discussed in Sect. 3.2. Conversely, the extrapolation procedure and melt parameterizations might overestimate heat delivery to the Filchner-Ronne and Ross cavities, which currently have little basal melting. Although some ocean models with CMIP3 atmospheric forcing have projected Weddell Sea warming (e.g., Hellmer et al., 2012;Timmermann and Hellmer, 2013), it is not clear that these large shelves will readily transition from cold, low-melt regimes to warm, high-melt regimes (Naughten et al., 2018). Compared to the WAIS, there is a relatively small SLR contribution from East Antarctica, apart from the Amery Glacier region in the slope-PIGL runs. .

Conclusions
Using the Community Ice Sheet Model, we ran an ensemble of ice sheet simulations based on the protocols for the ISMIP6 From the spun-up states, we ran each ensemble member forward for 500 years, applying ISMIP6 ocean thermal forcing for 1995-2100 and cycling repeatedly through the 2081-2100 forcing for the rest of the simulation. Using forcing from four CMIP5 models and two CMIP6 models, we ran 36 projection experiments. For each simulation we analyzed the ice mass 15 loss and associated sea level rise. In all cases, the Antarctic Ice Sheet loses mass, with total losses ranging from about 10 cm to nearly 2 m SLE. Mass loss begins slowly, accelerates in the late 21 st century, and continues steadily for the next several centuries without leveling off. The thinning of grounded ice (i.e., ice that contributes to SLR) is concentrated in the Amundsen, Filchner-Ronne, and Ross Basins of WAIS, all of which have reverse-sloping beds and are vulnerable to marine ice sheet instability. East Antarctica loses relatively little grounded ice. 20 These simulations are consistent with recent studies (e.g., Cornford et al., 2015;Pollard and DeConto, 2016;Larour et al., 2019) showing potential WAIS collapse, driven by ocean warming in regions with reverse-sloping beds. The novelty of this study lies in the application of the ISMIP6 Antarctic forcing protocols (with melt rates derived from the new thermal forcing data of Jourdain et al., in review) to an ensemble of multi-century simulations. The results suggest that, by the end of this century (assuming high-end emissions scenarios), sub-ice-shelf melt rates could be large enough, if sustained, to drive the 25 steady and possibly irreversible retreat of much of WAIS.
Ice loss is greater with the slope parameterization than with the local and nonlocal parameterizations, and greater with the PIGL calibration (which has higher values of the melt parameter γ 0 ) than with MeanAnt. Mass loss varies by a factor of about 5 between the most sensitive melt scheme (slope-PIGL) and the least sensitive (local-and nonlocal-MeanAnt), even when the ocean forcing comes from the same AOGCM. All the melt schemes leave out important physics, and it is not known whether 30 one parameterization or calibration is more accurate than the others. Similarly, ice loss varies by a factor of 2 to 4 between the warmest AOGCM (HadGEM2) and the coolest (MIROC), even when using the same melt scheme. Thus, while suggesting that much of WAIS is vulnerable to projected ocean warming, these results do not place firm bounds on the rate of future SLR.
Ice sheet retreat in these simulations is modest for the next several decades. It is not clear to what extent this inertia is real, rather than an artifact of the spin-up procedure. The model starts close to steady state, nominally in 1995, with no attempt to assimilate retreat that was already under way (e.g., Rignot et al., 2019). Several ISMs in the ISMIP6 Antarctic ensemble (Seroussi et al., in review) use similar spin-up methods that could underestimate near-term ice loss. (Other models use data assimilation techniques that are better equipped to capture ongoing trends.) Given this inertia, Antarctic projections ending 5 in 2100 (largely for practical reasons, since CMIP forcing data often is unavailable after 2100) do not give a full picture of sea-level commitment from ice sheet retreat. It is possible that ocean warming during this century could result in a commitment of several meters from WAIS, but with most of the SLR taking place in future centuries.
To these conclusions, we add several caveats. First, the inversion procedure gives large and possibly spurious negative temperature corrections for the Amundsen Sea basin, perhaps compensating for the absence of pinning points (e.g., for Thwaites 10 Glacier). For this reason, Amundsen Sea retreat-which already is proceeding vigorously-might be underestimated.
Second, the ocean data extrapolation transfers AOGCM heat anomalies from the open ocean to the distant grounding lines of the Filchner-Ronne and Ross shelves, where present-day melting is minimal. To the extent that the simple melt schemes deliver too much heat to grounding lines, the Filchner-Ronne and Ross retreat would be overestimated.
More generally, these simulations leave out many physical processes and feedbacks. Changes in atmospheric forcing are 15 omitted by design, to simplify the analysis; increased snowfall in the ISMIP6 multi-model ensemble can mitigate the mass loss from ocean forcing (Seroussi et al., in review). In terms of ice sheet physics, these simulations do not include hydrofracture or calving-front retreat, and thus are missing positive feedbacks associated with reduced buttressing of grounded ice by ice shelves (Sun et al., in review). Solid-Earth and sea-level feedbacks are missing; these have been found to delay (but not prevent) longterm ice retreat in the Thwaites basin . At 4-km grid resolution, processes such as grounding-line retreat 20 may be under-resolved. Moreover, the depth-based melt parameterizations ignore important process such as eddy heat transfer onto the continental shelf (Stewart and Thompson, 2015) and topographic steering. Without these processes, thermal forcing in cavities could be missing critical spatial structure. Finally, there is no interaction between the ice sheet and the cavity circulation or open ocean, so that the freshwater fluxes from increased melting are unable to modify the thermal forcing.
These uncertainties suggest several lines of research to further improve ice-sheet and sea-level projections:

25
-Running similar long-term simulations with multiple ice sheet models, to quantify structural ISM uncertainties.
-Extending more CMIP AOGCM simulations beyond 2100, to generate forcing for multi-century ice sheet simulations.
-Forcing ISMs with a greater variety of scenarios, including overshoot scenarios to study whether WAIS retreat, once begun, could be stopped or reversed.
-Adding more realistic physical processes and feedbacks, such as hydrofracture and solid-Earth/sea-level effects. Some

30
ISMs already simulate these processes, and a subset of the ISMIP6 Antartic projections included hydrofracture.
-Continuing to develop comprehensive data sets of sub-ice-shelf melt rates and Southern Ocean temperature and salinity, to better force models and validate parameterizations.
the Weddell Sea and Amundsen Sea sectors during the second phase of the Marine Ice Sheet-Ocean Model Intercomparison Project (MISOMIP2).
-Developing simpler models of ice-ocean cavities and sub-shelf melting, which can emulate the results of high-resolution models but are efficient enough to run for long timescales in global models.

5
-Incorporating interactive ice sheet-ocean coupling in the next generation of Earth system models.
Many of these efforts are already under way and could contribute to future intercomparison projects, including the anticipated ISMIP7.
Near the grounding line, GLPs typically compute the basal friction as an area-weighted average of the (possibly large) friction beneath grounded ice and zero friction beneath floating ice.

15
The GLP in CISM's Glissade dynamical core is similar to the PA_GB1 scheme of Gladstone et al. (2010), extended to two dimensions. Ice thickness H and bed elevation b (with b < 0 below sea level) are located at cell centers. Ice is deemed to be floating if it satisfies a flotation condition: where ρ i is the density of ice and ρ w is the density of seawater. Based on this condition, we define a flotation function which 20 depends linearly on H and b: which gives f float = 0 at the grounding line, f float < 0 for grounded ice, and f float > 0 for floating ice. When positive, f float is equal to the depth of the ocean cavity between the ice-shelf base and the seafloor.
Given f float at cell centers, Glissade uses bilinear interpolation to compute a grounded ice fraction φ g , which is located at We compute φ g as the fraction of the bounding box for which f float < 0. We write f float as a bilinear function of x and y: f float (x, y) = a + bx + cy + dxy, where x and y are scaled to vary between 0 and 1. With the southwest corner at (0, 0) and the northeast corner at (1, 1), the 10 coefficients are where we have dropped the subscript "float".
To give an example, suppose the southwest cell is grounded and the other three are floating. The grounding line then passes through the west and south edges of the bounding box. Along the south edge we have y = 0 and therefore f = a + bx. Setting The integrals for other configurations of floating and grounded cells can also be found analytically.
The grounded fraction φ g at each vertex enters the ice dynamics via the basal friction coefficient β, defined as the ratio of basal shear stress to velocity: