ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century

Abstract

Ice flow models of the Antarctic ice sheet are commonly used to simulate its future evolution in response to different climate scenarios and assess the mass loss that would contribute to future sea level rise. However, there is currently no consensus on estimates of the future mass balance of the ice sheet, primarily because of differences in the representation of physical processes, forcings employed and initial states of ice sheet models. This study presents results from ice flow model simulations from 13 international groups focusing on the evolution of the Antarctic ice sheet during the period 2015–2100 as part of the Ice Sheet Model Intercomparison for CMIP6 (ISMIP6). They are forced with outputs from a subset of models from the Coupled Model Intercomparison Project Phase 5 (CMIP5), representative of the spread in climate model results. Simulations of the Antarctic ice sheet contribution to sea level rise in response to increased warming during this period varies between $-\mathrm{7.8}$ and 30.0 cm of sea level equivalent (SLE) under Representative Concentration Pathway (RCP) 8.5 scenario forcing. These numbers are relative to a control experiment with constant climate conditions and should therefore be added to the mass loss contribution under climate conditions similar to present-day conditions over the same period. The simulated evolution of the West Antarctic ice sheet varies widely among models, with an overall mass loss, up to 18.0 cm SLE, in response to changes in oceanic conditions. East Antarctica mass change varies between $-\mathrm{6.1}$ and 8.3 cm SLE in the simulations, with a significant increase in surface mass balance outweighing the increased ice discharge under most RCP 8.5 scenario forcings. The inclusion of ice shelf collapse, here assumed to be caused by large amounts of liquid water ponding at the surface of ice shelves, yields an additional simulated mass loss of 28 mm compared to simulations without ice shelf collapse. The largest sources of uncertainty come from the climate forcing, the ocean-induced melt rates, the calibration of these melt rates based on oceanic conditions taken outside of ice shelf cavities and the ice sheet dynamic response to these oceanic changes. Results under RCP 2.6 scenario based on two CMIP5 climate models show an additional mass loss of 0 and 3 cm of SLE on average compared to simulations done under present-day conditions for the two CMIP5 forcings used and display limited mass gain in East Antarctica.

1 Introduction

Remote sensing observations of the Antarctic ice sheet have shown continuous ice mass loss over at least the past 4 decades , in response to changes in oceanic and atmospheric conditions. This overall mass loss has large spatial variations, as regions around Antarctica experience varying climate change patterns, and individual glaciers respond differently to similar forcings depending on their local geometry and internal dynamics . To date, the Amundsen and Bellingshausen sea sectors of West Antarctica and the Antarctic Peninsula have experienced significant mass loss, while East Antarctica has had a limited response to climate change .

Despite the rapid increase in the number of observations e.g., and the recent progresses of numerical ice flow models in capturing physical processes (e.g., grounding line migration, ice front evolution) and developing assimilation methods over the past decade , the uncertainty in the Antarctic ice sheet contribution to sea level over the coming centuries remains high . Understanding processes that caused past ice sheet changes and reproducing them is critical in order to improve and gain confidence in projections of ice sheet evolution over the next decades and centuries in response to climate change. Previous modeling studies showed variable Antarctic contribution to sea level rise over the coming century, depending on the physical processes included e.g.,, model initial states e.g.,, forcing used e.g., or model parameterizations e.g.,, leading to results varying between a few millimeters to more than a meter of sea level contribution by the end of the century . Model intercomparison efforts such as Ice2Sea and SeaRISE Sea-level Response to Ice Sheet Evolution, highlighted the large discrepancies in numerical ice flow model results, even when similar climate conditions are applied for model forcing. Furthermore, most of these experiments were carried out under extremely simplified climate forcings, limiting our understanding of how ice sheets may respond to realistic climate scenarios.

ISMIP6 Ice Sheet Model Intercomparison Project for CMIP6, is the primary effort of CMIP6 (Climate Model Intercomparison Project Phase 6) focusing on ice sheets and was designed to address these questions and improve our understanding of ice sheet–climate interactions. In a first stage, ice sheet model initialization experiments initMIP, focused on the role of initial conditions and model parameters in ice flow simulations. Antarctic experiments were based on simplified forcings: the surface mass balance (SMB) was averaged between several global and regional climate models and the ocean-induced basal melt was doubled compared to the amount of basal melt estimated from remote sensing observations . These experiments were used to assess the response of ice flow models to anomalies in these external forcings . Results showed that models respond similarly to changes in SMB, while changes in ocean-induced basal melt cause a large spread in model response. The initial ice shelf extent, which varies by a factor of 2.5 between the models with the smallest and largest ice shelf extents, as well as the treatment of sub-ice-shelf basal melt close to the grounding line and the model spatial resolution, were identified as the main sources of differences between the simulations .

In this study, we focus on projections of the Antarctic ice sheet forced by outputs from CMIP5 Atmosphere–Ocean General Circulation Models (AOGCMs), including both Climate Models and Earth System Models, under different climate conditions, as CMIP6 results were not available when the experimental protocol was designed . The ensemble of simulations focuses mostly on the 2015–2100 period and is based on 21 sets of ice flow simulations submitted by 13 international institutions. We investigate the relative role of climate forcings, Representative Concentration Pathway (RCP) scenarios, ocean-induced melt parameterizations and simulated physical processes on the Antarctic ice sheet contribution to sea level and the associated uncertainties. Most of the results are presented relative to simulations with a constant climate and therefore show the impact of climate warming relative to a scenario with a constant climate. We first describe the experiment setup and the forcings used for the simulations in Sect. 2. We then detail the ice flow models that took part in this intercomparison and summarize their main characteristics in Sect. 3. Section 4 analyzes the results and assesses the impact of the different scenarios and processes explored. Finally, we discuss the results, differences between models, most vulnerable regions and the main sources of uncertainties in Sect. 5.

2 Climate forcings and experiments

ISMIP6 is an endorsed MIP (Model Intercomparison Project) of CMIP6, and experiments performed as part of ISMIP6 projections are therefore based on outputs from AOGCMs taking part in CMIP. As results from CMIP6 were not available at the time the experimental protocol was determined , it was decided to rely primarily on available CMIP5 outputs to assess the future evolution of the Greenland and Antarctic ice sheets. This choice allowed an in-depth analysis of CMIP5 AOGCM outputs and the selection of a subset of CMIP5 models that would capture the spread of climate evolution. The choice of using only a subset of AOGCMs limits the number of simulations required from each ice sheet modeling group, while still sampling the uncertainty in future ice sheet evolution associated with variations in climate models . Additional simulations based on CMIP6 are ongoing and will be the subject of a forthcoming publication.

In this section, we summarize the experimental protocol for ISMIP6-Antarctica projections, including the choice of CMIP5 climate and Earth system models, the processing of their outputs in order to derive atmospheric and oceanic forcings applicable to ice sheet models, and the processes included in the experiments. We then list the experiments analyzed in the present work. More details on the experimental protocol can be found in , while the selection of the CMIP5 model ensemble is explained in . A detailed description of the ocean melt parameterization and calibration is available in .

2.1 Selection of CMIP5 climate models

The forcings applied to ISMIP6-Antarctica projections are derived from both RCP 8.5 and RCP 2.6 scenarios, with most experiments based on RCP 8.5, in order to estimate the full extent of changes possible by 2100 with varying climate forcings. A few RCP 2.6 scenarios are used to assess the response of the ice sheet to more moderate climate changes.

After selecting CMIP5 climate and Earth system models that performed both RCP 8.5 and RCP 2.6 scenarios, they were first assessed on their ability to represent present climate conditions around the Antarctic ice sheet. A historical bias metric was computed, incorporating atmosphere and surface oceanic conditions south of 40${}^{\circ }$ S and oceanic conditions in six ocean sectors shallower than 1500 m around Antarctica. Atmospheric and surface metrics were evaluated against the European Centre for Medium-Range Weather Forecasts “Interim” reanalysis ERA-Interim,. Ocean metrics were compared to a reference climatology combining the 2018 World Ocean Atlas , EN4 ocean climatology and temperature profiles from Logger-equipped seals . Following this assessment of AOGCMs, we analyzed the changes projected between 1980–2000 and 2080–2100 in oceanic and atmospheric conditions under the RCP 8.5 scenario. We chose six CMIP5 models that performed better than the median at capturing present-day conditions and represented a large diversity in projected changes. These climate and Earth system models are CCSM4, MIROC-ESM-CHEM, and NorESM1-M for the core experiments and CSIRO-Mk3-6-0, HadGEM2-ES, and IPSL-CM5A-M for the CMIP5 Tier 2 experiments (see Sect. ). Two of these models, NorESM1-M and IPSL-CM5A-M, were also chosen to provide forcings for the RCP 2.6 scenario. We refer to for a detailed description of the model evaluation and selection.

This choice of CMIP5 models was designed both to select models that best capture the variables relevant to ice sheet evolution and to maximize the diversity in projected 21st century climate evolution, while limiting the number of simulations. CMIP5 model choices were made independently for Greenland and Antarctica, to focus on the specificities of each ice sheet and region. We derived external forcings for the Antarctic ice sheet from these CMIP5 model outputs and provided yearly forcing anomalies for participating models.

2.2 Atmospheric forcing

Using the CMIP5 models selected, atmospheric forcings were derived in the form of yearly averaged surface mass balance anomalies and surface temperature anomalies compared to the 1980–2000 period. The SMB anomalies include changes in precipitation, evaporation, sublimation and runoff and are presented in the form of water-equivalent quantities. These anomalies are then added to reference surface mass balance and surface temperature fields that are used as a baseline in the ice flow models, similar to the approach used in .

SMB conditions are often estimated using Regional Climate Models (RCMs), such as the Regional Atmospheric Climate Model RACMO, and Modèle Atmosphérique Régional MAR,, forced at their boundaries with AOGCMs outputs. As high-resolution RCM integrations for the full Antarctic ice sheet are complex and typically require additional boundary forcing and considerable time and computational resources, it was decided not to follow this approach for ISMIP6-Antarctica Projections but to use AOGCM outputs directly. Further details on the derivation of atmospheric forcing can be found in .

2.3 Oceanic forcing

Melting at the base of ice shelves is caused by the underlying circulation of ocean waters, with warmer waters and stronger currents increasing the amount of basal melt. However, converting ocean properties into basal melt forcing under the ice shelves remains challenging . Similar to what is done for the atmospheric forcing, the ocean forcing is derived from the CMIP5 AOGCMs outputs. However, the CMIP5 models do not always resolve the Antarctic continental shelf, and none include ice shelf cavities. The first task to prepare the ocean forcing was therefore to extrapolate relevant oceanic conditions (temperature and salinity) to areas not included in CMIP5 ocean models, including areas currently covered by ice that could become ice-free in the future. These areas include sub-ice-shelf cavities and areas beneath the grounded ice sheet that could be exposed to the ocean following ice thinning and grounding line retreat. Three-dimensional fields of ocean salinity, temperature and thermal forcing were then computed as annual mean values over the 1995–2100 period. We refer to for more details on the extrapolation of oceanic fields and computation of ocean thermal forcing.

Converting ocean conditions into ocean-induced melt at the base of ice shelves is an active area of research, and several parameterizations with different levels of complexity have recently been proposed for converting ocean conditions into ice shelf melt rates e.g.,. As only a limited number of direct observations of ocean conditions e.g., and ice shelf melt rates e.g., exist, these parameterizations are difficult to calibrate and evaluate. Some parameterizations are relatively complex and based on nonlocal quantities and can therefore be difficult to implement in continental-scale parallel ice sheet models. Furthermore, such parameterizations do not account for feedbacks between the ice and ocean dynamics, which are likely only captured by coupled ice–ocean models .

For these reasons, ISMIP6-Antarctica Projections includes two options that can be adopted for the sub-ice-shelf melt parameterization: (1) a standard parameterization based on a prescribed relation between ocean thermal forcing and ice shelf melting rates and (2) an open parameterization left to the discretion of the ice sheet modeling groups. Such a framework allows us to evaluate the response to a wide spectrum of melt parameterizations with the open framework while also capturing the uncertainty related to the ice sheet response under a more constrained setup in the standard framework. The standard parameterization was chosen as a trade-off between a simple parameterization that most modeling groups could implement in a limited time while capturing melt rate patterns as realistically as possible. Results from an idealized case comparing coupled ice–ocean models with different melt parameterizations suggested that a nonlocal, quadratic melt parameterization was best able to mimic the coupled ice–ocean results over a broad range of ocean forcing . These results were performed on an idealized case similar to the Marine Ice Sheet Ocean Model Intercomparison Project MISOMIP, and have not yet been tested on realistic geometries. The non-quadratic melt parameterization suggested in is as follows:

1$$\begin{array}{rl}m\left(x,y\right)& ={\mathit{\gamma }}_{\mathrm{0}}\times {\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{sw}}{c}_{\mathrm{pw}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}}}\right)}^{\mathrm{2}}\times \left(\mathrm{TF}\left(x,y,z_{\text{draft}}\right)\right.\\ & \left.+\mathit{\delta }{T}_{\text{sector}}\right)\times \left|⟨\mathrm{TF}{⟩}_{\text{draft}\in \text{sector}}+\mathit{\delta }{T}_{\text{sector}}\right|,\end{array}$$ where ${\mathit{\gamma }}_{\mathrm{0}}$ is a coefficient similar to an exchange velocity, ${\mathit{\rho }}_{\mathrm{sw}}$ the ocean density, $c_{\mathrm{pw}}$ the specific heat of sea water, ${\mathit{\rho }}_{\mathrm{i}}$ the ice density, $L_{\mathrm{f}}$ the ice latent heat of fusion, TF $(x,y,z_{\text{draft}})$ the local ocean thermal forcing at the ice shelf base, $\left|⟨\mathrm{TF}{⟩}_{\text{draft}\in \text{sector}}\right|$ the ocean thermal forcing averaged over a sector and $\mathit{\delta }{T}_{\text{sector}}$ the temperature correction for each sector. The values for ${\mathit{\gamma }}_{\mathrm{0}}$ and $\mathit{\delta }{T}_{\text{sector}}$ in this equation were calibrated combining observations of ocean conditions and remote sensing estimates of melt rates . Two calibrations based either on circum-Antarctic observations (the “MeanAnt” method) or on observations close to the grounding line of Pine Island Glacier (the “PIGL” method) were performed in a two-step process. The coefficient ${\mathit{\gamma }}_{\mathrm{0}}$ is first calibrated assuming $\mathit{\delta }T$ equal to zero and using 10${}^{\mathrm{5}}$ random samplings of melt rate and ocean temperature, so that the total melt produced under the ice shelves is similar to melt rates estimated in and . This process provides a distribution of possible ${\mathit{\gamma }}_{\mathrm{0}}$ values. The $\mathit{\delta }{T}_{\text{sector}}$ values are then calibrated for each of the 16 sectors of Antarctica seefor details, so that the melt in each basin agrees with average estimated melt in this sector. The median value of ${\mathit{\gamma }}_{\mathrm{0}}$ is used for all but two runs. These two experiments assess the impact of uncertainty in ${\mathit{\gamma }}_{\mathrm{0}}$ by using the 5th and 95th percentile values from the distribution. The second calibration, “PIGL”, uses the same process but is constrained with only a subset of observations under Pine Island ice shelf and close to its grounding line, since these values are the most relevant for highly dynamic ice streams that have the highest sub-shelf melt . This calibration leads to higher values of ${\mathit{\gamma }}_{\mathrm{0}}$, corresponding to a greater sensitivity of melt rates to changes in ocean temperature.

The choice of melt parameterization and its calibration with observations is described in detail in . For models that could not implement such a nonlocal parameterization, a local quadratic parameterization similar to Eq. (), with the nonlocal thermal forcing replaced by local thermal forcing, was also designed and calibrated to provide similar results .

2.4 Ice shelf collapse forcing

Several ice shelves in the Antarctic Peninsula have collapsed over the past 3 decades . One mechanism proposed to explain the collapse of these ice shelves is the presence of significant amounts of liquid water on their surface, which causes hydrofracturing and ultimately leads to their collapse . Other mechanisms, such as ocean surface waves, rheological weakening, surface load shifts due to water movement or basal melting , have also been proposed to explain these ice shelf collapse but are not investigated in this study. Ice shelf collapse reduces the buttressing forces provided to the upstream grounded ice and leads to acceleration and increased mass loss of the glaciers feeding them , but more dramatic consequences have been envisioned if ice shelves were to collapse in front of thick glaciers resting on retrograde bed slopes . As the presence of liquid water at the surface of Antarctic ice shelves is expected to increase in a warming climate , we propose experiments that include ice shelf collapse. The response of grounded ice streams to such a collapse is not imposed but arises from the various model representations of boundary conditions and transitions from grounded to floating ice. Apart from these experiments testing the impact of ice shelf collapse, the other experiments should not include ice shelf collapse.

Ice shelf collapse forcing is described as a yearly mask that defines the regions and times of collapse. The criteria for ice shelf collapse are based on the presence of mean annual surface melting above 725 mm over a decade, similar to numbers proposed in , and corresponding to the average melt simulated by RACMO2 over the Larsen A and B ice shelves in the decade before their collapse. The amount of surface melting was computed from CMIP5 modeled surface air temperature using the methodology described in .

2.5 List of experiments

The list of experiments for ISMIP6-Antarctica Projections is described and detailed in . It includes a historical experiment (historical), control runs (ctrl and ctrl_proj), simple anomaly experiments similar to initMIP-Antarctica (asmb and abmb), 13 core (Tier 1) experiments, and 8 Tier 2 experiments based on CMIP5 forcing. The list is repeated in Table  for completeness. In summary, these experiments include the following variations:

• 12 experiments based on RCP 8.5 scenarios from 6 CMIP5 models (open and standard melt parameterizations);

• 4 experiments based on RCP 2.6 scenarios from 2 CMIP5 models (open and standard melt parameterizations);

• 2 experiments including ice shelf collapse (open and standard melt parameterizations);

• 2 experiments testing the uncertainty in the melt parameterization (standard melt parameterization only);

• 2 experiments testing the uncertainty in the melt calibration (standard melt parameterizations only).

All experiments start in 2015, except for the historical, ctrl, asmb and abmb experiments, which start at the model initialization time. The historical experiment runs from the initialization time until the beginning of 2015, while the ctrl, asmb and abmb experiments run for either 100 years or until 2100, whichever is longer. All the other experiments run from January 2015 to the end of 2100. The ctrl_proj run is a control run similar to ctrl: a simulation under constant climate conditions representative of the recent past. The only difference is that ctrl_proj starts in 2015 and lasts until 2100, while ctrl starts from the ice models' initial state (which varies between 1850 and 2015 for the various models) and lasts at least 100 years.

Table 1

List of ISMIP6-Antarctica projections for the core (Tier 1) and Tier 2 experiments based on CMIP5 AOGCMs.

			Ocean	Ocean	Ice shelf
Experiment	AOGCM	Scenario	forcing	sensitivity	fracture	Tier
historical	None	None	Free	Medium	No	Tier 1 (Core)
ctrl	None	None	Free	Medium	No	Tier 1 (Core)
ctrl_proj	None	None	Free	Medium	No	Tier 1 (Core)
asmb	None	None	Same as ctrl $+$	Medium	No	Tier 1 (Core)
			SMB anomaly
abmb	None	None	Same as ctrl $+$	Medium	No	Tier 1 (Core)
			melt anomaly
exp01	NorESM1-M	RCP8.5	Open	Medium	No	Tier 1 (Core)
exp02	MIROC-ESM-CHEM	RCP8.5	Open	Medium	No	Tier 1 (Core)
exp03	NorESM1-M	RCP2.6	Open	Medium	No	Tier 1 (Core)
exp04	CCSM4	RCP8.5	Open	Medium	No	Tier 1 (Core)
exp05	NorESM1-M	RCP8.5	Standard	Medium${}^{\ast }$	No	Tier 1 (Core)
exp06	MIROC-ESM-CHEM	RCP8.5	Standard	Medium${}^{\ast }$	No	Tier 1 (Core)
exp07	NorESM1-M	RCP2.6	Standard	Medium${}^{\ast }$	No	Tier 1 (Core)
exp08	CCSM4	RCP8.5	Standard	Medium${}^{\ast }$	No	Tier 1 (Core)
exp09	NorESM1-M	RCP8.5	Standard	High${}^{\ast }$	No	Tier 1 (Core)
exp10	NorESM1-M	RCP8.5	Standard	Low${}^{\ast }$	No	Tier 1 (Core)
exp11	CCSM4	RCP8.5	Open	Medium	Yes	Tier 1 (Core)
exp12	CCSM4	RCP8.5	Standard	Medium${}^{\ast }$	Yes	Tier 1 (Core)
exp13	NorESM1-M	RCP8.5	Standard	PIGL	No	Tier 1 (Core)
expA1	HadGEM2-ES	RCP8.5	Open	Medium	No	Tier 2
expA2	CSIRO-MK3	RCP8.5	Open	Medium	No	Tier 2
expA3	IPSL-CM5A-MR	RCP8.5	Open	Medium	No	Tier 2
expA4	IPSL-CM5A-MR	RCP2.6	Open	Medium	No	Tier 2
expA5	HadGEM2-ES	RCP8.5	Standard	Medium${}^{\ast }$	No	Tier 2
expA6	CSIRO-MK3	RCP8.5	Standard	Medium${}^{\ast }$	No	Tier 2
expA7	IPSL-CM5A-MR	RCP8.5	Standard	Medium${}^{\ast }$	No	Tier 2
expA8	IPSL-CM5A-MR	RCP2.6	Standard	Medium${}^{\ast }$	No	Tier 2

${}^{\ast }$ For the “standard” parameterization, the low, medium and high ocean sensitivity correspond to the 5th, 50th and 95th percentile values of the “MeantAnt” ${\mathit{\gamma }}_{\mathrm{0}}$ distribution .

Most analyses presented in this study follow an “experiment minus ctrl_proj” approach, so the results provide the impact of change in climatic conditions relative to ice sheets forced with present-day conditions until 2100. We know that ice sheets respond nonlinearly to changes in climate conditions, but such an approach is necessary as ice flow model simulations often do not accurately capture the trends observed over the recent past .

3 Ice flow models

3.1 Model setups

Similar to the philosophy adopted for initMIP-Antarctica, there are no constraints on the method or datasets used to initialize ice sheet models. The exact initialization date is also left to the discretion of individual modeling groups, thus the historical experiment length varies among groups (some groups start directly at the beginning of 2015 and therefore did not submit a historical run). The resulting ensemble includes a variety of model resolutions, stress balance approximations and initialization methods, representative of the diversity of the ice sheet modeling community (see Sect.  for more details on participating models).

The only constraints imposed on the ice sheet models are that (1) models have to simulate ice shelves and the evolution of grounding lines and that (2) models have to use the atmospheric and oceanic forcings varying in time and based on CMIP5 model outputs provided. The inclusion of ice cliff failure, on the other hand, was not allowed, except in the ice shelf collapse experiments. Groups were invited to submit one or several sets of experiments, and modelers were asked to submit the full suite of open (with the melt parameterization of their choice; see Table ) and/or standard core experiments if possible. Unlike what was imposed for initMIP-Antarctica, models were free to include additional processes not specified here (e.g., changes in bedrock topography in response to changes in ice load, feedback between SMB and surface elevation).

Annual values for both scalar and two-dimensional outputs were reported on standard grids with resolutions of 4, 8, 16 or 32 km. Scalar quantities were recomputed from the two-dimensional fields submitted for consistency and in order to create regional scalars used for the regional analysis. The two-dimensional fields were also conservatively regridded onto the standard 8 km grid to facilitate spatial comparison and analysis. The outputs requested are listed in Appendix . Each group also submitted a README file summarizing the model characteristics.

3.2 Participating models

A total of 16 sets of simulations from 13 groups were submitted to ISMIP6-Antarctica projections. The groups and ice sheet modelers who ran the simulations are listed in Table . Simulations are performed using various ice flow models, a range of grid resolutions, different approximations of the stress balance equation, varying basal sliding laws, and multiple external forcings; a diverse set of processes were included in the simulations. Table  summarizes the main characteristics of the 16 sets of simulations. Short descriptions of the initialization method and main model characteristics are also provided in Appendix .

Table 2

List of participants, modeling groups and ice flow models in ISMIP6-Antarctica projections.

Contributors	Group ID	Ice flow model	Group
Thomas Kleiner,	AWI	PISM	Alfred Wegener Institute for Polar and Marine Research,
Angelika Humbert			Bremerhaven, Germany
Matthew Hoffman,	DOE	MALI	Los Alamos National Laboratory, Los Alamos, NM, USA
Tong Zhang,
Stephen Price
Ralf Greve,	ILTS_PIK	SICOPOLIS	Institute of Low Temperature Science,
			Hokkaido University, Sapporo, Japan
Reinhard Calov			Potsdam Institute for Climate Impact Research, Potsdam, Germany
Heiko Goelzer,	IMAU	IMAUICE	Institute for Marine and Atmospheric research,
Roderik van de Wal			Utrecht, The Netherlands
Nicole-Jeanne Schlegel,	JPL	ISSM	Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA
Hélène Seroussi
Christophe Dumas,	LSCE	Grisli	Laboratoire des Sciences du Climat et de l'Environnement
Aurelien Quiquet			Université Paris-Saclay, Gif-sur-Yvette, France
Gunter Leguy,	NCAR	CISM	National Center for Atmospheric Research, Boulder, CO, USA
William Lipscomb
Ronja Reese,	PIK	PISM	Potsdam Institute for Climate Impact Research, Potsdam, Germany
Torsten Albrecht,
Ricarda Winkelmann
Tyler Pelle,	UCIJPL	ISSM	University of California, Irvine, CA, USA
Mathieu Morlighem,
Hélène Seroussi			Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA
Frank Pattyn,	ULB	f.ETISh	Université libre de Bruxelles, Brussels, Belgium
Sainan Sun
Chen Zhao,	UTAS	Elmer/Ice	Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, Australia
Rupert Gladstone,			Arctic Centre, University of Lapland, Rovaniemi, Finland
Thomas Zwinger			CSC-IT Center for Science, Espoo, Finland
Jonas Van Breedam,	VUB	AISMPALEO	Vrije Universiteit Brussel, Brussels, Belgium
Philippe Huybrechts
Nicholas Golledge,	VUW	PISM	Antarctic Research Centre, Victoria University of Wellington,
Daniel Lowry			and GNS Science, Wellington, New Zealand

Table 3

List of ISMIP6-Antarctica projection simulations and main model characteristics. Numerics are defined as follows: finite difference (FD), finite elements (FE) and finite volumes (FV). Initialization methods used are as follows: spin-up (SP), spin-up with ice thickness target values SP+; see, data assimilation (DA), data assimilation with relaxation (DA+), data assimilation of ice geometry only (DA${}^{\ast }$) and equilibrium state (Eq). Melt in partially floating cells is listed as follows: melt either applied or not applied over the entire cell based on a floating condition (floating condition) and melt applied based on a sub-grid scheme (sub-grid); N/A refers to models that do not have partially floating cells. Ice front migration schemes are based on strain rate StR,, retreat only (RO), fixed front (fix), minimum thickness height (MH), and divergence and accumulated damage div,. Basal melt rate parameterization in open framework are listed as follows: linear function of thermal forcing lin,, quadratic local function of thermal forcing quad,, PICO parameterization PICO,, PICOP parameterization PICOP,, plume model Plume, and nonlocal parameterization with slope dependence of the melt Nonl4ocal + Slope,. Basal melt rate parameterization in standard framework is listed as follows: local or nonlocal quadratic function of thermal forcing and local or nonlocal anomalies .

Model name	Numerics	Stress	Resolution	Init.	Initial	Melt in partially	Ice	Open melt	Standard melt
		balance	(km)	Method	Year	floating cells	front	parameterization	parameterization
AWI_PISM	FD	Hybrid	8	Eq	2005	Sub-grid	StR	Quad	Nonlocal
DOE_MALI	FE/FV	HO	2–20	DA+	2015	Floating condition	Fix	N/A	Nonlocal anom.
ILTS_PIK_SICOPOLIS	FD	Hybrid	8	SP+	1990	Floating condition	MH	N/A	Nonlocal
IMAU_IMAUICE1	FD	Hybrid	32	Eq	1978	No	Fix	N/A	Local anom.
IMAU_IMAUICE2	FD	Hybrid	32	SP	1978	No	Fix	N/A	Local anom.
JPL1_ISSM	FE	SSA	2–50	DA	2007	Sub-grid	Fix	N/A	Nonlocal
LSCE_GRISLI	FD	Hybrid	16	SP+	1995	N/A	MH	N/A	Nonlocal
NCAR_CISM	FE/FV	L1L2	4	SP+	1995	Sub-grid	RO	Nonlocal	Nonlocal
								+ Slope
PIK_PISM1	FD	Hybrid	8	SP	1850	Sub-grid	StR	PICO	N/A
PIK_PISM2	FD	Hybrid	8	SP	2015	Sub-grid	StR	PICO	N/A
UCIJPL_ISSM	FE	HO	3–50	DA	2007	Sub-grid	Fix	PICOP	Nonlocal
ULB_FETISH_16km	FD	Hybrid	16	DA${}^{\ast }$	2005	N/A	Div	Plume	Nonlocal
ULB_FETISH_32km	FD	Hybrid	32	DA${}^{\ast }$	2005	N/A	Div	Plume	Nonlocal
UTAS_ElmerIce	FE	Stokes	4–40	DA	2015	Sub-grid	Fix	N/A	Local
VUB_AISMPALEO	FD	SIA+SSA	20	SP	2000	N/A	MH	N/A	Nonlocal anom.
VUW_PISM	FD	Hybrid	16	SP	2015	No	StR	Lin	N/A

The 16 sets of submitted simulations have been performed using 10 different ice flow models. Amongst the simulations, 3 use the finite-element method, 2 use a combination of finite element and finite volume, and the remaining 11 the finite-difference method. One simulation is based on a Full-Stokes stress balance, two use the 3D higher-order approximations HO,, one is based on the L1L2 approximation and one is based on the shelfy-stream approximation SSA,, while the other simulations combine the SSA with the shallow ice approximation SIA,. The model resolutions range between 4 and 20 km for models that use regular grids but can be as low as 2 km in specific areas, such as close to the grounding line or shear margins for models with spatially variable resolution .

As in initMIP-Antarctica , the initialization procedure reflects the broad diversity in the ice sheet modeling community: two simulations start from an equilibrium state, five models start from a long spin-up and three simulations from data assimilation of recent observations. The remaining simulations combine the latter two approaches by either adding constraints to their spin-up (three simulations) or running short relaxations after performing data assimilation (three simulations). The initialization year varies between 1850 and 2015, therefore the length of the historical experiment varies between 0 and 115 years.

All submissions are required to include grounding line evolution (see Sect. ), but the treatment of grounding line evolution and ocean melt in partially floating grid cells is left to the discretion of the modeling groups. Simulating ice front evolution (i.e., calving) in the simulations is also encouraged but not required, and the choice of ice front parameterization is free. Six models use a fixed ice front (except for the ice shelf collapse experiments, for which retreat is imposed), while the other models rely on a combination of minimum ice thickness, strain rate values and stress divergence to evolve their ice front position.

The simulations were performed using the open and/or standard melt parameterizations: five sets of simulations include results based on both the open and standard framework, leading to a total of 21 sets of simulations in total when the open and standard parameterizations are analyzed separately; this parameterization affects the results significantly, therefore the open and standard parameterizations are analyzed separately from now on. Ocean-induced melt rates under ice shelves follow the standard melt framework described in Sect.  for 13 sets of simulations: 10 submissions use the nonlocal form, while 3 are based on the local form, and three of these 13 sets of simulations are based on the nonlocal or local anomaly forms . The open melt framework was used by eight sets of simulations that rely on a linear melt dependence of thermal forcing , a quadratic local melt parameterization with a calibration different than the standard framework, a plume model , a box model , a combination of box and plume models , or a nonlocal quadratic melt parameterization combined with ice shelf basal slope .

The modeling groups were asked to submit a full suite of core experiments based on the standard melt parameterization, the open one or both. Most groups were able to do so, but several groups did not submit the ice shelf collapse experiments, and one group (UTAS_ElmerIce) ran only a subset of experiments due to the high cost of running a Full-Stokes model of the Antarctic continent. Simulations that initialize their model in January 2015 (see Table ) do not have a historical run, and their ctrl and ctrl_proj are therefore identical. Seven submissions also performed some or all of the Tier 2 experiments (expA1–A8). Table  lists all the experiments done by the modeling groups.

Table 4

List of experiments performed as part of ISMIP6-Antarctica projections by the modeling groups.

Experiment	AWI_PISM	DOE_MALI	ILTS_PIK_SICOPOLIS	IMAU_IMAUICE1	IMAU_IMAUICE2	JPL1_ISSM	LSCE_GRISLI	NCAR_CESM	PIK_PISM1	PIK_PISM2	UCIJPL_ISSM	ULB_fETISh_16	ULB_fETISh_32	UTAS_ElmerIce	VUB_AISMPALEO	VUW_PISM
historical	X		X	X	X	X	X	X	X		X	X	X		X	X
ctrl	X	X	X	X	X	X	X	X	X	X	X	X	X	X	X	X
ctrl_proj	X	X${}^{\ast }$	X	X	X	X	X	X	X	X${}^{\ast }$	X	X	X	X${}^{\ast }$	X	X
asmb	X	X	X	X	X	X	X	X	X	X	X	X	X		X	X
abmb	X	X	X	X	X	X	X	X	X	X	X	X	X	X	X	X
exp01	X							X	X	X	X	X	X			X
exp02	X							X	X	X	X	X	X			X
exp03	X							X	X	X	X	X	X			X
exp04	X							X	X	X	X	X	X			X
exp05	X	X	X	X	X	X	X	X			X	X	X	X	X
exp06	X	X	X	X	X	X	X	X			X	X	X	X	X
exp07	X	X	X	X	X	X	X	X			X	X	X		X
exp08	X	X	X	X	X	X	X	X			X	X	X		X
exp09	X	X	X	X	X	X	X	X			X	X	X		X
exp10	X	X	X	X	X	X	X	X			X	X	X		X
exp11	X										X	X	X
exp12	X	X	X	X	X	X	X				X	X	X
exp13	X	X	X	X	X	X	X	X			X	X	X	X	X
expA1	X							X				X	X
expA2	X							X				X	X
expA3	X							X				X	X
expA4	X							X				X	X
expA5	X		X		X	X	X	X			X	X	X		X
expA6	X		X		X	X	X	X			X	X	X		X
expA7	X		X		X	X	X	X			X	X	X		X
expA8	X		X		X	X	X	X			X	X	X

${}^{\ast }$ Indicates simulations initialized directly at the beginning of 2015 for which ctrl and ctrl_proj experiments are identical.

4 Results

We detail the simulation results here. We start by describing the initial state and the historical and control runs. We then analyze the NorESM1-M RCP 8.5 runs, and the RCP 8.5 simulations based on the six different CMIP5 model forcings. Next, we compare the RCP 8.5 and RCP 2.6 results for the two CMIP5 models selected to provide RCP 2.6 scenario forcings. We then investigate the effect of uncertainty in the melt parameterization and calibration. Finally, we explore the role of ice shelf collapse.

Results based on the open and standard melt parameterizations are combined, except in Sect. , where we investigate the difference between these approaches. This means that 21 independent sets of results are extracted from the 16 submissions (8 based on the open melt framework and 13 based on the standard framework). No weighting based on the number of submissions or agreement with observations is applied.

4.1 Historical run and 2015 conditions

As the initialization date for different models varies, all models run a short historical simulation until 2015. The length of this simulation varies between 165 years for PIK_PISM1, which starts in 1850, and 0 years for the three models (DOE_MALI, PIK_PISM2 and UTAS_ElmerIce) that start directly in 2015. During the historical run, simulations are forced with oceanic and atmospheric conditions representative of the conditions estimated during this period. The total annual SMB over Antarctica varies between 2140 and 3230 Gt yr${}^{-\mathrm{1}}$, with large interannual variations of up to 600 Gt yr${}^{-\mathrm{1}}$ (see Fig. a). The total annual ocean-induced basal melt rates under ice shelves during the historical period varies between 0 and 4200 Gt yr${}^{-\mathrm{1}}$, with large interannual variations up to 500 Gt yr${}^{-\mathrm{1}}$. The ice volume above floatation, however, experiences limited variations during the historical period, up to a 6000 Gt change (Fig. b).

Figure 1

Evolution of surface mass balance (a, in Gt yr${}^{-\mathrm{1}}$), basal melt rate (b, in Gt yr${}^{-\mathrm{1}}$) and volume above floatation (c, in Gt) during the historical and ctrl_proj experiments for all the simulations performed with the open and standard framework. Note the different scale on the time axis prior to 1950.

[Figure omitted. See PDF]

All historical simulations end in December 2014, at which point the projection experiments start. Figure  shows the total ice and floating ice extent for all submissions at the beginning of the experiments. The simulated ice-covered area varies between 1.36 and $\mathrm{1.45}\times {\mathrm{10}}^{\mathrm{7}}$ km${}^{\mathrm{2}}$, or 6.0 %. There is good agreement between the modeled ice extent and the observed ice front around the entire continent and a smaller spread compared to the initMIP-Antarctica submissions, in which the ice extent varied between 1.35 and $\mathrm{1.50}\times {\mathrm{10}}^{\mathrm{7}}$. The extent of ice shelves shown in Fig. b varies between 1.19 and $\mathrm{1.92}\times {\mathrm{10}}^{\mathrm{6}}$ km${}^{\mathrm{2}}$, which is a much smaller spread in the results than in the initMIP-Antarctica experiments (between 0.92 and $\mathrm{2.51}\times {\mathrm{10}}^{\mathrm{6}}$) and a better agreement with observations . Not only the large ice shelves but also the smaller ice shelves of the Amundsen and Bellingshausen sea sectors, the Antarctic Peninsula, and Dronning Maud Land have a location and extent that is usually within several tens of kilometers of observations. A few models have ice shelves that extend slightly farther than the present-day ice over large parts of the continent, but they extend only a few tens of kilometers past the observed ice front location. Finally, the location of the grounding line on the Ross ice streams fluctuates by several hundred kilometers between the models, which is not surprising as the Ross ice streams rest over relatively flat bedrock, and thus small changes in model configuration lead to large variations in the grounding line position. The 2015 ice volume and ice volume above floatation are reported in Table  and on Fig. c. They indicate a variation of 6.8 % of the total ice mass among the simulations, between 2.31 and $\mathrm{2.49}\times {\mathrm{10}}^{\mathrm{7}}$ Gt, and a variation of 7.7 % in the total ice mass above floatation, between 1.99 and $\mathrm{2.15}\times {\mathrm{10}}^{\mathrm{7}}$ Gt or between 55.0 and 59.4 m of sea level equivalent (SLE), when the latest estimate is $\mathrm{57.9}\pm \mathrm{0.9}$ m . Figure  shows the root-mean-square error (RMSE) between modeled and observed thickness and velocity at the beginning of the experiments. The RMSE thickness varies between 92 and 396 m, while the RMSE velocity varies between 79 and 446 m yr${}^{-\mathrm{1}}$, which is comparable to values reported for initMIP-Antarctica .

4.2

Control experiment ctrl_proj

All the experiments start from the 2015 configuration and are run with varying atmospheric and oceanic forcings until 2100. The ctrl_proj experiment also starts from this configuration but is run with constant climate conditions (no oceanic or atmospheric anomalies added), similar to those observed over the past several decades. The exact choice of forcing conditions for this run was not imposed and therefore varies between the simulations. Figure  shows that, similarly to the historical run, the SMB and basal melt vary significantly between the simulations. The SMB varies between 2320 and 3090 Gt yr${}^{-\mathrm{1}}$, while the basal melt varies between 0 and 3740 Gt yr${}^{-\mathrm{1}}$. However, unlike what is observed in the historical run, there is limited interannual fluctuation, since a mean climatology is used for this run.

During the 86 years of the ctrl_proj experiment, the simulated evolution of ice mass above floatation varies between $-\mathrm{51}$ 500 and $\mathrm{46}$ 700 Gt (between $-\mathrm{130}$ and 142 mm SLE; see Table ). The trend in the ctrl_proj mass above floatation is significant in several models and negligible in others. As in initMIP-Antarctica, models initialized with a steady state or a spin-up tend to have smaller trends than models initialized with data assimilation. Since constant climate conditions are applied, trends cannot be considered a physical response of the Antarctic ice sheet but rather highlight the effect of model choices to initialize the simulation and represent ice sheet evolution, the lack of physical processes , the limited number or inaccuracy of observations , and the need to better integrate observations in ice flow models .

Figure 2

Total (a) and floating (b) ice extent at the beginning of the experiments (January 2015). Colors indicate the number of models simulating total ice (a) and floating ice (b) extent at every point of the 8 km grid. Black lines are observations of the total and floating ice extent, respectively .

[Figure omitted. See PDF]

Figure 3

Root-mean-square error in ice thickness (a, in m) and ice velocity (b, in m yr${}^{-\mathrm{1}}$) between modeled and observed values at the beginning of the experiments (January 2015).

[Figure omitted. See PDF]

All the results presented in the remainder of the paper are shown relative to the outputs from the ctrl_proj experiment. As a consequence, these results should be interpreted as the models' simulated response to additional climate change compared to a scenario where the climate remains constant and similar to the past few decades. Submissions that include both open and standard experiment results can have significant variations in their historical and ctrl_proj depending on whether the open or standard melt parameterization is used (see Fig.  and Tables  and ). We therefore remove the trends from the ctrl_proj open or standard melt parameterization from the experiments based on the open or standard framework, respectively.

4.3 Projections under RCP 8.5 scenario with NorESM1 forcing

The NorESM1-M RCP 8.5 scenario (exp01 and exp05; see Table ) produces mid-to-high changes in the ocean and low changes in the atmosphere over the 21st century compared to other CMIP5 AOGCMs . The effects of these changes on the simulated evolution of the Antarctic ice sheet are summarized in Figs. , and . Figure  shows that under this forcing, the ice sheet loses a volume above floatation varying between $-\mathrm{26}$ and 166 mm of SLE between 2015 and 2100, relative to ctrl_proj experiments. The impact of the forcing remains limited until 2050, with changes between $-\mathrm{2}$ and 27 mm. It quickly increases after 2050, at which point the simulations start to diverge strongly.

Figure 4

Evolution of ice volume above floatation (in mm SLE) over 2015–2100 from the NorESM1-M RCP 8.5 scenario (exp01 and exp05) relative to ctrl_proj.

[Figure omitted. See PDF]

Figure  shows that the sea level contribution and the mechanisms at play vary significantly for the West Antarctic ice sheet (WAIS), East Antarctic ice sheet (EAIS) and the Antarctic Peninsula. In the WAIS, the additional SMB is limited to a few millimeters (between $-\mathrm{2}$ and 2 mm SLE), and all models predict a mass loss varying between 0 and 154 mm SLE relative to ctrl_proj. EAIS experiences a significant increase in SMB, with a cumulative additional SMB causing between 20 and 25 mm SLE of mass gain relative to ctrl_proj. This mass gain is partially offset by the dynamic response of outlet glaciers in the EAIS, resulting in a total volume change varying between a 24 mm SLE mass gain and 38 mm SLE mass loss. The small size of the Antarctic Peninsula and limited mass of its glaciers make it a smaller contributor to sea level change compared to WAIS and EAIS: the contribution to sea level varies between $-\mathrm{6}$ and 1 mm SLE relative to ctrl_proj, with a signal split between the additional SMB (between 0 and 3 mm SLE mass gain) and dynamic response. These results therefore highlight the contrast between the EAIS and the Antarctic Peninsula, which are projected to either gain or lose mass and where SMB changes are relatively large, and the WAIS, which is dominated by a dynamic mass loss caused by the changing ocean conditions.

Figure 5

Regional change in volume above floatation (in mm SLE) and integrated SMB changes over the grounded ice (diamond shapes, in mm SLE) for the 2015–2100 period under medium RCP 8.5 forcing from NorESM1-M RCP 8.5 scenario (exp01 and exp05) relative to ctrl_proj.

[Figure omitted. See PDF]

Regions with the largest simulated changes can also be seen in Fig. , which shows the mean change in thickness and velocity between 2015 and 2100 for the 21 NorESM1-M simulations relative to ctrl_proj. Most Antarctic ice shelves thin by 20 m or more over the 86-year simulation, with the Ross ice shelf experiencing the largest thinning of about 75 m on average (Fig. a). This thinning does not propagate to the ice streams feeding the ice shelves, except for Thwaites Glacier in the Amundsen Sea sector and Totten Glacier in Wilkes Land. Many coastline regions, on the other hand, experience a small thickening, as is the case for the Antarctic Peninsula, Dronning Maud Land and Kemp Land, where the relative thickening is about 6 m next to the coast. Variations between the simulation are large and dominate the signal in many places (Fig. c). Changes in velocity (Fig. b) over ice shelves are more limited and not homogeneous, with acceleration close to the grounding line areas and slowdown close to the ice front, as observed for the Ross and Ronne-Filchner ice shelves. Some accelerations are observed on grounded parts of Thwaites, Pine Island and Totten glaciers as well. However, there is a large discrepancy in velocity changes among the simulations, and the standard deviation in velocity change is larger than the mean signal over most of the continent (Fig. d).

Figure 6

Mean (a and b) and standard deviation (c and d) of simulated thickness change (a and c, in m) and velocity change (b and d, in m yr${}^{-\mathrm{1}}$) between 2015 and 2100 under medium forcing from the NorESM1-M RCP 8.5 scenario (exp01 and exp05) relative to ctrl_proj.

[Figure omitted. See PDF]

4.4 Projections under RCP 8.5 scenario with various forcings

Outputs from six CMIP5 AOGCMs were used to perform RCP 8.5 experiments (see Table ). Figure  shows the evolution of the simulated ice volume above floatation relative to ctrl_proj for all the individual RCP 8.5 simulations performed, as well as the mean values for each AOGCM. As seen above for NorESM1-M, changes are small for most simulations until 2050, after which differences between AOGCMs and ice flow simulations start to emerge. Runs with HadGEM2-ES lead to significant sea level rise, with a mean ice mass loss of 96 mm SLE (standard deviation: 72 mm SLE) for the 15 submissions of expA1 and expA5. Runs performed with CCSM4 show the largest ice mass gain, with a mean gain of 37 mm SLE (standard deviation: 34 mm SLE) for the 21 submissions of exp04 and exp08. Results for CSIRO-MK3 and IPSL-CM5A-MR are similar to CCSM4 at a continental scale but with slightly lower mass gain on average, while results from MIROC-ESM-CHEM simulate very little change, with a mean mass loss of 3 mm SLE.

Figure 7

Evolution of ice volume above floatation (in mm SLE) over the 2015–2100 period with medium forcing from the six CMIP5 models and RCP 8.5 scenario relative to ctrl_proj. Thin lines show results from individual ice sheet model simulations, and thick lines show mean values averaged for each CMIP5 model forcing. Bars on the right show the spread of results in ice flow models and mean values for the six CMIP5 forcings in 2100.

[Figure omitted. See PDF]

Figure  shows the regional differences in these contributions relative to ctrl_proj. Simulations suggest that WAIS will lose mass on average with four of the CMIP5 model forcings and gain mass with CSIRO-MK3 and IPSL-CM5A-MR. For the EAIS, results from five out of six CMIP5 model forcings lead to a mass gain on average, while HadGEM2-ES forcing causes a mass loss in the EAIS, with $\mathrm{23}\pm \mathrm{26}$ mm SLE. Uncertainties are larger for WAIS than EAIS and larger for CMIP5 models that experience larger changes in ocean conditions. This is similar to what was observed in initMIP-Antarctica : in that study, changes in oceanic conditions (based on a forcing much simpler than is used in the current study) lead to a much larger spread in ice sheet evolution than changes in SMB. Changes in the Antarctic Peninsula lead to mass change between $-\mathrm{6}$ and 6 mm SLE on average.

Figure 8

Regional change in volume above floatation (in mm SLE) for 2015–2100 from six CMIP5 model forcings under the RCP 8.5 scenario with median forcing, relative to ctrl_proj. Black lines show standard deviations.

[Figure omitted. See PDF]

4.5 Projections under RCP 8.5 and RCP 2.6 scenarios

Two CMIP5 models were chosen to run both RCP 8.5 and RCP 2.6 experiments: NorESM1-M and IPSL-CM5A-MR. Figure  shows the evolution of the Antarctic ice sheet under these two scenarios relative to ctrl_proj for both models. Only ice flow models that performed both RCP 8.5 and RCP 2.6 experiments were used to compare these scenarios, so two RCP 8.5 runs were not included, leading to the analysis of 20 NorESM1-M and 13 IPSL-CM5A-MR pairs of experiments.

Figure 9

Impact of RCP scenario on projected evolution of ice volume above floatation for the NorESM1-M (a) and IPSL (b) models. Red and blue curves show mean evolution for RCP 8.5 and RCP 2.6, respectively, and the shaded background shows the standard deviation.

[Figure omitted. See PDF]

Results from NorESM show no significant change between the two scenarios in terms of simulated ice volume above floatation by 2100 (Fig. a). Both scenarios lead to a mean sea level contribution of about 25 mm SLE in 2100, with a higher standard deviation for the RCP 8.5 scenario (49 mm for RCP 8.5 and 37 mm for RCP 2.6). However, the overall similar behavior hides large regional differences revealed in Fig. a. The WAIS loses more mass while the EAIS gains more ice mass in RCP 8.5 compared to RCP 2.6. The additional SMB is greater for all regions under RCP 8.5 compared to RCP 2.6 (18 mm additional SLE in the EAIS and 2 mm additional SLE for the WAIS and Antarctic Peninsula) but is compensated for by a large dynamic response to ocean changes in both WAIS and EAIS.

Figure 10

Regional change in volume above floatation (in mm SLE) and integrated SMB changes over the grounded ice (diamond shapes, in mm SLE) for 2015–2100 under the RCP 8.5 (red) and RCP 2.6 (blue) scenario forcings from NorESM1-M (a) and IPSL (b) relative to ctrl_proj from individual model simulations.

[Figure omitted. See PDF]

Simulations based on IPSL-CM5A-MR forcing, on the other hand, show significant differences in ice contribution to sea level at a continental scale. Ice contributes to $-\mathrm{33}\pm \mathrm{15}$ mm SLE for the RCP 8.5 scenario and $\mathrm{1}\pm \mathrm{9}$ mm SLE for the RCP 2.6 scenario (Fig. ). For RCP 2.6, the overall mass loss in the WAIS is compensated for by mass gain in the EAIS, leading to an overall ice mass that is nearly constant (Fig. ). For RCP 8.5, there are large mass gains in all ice sheet regions as SMB increases significantly. Only a few simulations show mass loss of the WAIS relative to ctrl_proj. Similar to what is observed for NorESM1-M, the uncertainty is larger for RCP 8.5, as oceanic changes are more pronounced in this scenario.

Overall, these two CMIP5 models respond very differently to increased carbon concentrations, which is reflected in the differences in ice sheet evolution.

4.6 Impact of ice shelf basal melt parameterization

All of the RCP 8.5 experiments were simulated with the open (exp01-04) and standard (exp05-08) melt frameworks (Table ). The standard framework allows us to assess the uncertainty associated with ice flow models when the processes controlling ice–ocean interactions are fixed. The open framework, in contrast, allows for additional uncertainties due to the physics of ice–ocean interactions that remain a subject of active research . We now investigate the effects of these different approaches on simulation results.

Figure  shows the cumulative ocean-induced basal melt and the change in ice volume above floatation between 2015 and 2100 and relative to ctrl_proj for the six RCP 8.5 experiments and for the 8 and 14 submissions using the open and standard melt frameworks, respectively. The basal melt applied in the standard framework is higher than the basal melt resulting from the open framework for about half of the experiments and Antarctic regions and lower for the other half. However, despite the similar melt rates applied, the sea level contribution relative to ctrl_proj is higher (either more mass loss or less mass gain) in the open framework than in the standard framework in the WAIS and EAIS, except for IPSL in the WAIS. Numbers are small and similar in the Antarctic Peninsula. The mean additional sea level contribution (either more mass loss or less mass gain) simulated in the open framework is 25 mm SLE for WAIS and 20 mm for EAIS. The standard deviation of both basal melt and sea level contribution is larger in the open melt framework (see Fig. ), which is expected given the additional flexibility in the melt parameterization and the wide range of melt parameterizations used in the open framework (see Table ).

Figure 11

Regional change in integrated basal melt (a, in Gt) and volume above floatation (b, in mm SLE) for 2015–2100 under medium forcing from the six CMIP5 AOGCMs using RCP 8.5 forcing, relative to ctrl_proj for the open (solid patterns) and standard basal melt (dashed patterns) frameworks. Black lines show the standard deviations.

[Figure omitted. See PDF]

4.7 Impact of ice shelf melt uncertainties

The effect of uncertainties in the melt rate parameterization is assessed exclusively for the standard melt parameterization framework, for which different choices of parameters can be used in a similar way by all models (exp05, exp09, exp10 and exp13 in Table ). Here we assess the effect of two sources of uncertainty that affect the choice of ${\mathit{\gamma }}_{\mathrm{0}}$ and the regional ${\mathit{\delta }}_T$ values. The melt parameterization provides a distribution of ${\mathit{\gamma }}_{\mathrm{0}}$, and the median value is used for most experiments (see Table ). Two experiments (exp09 and exp10) use the 5th and 95th percentile values of the distribution to estimate the effect of parameter uncertainty on basal melt and ice mass loss. A third experiment investigates the effect of the dataset used to calibrate the melt parameterization (exp13): instead of using all the melt rates and ocean conditions around Antarctica, it uses only the high melt values near the Pine Island ice shelf grounding line (“PIGL” coefficient; see Sect. ), which results in ${\mathit{\gamma }}_{\mathrm{0}}$ that is an order of magnitude higher . All of these experiments are based on NorESM1-M and RCP 8.5, so the applied SMB is similar in all experiments (only the basal melt differs). The initial basal melt is calibrated to be equal to observed values in each case and for each Antarctic basin, so only the initial distribution of melt and its evolution in time vary while its total initial magnitude is similar.

Figure a shows the effect of using the 5th, 50th and 95th percentile values of the ${\mathit{\gamma }}_{\mathrm{0}}$ distribution for models that performed these three experiments. The total melt starts from similar values but diverges quickly as ocean conditions change. By 2100, the additional mean total melt applied relative to the control experiment is 3030 Gt yr${}^{-\mathrm{1}}$ for the median value, while it is 2540 and 3460 Gt yr${}^{-\mathrm{1}}$, respectively, for the 5th and 95th percentile values of the ${\mathit{\gamma }}_{\mathrm{0}}$ distribution. While these differences represent about 15 % of the additional melt applied, they fall largely within the spread of basal melt values applied for the median ${\mathit{\gamma }}_{\mathrm{0}}$ for the different simulations (caused by the different model geometries) and are smaller than interannual variations. Impacts of these changes on ice dynamics are shown in Fig. c. The mean sea level contribution with the median ${\mathit{\gamma }}_{\mathrm{0}}$ is 6.9 mm SLE, while it is $-\mathrm{0.7}$ and 12.0 mm SLE in 2100 for the 5th and 95th percentile compared to the ctrl_proj experiment. The overall evolution of Antarctica remains similar only for a couple of decades, at which point the three experiments start to diverge.

Figure 12

Impact of basal melt parameterization (a and c; 5th, 50th and 95th percentile values of ${\mathit{\gamma }}_{\mathrm{0}}$ distribution) and calibration (b and d; “MeanAnt” and “PIGL” calibrations) on basal melt evolution (a and b, in Gt yr${}^{-\mathrm{1}}$) and ice volume above floatation (c and d, in mm SLE) relative to ctrl_proj over 2015–2100. Lines show the mean values, and the shaded background shows the simulation spreads. Note that the $y$ axes differ in all plots.

[Figure omitted. See PDF]

Figure  also highlights the role of the calibration method. The “MeanAnt” and “PIGL” experiments (exp05 and exp13) start with similar total melt values and are both calibrated to be in agreement with current observations of melt (because models have initial geometries that differ from observations, they can have some differences in the amount of total initial melt). The total melt diverges between the two experiments after just a few years, and continues to diverge during the 21st century as ocean conditions and ice shelf configurations change, reaching 3030 and 5790 Gt yr${}^{-\mathrm{1}}$ on average in 2100 for the “MeanAnt” and “PIGL” experiments relative to the ctrl_proj experiment (Fig. b), respectively. The effect on ice dynamics and sea level is large, with a 12 times larger mean contribution to sea level by 2100 relative to ctrl_proj for the “PIGL” experiment, reaching a mean SLE contribution of 30 mm; see Fig. d. This is the simulation with the greatest amounts of ice loss, with models predicting mass loss of up to 30 cm SLE by 2100 compared to the ctrl_proj experiment. This melt parameterization causes larger melt rates close to grounding lines and higher sensitivity to ocean warming, as ${\mathit{\gamma }}_{\mathrm{0}}$ is an order of magnitude larger for the “PIGL” parameterization than for the “MeanAnt” parameterization. Thus, this run represents an upper end to plausible values for sub-shelf melting, yet it is calibrated to simulate initial basal melting in agreement with present-day observations. It also highlights the nonlinear ice sheet response to submarine melt forcing: the doubling of basal melt relative to the ctrl_proj experiment leads to a 10 times greater ice mass loss relative to the ctrl_proj results.

4.8 Impact of ice shelf collapse

The effect of ice shelf collapse is tested with exp11 and exp12 for the open and standard frameworks, respectively (Table ). These experiments are based on outputs from CCSM4 and are similar to exp04 and exp08: the SMB and ocean thermal forcing are similar, so the two sets of experiments only differ by the inclusion of ice shelf collapse. As mentioned in Sect. , the processes included in the response of the tributary ice streams feeding into these ice shelves is left to the judgment of modeling groups. However, no group included the marine ice cliff instability following ice shelf collapse. Only the 14 simulations (including 4 open and 10 standard melt parameterizations) that performed the ice shelf collapse experiments are included in the analysis of ice shelf collapse. Results from 7 simulations of exp04 and exp08 were therefore excluded from the ensemble with no ice shelf collapse.

As shown in , the presence of significant liquid water on the surface of ice shelves is limited to less than 60 000 km${}^{\mathrm{2}}$ until 2040, so ice shelf collapse is marginal. Starting in 2040, it rapidly increases, reaching 460 000 km${}^{\mathrm{2}}$ by 2100. The evolution of ice shelf extent in the ice sheet simulations reflects this evolution: Fig. a shows the evolution of ice shelf extent for the CCSM4 simulations with and without ice shelf collapse. As the external forcings are similar in both runs, the difference comes from the ice shelf collapse and the response to this collapse. In the simulations without collapse, ice shelf extent remains relatively constant, with 11 000 km${}^{\mathrm{2}}$ change on average compared to ctrl_proj on average. When ice shelf collapse is included, ice shelf extent is reduced by 66 000 km${}^{\mathrm{2}}$ between 2015 and 2100 compared to the ctrl_proj runs on average for the 14 ice sheet simulations.

While ice shelf collapse does not directly contribute to sea level rise, the dynamic response of the ice streams to the collapse leads to an average of 28 mm SLE difference between the two scenarios relative to the ctrl_proj experiment (Fig. a). These changes occur largely over the Antarctic Peninsula, next to George VI ice shelf, but also on Totten Glacier (see Fig. a). Including ice shelf collapse leads to a concurrent acceleration of up to 100 m yr${}^{-\mathrm{1}}$ in these same regions (see Fig. b). However, large uncertainties dominate these model responses.

Figure 13

Evolution of basal melt (a, in Gt yr${}^{-\mathrm{1}}$) and ice volume above floatation relative to ctrl_proj (b, in mm SLE) without (red) and with (cyan) ice shelf collapse over the 2015–2100 period under the CCSM4 RCP 8.5 forcing. Lines show the mean values, and the shaded background shows the standard deviations. Note the negative values of sea level contribution and therefore mass gain in panel (b).

[Figure omitted. See PDF]

Figure 14

Mean simulated thickness change (a, in m) and velocity change (b, in m yr${}^{-\mathrm{1}}$) between 2015 and 2100 with ice shelf collapse under CCSM4 RCP 8.5 scenario (exp11 and exp12) relative to similar experiments without ice shelf collapse (exp04 and exp08). Hatched areas show areas experiencing ice shelf collapse by 2100.

[Figure omitted. See PDF]

The ice shelf collapse experiments are based on CCSM4, as this model shows the largest potential for ice shelf collapse out of the six AOGCMs selected . Similar experiments performed with other AOGCMs are therefore expected to show a lower response to ice shelf collapse.

5 Discussion

ISMIP6-Antarctica projections under the RCP 8.5 scenario show a large spread of Antarctic ice sheet evolution over 2015–2100, depending on the ice flow model adopted, the CMIP5 forcings applied, the ice sheet model processes included, and the form and calibration of the basal melt parametrization. The results presented here suggest the contribution to sea level with the “MeanAnt” calibration in response to this scenario varies between a sea level drop of 7.6 cm and a sea level increase of over 27 cm, compared to a constant climate similar to that of the past few decades. Contributions up to 30 cm are also simulated when the melt parameterization is calibrated to produce high melt rates near Pine Island's grounding line (see Sect. ). The latter parameterization is calibrated with the same present-day observations but has a much stronger sensitivity to ocean forcing , leading to more rapid increases in basal melting as ocean waters in ice shelf cavities warm. As observations of ocean conditions within ice shelf cavities and resulting ice shelf melt rates remain limited, these numbers cannot be excluded from consideration.

All the simulations results reported here describe Antarctic mass loss relative to that from a constant climate, so the mass loss trend over the past few decades needs to be added to obtain a total Antarctic contribution to sea level through 2100. The recent IMBIE assessment estimated the Antarctic mass loss to be between 38 and 219 Gt yr${}^{-\mathrm{1}}$, depending on the time period considered , which corresponds to a cumulative mass loss of 9 and 52 mm over 2015–2100. Adding this to the range of Antarctic mass loss simulated as part of ISMIP6 gives a range of between $-\mathrm{6.7}$ and 35 cm SLE. These numbers cover the wide range of results previously published e.g., but do not reproduce the highest contributions up to 1 m previously reported. These numbers show less spread than the simulations performed under the SeaRISE experiments, mostly due to the lower basal melt anomalies applied under ice shelves . They are also similar to numbers presented by , where the likely range (5 %–95 % of model range) of Antarctic contribution to global mean sea level rise between the 1986–2005 period and 2100 under the RPC 8.5 scenario was between $-\mathrm{8}$ and 14 cm.

The simulated response of the ice sheet to changes in ocean forcings has significant spatial variation, suggesting that some sectors of the ice sheet are significantly more vulnerable to changes in ocean circulation than others. Figure  shows the sensitivity of the 18 Antarctic basins to changes in oceanic conditions using all the RCP 8.5 experiments performed by all the ice sheet models based on medium ocean conditions. The dynamic mass loss (total ice above floatation mass loss minus SMB change) between 2015 and 2100 is represented as a function of the cumulative ocean-induced melt over the same period, both relative to ctrl_proj. The Amundsen Sea sector and Wilkes Land show the largest dynamic response and therefore sensitivity to increase in ocean-induced basal melting. Glaciers feeding the west side of the Ross ice shelf show very small response despite relatively large increased basal melt, as only very narrow glaciers protected by wide stabilizing ridges cross the Transantarctic Mountains to enter this area (Morlighem et al., 2020). The Ross ice streams and glaciers feeding the Ronne ice shelf also experience limited dynamic response to increased basal melt. For the other regions, none of the CMIP5 forcing used predicted a large increase in ocean-induced melt by 2100, so we cannot conclude on the sensitivity of these sectors to oceanic forcings.

Figure 15

Sensitivity of individual basins to increased ocean basal melt over the 2015–2100 period: (a) the Antarctic Peninsula, (b) WAIS and (c) EAIS. The dynamic mass loss is approximated as to the total mass loss minus the cumulative anomaly in surface mass balance. It is shown as a function of cumulative ocean-induced basal melt anomaly over the same period for each of the 18 main Antarctic basins and for all RCP 8.5 experiments with medium ocean forcing. Dynamic change and basal melt are both relative to ctrl_proj experiment. Antarctic maps show the location of the 18 Antarctic basins.

[Figure omitted. See PDF]

The large spread in Antarctic ice sheet projections reported here contrasts with the relatively narrow range of projections reported as part of ISMIP6 in for the Greenland ice sheet. We attribute this difference to the dominant role of SMB in driving future evolution of Greenland and the more constrained forcing applied for ice front retreat in Greenland, in which most models used a prescribed a retreat rate.

For Antarctica, we find that uncertainties in the sea level estimates come from the spread in AOGCM forcing (see Sect. ), the melt parameterization adopted and its calibration (see Sects.  and ), and the spread caused by the choices made by the ice flow models for their initialization and the physical processes included (see Sect.  and ). All of these sources of uncertainty affect the results, and uncertainties in ocean conditions and their conversion into basal melt rates through parameterization lead to the largest spread of results, especially when different datasets are used for parameter calibration. Additional Antarctic mass losses of more than 20 cm SLR by 2100 under RCP 8.5 compared to constant climate conditions are reached only for the simulations based on the PIGL calibration (Fig. ) or as part of the open melt framework. Furthermore, not only does the magnitude of basal melt influence Antarctic dynamics, but the spatial distribution of melt rates has a strong effect on the results, as observed when comparing the open and standard experiments (Sect. ). These findings are similar to those described by based on idealized model configurations and highlight the need to acquire more observations and to use coupled ice–ocean models to better understand ice–ocean interactions and represent them in ice flow models .

The results presented here do not include any weighting of the ice flow models based on their agreements with observations or the number of simulations submitted. As explained in previous studies , the range of initialization techniques adopted by models leads to various biases. Some models are initialized with a long paleoclimate spin-up, giving limited spurious trends but an initial configuration further from the observed state, whereas other models initialized with data assimilation of present-day observations can capture these conditions accurately but often have nonphysical trend in their evolution. Assigning weights to different models is therefore a complicated question that is not addressed in the present study. This choice might lead to an overrepresentation of the models that submitted several contributions but is similar to that adopted within the larger CMIP framework.

The simulations performed as part of ISMIP6-Antarctica projections represent a significant improvement compared to previous intercomparisons of Antarctic evolution, especially in terms of the treatment of ice shelves, grounding line evolution, and ocean-induced basal melt that were not always included in previous continental Antarctic models . This progress is representative of improvements made to ice flow models over the past decade . Ice shelf melt parameterizations have been improved to reproduce the main features of basal melt simulated in ocean models and captured in observations. They are based on simulated ocean conditions extrapolated in ice shelf cavities, while uniform prescribed values were used in previous efforts . Grounding line migration and model resolution have been significantly improved (see Table ) and an increasing number of models are simulating ice front migrations. However, several limitations remain, regarding both external forcings and ice flow models . SMB forcing from AOGCMs generally has a coarse resolution, and no regional model was used to downscale the forcing, unlike what was done for Greenland , so SMB might not be well captured in regions with steep surface slopes. The inclusion of surface elevation feedbacks was left to the discretion of ice modeling groups, and no model included one, so this positive feedback was neglected in the present simulations. Because CMIP5 AOGCMs do not include ocean circulation under ice shelves, several simplifying assumptions must be made to estimate ocean conditions in ice shelf cavities . Ice–ocean interactions in ice shelf cavities are poorly observed and constrained , leading to additional limitations on the representation of ocean-induced sub-shelf melt. While pan-Antarctic estimates of basal melt have been produced , we are missing time series of basal melt at that scale and coinciding observations of oceanic conditions. Despite the progress in ice sheet numerical modeling over the last decade , significant limitations remain in our understanding of basal sliding , basal hydrology , calving and interaction with solid Earth . Finally, there was no incentive for models to represent the changes recently observed in Antarctica. However, as a variety of remote sensing observations are starting to provide time series of ice sheet changes over the recent past, it is becoming increasingly important to assess the ability of models to reproduce such observations in order to gain confidence in the projections.

The analysis of the simulations conducted as part of ISMIP6-Antarctic are projections that are presented here as relative to the ctrl_proj control experiments and therefore represent estimates of mass loss caused by variations in climate compared to a scenario with a constant climate. It was decided that using results of ice flow simulations directly, without subtracting the trend from a control run, is not yet appropriate given the large trend in the historical simulations and control experiments (Fig. ). Such a trend does not represent recent physical changes but rather limitations in observations , external forcings , ice flow models , and procedures used to initialize ice flow models . As ice sheets respond nonlinearly to changes, such an approach introduces a bias in the ice response, but this approach was deemed to be the most appropriate approach given current limitations. This same approach has been adopted in other recent ice flow modeling studies e.g.,. The choice of AOGCMs was made to cover a large range of responses to RCP scenarios but is not representative of the mean changes exhibited by CMIP5 AOGCMs . As a result, we expect that the spread of model response represented here covers the diversity of AOGCM outputs. However, computing mean values using different AOGCMs should be avoided, as only a few AOGCMs were sampled. Finally, all the results presented here are based on CMIP5 AOGCMs. Additional results based on CMIP6 AOGCMs will be presented in following publications.

6 Conclusions

Here we present simulations of the Antarctic ice sheet evolution between 2015 and 2100 from a multi-model ensemble, as part of the ISMIP6 framework. Ice sheet models from 13 international ice sheet modeling groups are forced with outputs from AOGCMs chosen to represent a large spread of possible evolution of oceanic and atmospheric conditions around Antarctica over the 21st century. Simulation results suggest that the Antarctic ice sheet could contribute between $-\mathrm{7.6}$ and 30.0 cm of SLE under the RCP 8.5 scenario compared to a scenario of constant conditions representative of the past decade. Climate models suggest significant increases in surface mass balance that are partially balanced by dynamic changes in response to ocean warming. Simulations suggest strong regional differences: WAIS loses mass under most scenarios and for all ice sheet models, as the increase in surface mass balance remains limited but the increase in ice discharge are large. EAIS, on the other hand, gains mass in many simulations, as dynamic mass loss is too limited to compensate for the large increase in surface mass balance. The regions most vulnerable to changes in the simulations are the Amundsen Sea sector in West Antarctica and Wilkes Land in East Antarctica. Simulations of the Antarctic ice sheet evolution under the RCP 2.6 scenario contribute less to sea level rise and have a smaller spread in SLE contribution between $-\mathrm{1.4}$ and 15.5 cm relative to a constant forcing, with less surface mass balance increase and a smaller dynamic response. The main sources of uncertainties highlighted in this study are the physics of ice flow models, the climate conditions used to force the ice sheet and the representation of ocean-induced melt at the base of ice shelves.

Appendix A Requested outputs

The model outputs requested as part of ISMIP6 are listed in Table . Annual values were submitted for both scalar and two-dimensional variables. Flux variables reported are averaged over calendar years, while state variables are reported at the end of calendar years.

Table A1

Data requests for Antarctica projections. ST stands for state variable, FL stands for flux variable and CST stands for constant.

Variable name	Type	Standard name	Unit
Ice sheet thickness	ST	land_ice_thickness	m
Ice sheet surface elevation	ST	surface_altitude	m
Ice sheet base elevation	ST	base_altitude	m
Bedrock elevation	ST	bedrock_altitude	m
Geothermal heat flux	CST	upward_geothermal_heat_flux_at_ground_level	W m${}^{-\mathrm{2}}$
Surface mass balance flux	FL	land_ice_surface_specific_mass_balance_flux	kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$
Basal mass balance flux	FL	land_ice_basal_specific_mass_balance_flux	kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$
Ice thickness imbalance	FL	tendency_of_land_ice_thickness	m s${}^{-\mathrm{1}}$
Surface velocity in $x$ direction	ST	land_ice_surface_x_velocity	m s${}^{-\mathrm{1}}$
Surface velocity in $y$ direction	ST	land_ice_surface_y_velocity	m s${}^{-\mathrm{1}}$
Surface velocity in $z$ direction	ST	land_ice_surface_upward_velocity	m s${}^{-\mathrm{1}}$
Basal velocity in $x$ direction	ST	land_ice_basal_x_velocity	m s${}^{-\mathrm{1}}$
Basal velocity in $y$ direction	ST	land_ice_basal_y_velocity	m s${}^{-\mathrm{1}}$
Basal velocity in $z$ direction	ST	land_ice_basal_upward_velocity	m s${}^{-\mathrm{1}}$
Mean velocity in $x$ direction	ST	land_ice_vertical_mean_x_velocity	m s${}^{-\mathrm{1}}$
Mean velocity in y direction	ST	land_ice_vertical_mean_y_velocity	m s${}^{-\mathrm{1}}$
Ice surface temperature	ST	temperature_at_ground_level_in_snow_or_firn	K
Ice basal temperature	ST	land_ice_basal_temperature	K
Magnitude of basal drag	ST	magnitude_of_land_ice_basal_drag	Pa
Land ice calving flux	FL	land_ice_specific_mass_flux_due_to_calving	kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$
Grounding line flux	FL	land_ice_specific_mass_flux_due_at_grounding_line	kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$
Land ice area fraction	ST	land_ice_area_fraction	1
Grounded ice sheet area fraction	ST	grounded_ice_sheet_area_fraction	1
Floating ice sheet area fraction	ST	floating_ice_sheet_area_fraction	1
Total ice sheet mass	ST	land_ice_mass	kg
Total ice sheet mass above floatation	ST	land_ice_mass_not_displacing_sea_water	kg
Area covered by grounded ice	ST	grounded_land_ice_area	m${}^{\mathrm{2}}$
Area covered by floating ice	ST	floating_ice_shelf_area	m${}^{\mathrm{2}}$
Total SMB flux	FL	tendency_of_land_ice_mass_due_to_surface_mass_balance	kg s${}^{-\mathrm{1}}$
Total BMB flux	FL	tendency_of_land_ice_mass_due_to_basal_mass_balance	kg s${}^{-\mathrm{1}}$
Total calving flux	FL	tendency_of_land_ice_mass_due_to_calving	kg s${}^{-\mathrm{1}}$
Total grounding line flux	FL	tendency_of_grounded_ice_mass	kg s${}^{-\mathrm{1}}$

Appendix B Summary of initial state and control run evolution

We report here the scalar values of simulated Antarctic ice mass, ice mass above floatation, ice extent, and ice shelf extent in Tables  and . Values are reported at the beginning of January 2015, when the experiments start in Table B1. We also report the evolution of ice mass, ice mass above floatation, ice extent and ice shelf extent during the ctrl_proj simulation (between 2015 and 2100) in Table .

Table B1

Simulated Antarctic ice mass, ice mass above floatation, total ice extent and floating ice extent at the beginning of the experiments (January 2015).

		Ice mass	Total	Floating
Model name	Ice mass	above floatation	ice extent	ice extent
	(10${}^{\mathrm{7}}$ Gt)	(10${}^{\mathrm{7}}$ Gt)	(10${}^{\mathrm{7}}$ km${}^{\mathrm{2}}$)	(10${}^{\mathrm{6}}$ km${}^{\mathrm{2}}$)
AWI_PISM_open	2.49	2.14	1.43	1.25
AWI_PISM_std	2.49	2.14	1.43	1.25
DOE_MALI_std	2.44	2.10	1.38	1.47
ILTS_PIK_SICOPOLIS_std	2.45	2.12	1.40	1.64
IMAU_IMAUICE1_std	2.32	1.99	1.41	1.51
IMAU_IMAUICE2_std	2.31	1.99	1.41	1.52
JPL1_ISSM_std	2.44	2.10	1.39	1.45
LSCE_GRISLI_std	2.47	2.13	1.40	1.46
NCAR_CISM_open	2.41	2.08	1.38	1.30
NCAR_CISM_std	2.41	2.08	1.38	1.30
PIK_PISM1_open	2.48	2.15	1.38	1.43
PIK_PISM2_open	2.49	2.15	1.39	1.44
UCIJPL_ISSM_open	2.40	2.08	1.36	1.47
UCIJPL_ISSM_std	2.40	2.08	1.36	1.47
ULB_fETISh_16_open	2.42	2.07	1.45	1.89
ULB_fETISh_16_std	2.42	2.07	1.45	1.92
ULB_fETISh_32_open	2.43	2.09	1.41	1.63
ULB_fETISh_32_std	2.43	2.09	1.42	1.70
UTAS_ELmerIce_std	2.43	2.09	1.41	1.35
VUB_AISMPALEO_std	2.49	2.14	1.42	1.19
VUW_PISM_open	2.43	2.08	1.39	1.34

Table B2

Simulated Antarctic ice mass, ice mass above floatation, total ice extent and floating ice extent change during the ctrl_proj experiment (between 2015 and 2100).

	Ice mass	Ice mass above	Total ice	Floating ice
Model name	change	floatation change	extent change	extent Change
	(Gt)	(Gt)	(10${}^{\mathrm{3}}$ km${}^{\mathrm{2}}$)	(10${}^{\mathrm{4}}$ km${}^{\mathrm{2}}$)
AWI_PISM_open	3394	$-\mathrm{1486}$	16.7	1.48
AWI_PISM_std	3394	$-\mathrm{1486}$	16.7	1.48
DOE_MALI_std	$-\mathrm{70}$ 394	$-\mathrm{51}$ 458	12.2	0.08
ILTS_PIK_SICOPOLIS_std	578	$-\mathrm{120}$	$-\mathrm{1.0}$	$-\mathrm{0.57}$
IMAU_IMAUICE1_std	$-\mathrm{10}$	$-\mathrm{22}$	0.0	0.21
IMAU_IMAUICE2_std	$-\mathrm{25}$ 564	$-\mathrm{17}$ 836	0.0	1.04
JPL1_ISSM_std	$-\mathrm{34}$ 450	$-\mathrm{33}$ 210	0.0	2.87
LSCE_GRISLI_std	3904	$-\mathrm{8972}$	56.2	8.25
NCAR_CISM_open	$-\mathrm{9126}$	$-\mathrm{4950}$	$-\mathrm{0.9}$	0.75
NCAR_CISM_std	548	122	$-\mathrm{0.2}$	$-\mathrm{0.00}$
PIK_PISM1_open	$-\mathrm{22}$ 374	$-\mathrm{5324}$	$-\mathrm{31.9}$	$-\mathrm{1.15}$
PIK_PISM2_open	2432	1826	4.5	0.27
UCIJPL_ISSM_open	12 594	$-\mathrm{9308}$	0.0	5.47
UCIJPL_ISSM_std	43 258	9080	0.0	9.03
ULB_fETISh_16_open	$-\mathrm{83}$ 960	$-\mathrm{39}$ 872	$-\mathrm{92.7}$	$-\mathrm{0.55}$
ULB_fETISh_16_std	$-\mathrm{22}$ 352	$-\mathrm{9850}$	5.5	$-\mathrm{7.45}$
ULB_fETISh_32_open	$-\mathrm{84}$ 112	$-\mathrm{12}$ 830	$-\mathrm{83.8}$	$-\mathrm{7.94}$
ULB_fETISh_32_std	52 896	47 080	13.3	$-\mathrm{9.27}$
UTAS_ELmerIce_std	58 810	13 380	0.0	$-\mathrm{16.97}$
VUB_AISMPALEO_std	$-\mathrm{20}$ 124	$-\mathrm{7970}$	$-\mathrm{2.4}$	0.89
VUW_PISM_open	$-\mathrm{1680}$	$-\mathrm{5102}$	141.8	14.30

Appendix C Ice flow model initialization and characteristics

The descriptions below summarize the initialization procedure and main characteristics of the different ice flow modeling groups.

C1

AWI_PISM

The AWI_PISM ice sheet model is based on the Parallel Ice Sheet Model PISM, version 1.1.4 with modifications for ISMIP6. PISM solves a hybrid combination of the non-sliding shallow ice approximation (SIA) and the shallow shelf approximation (SSA) for grounded ice, where the SSA solution acts as a sliding law, and only the SSA for floating ice. PISM also solves for enthalpy to account for the temperature and water content of the ice in the rheology. The model uses a structured rectangular grid with a uniform horizontal resolution of 8 km (16 km early in the spin-up) and 81 vertical $z$-coordinate levels that are refined towards the base. The total ice domain height is 6000 m with an additional heat conducting bedrock layer of 2000 m thickness (21 equal levels). The calving front can evolve freely on a sub-grid scale . In addition to calving below a certain thickness threshold (here 150 m), a kinematic first-order calving law, called eigen-calving , is utilized with the calving parameter $K={\mathrm{10}}^{\mathrm{17}}\mathrm{m}\mathrm{s}$. Floating ice that extends far into the open ocean (seafloor elevation reaches 2000 m below sea level) is also calved off. The grounding line position is determined using hydrostatic equilibrium. Basal friction in partially grounded cells is weighted according to the grounded area fraction . The nonlocal quadratic melt scheme and the related datasets provided by ISMIP6 are used to compute the ice shelf basal melt in the spin-up and all “standard” experiments. For the “open” experiments, the local quadratic melt scheme is used. Ice shelf basal melt is applied on a sub-grid scale.

To initialize the model, an equilibrium-type spin-up based on steady present-day climate has been performed. Atmospheric forcing (2 m air temperature and precipitation) is the multi-annual mean 1995–2014 (ISMIP6 reference period) from RACMO2.3p2 . For the surface mass balance, a positive degree-day scheme is used. Geothermal heat flux is from and the bedrock elevation is fixed in time. The ocean is forced with the present-day ocean forcing field provided by ISMIP6. The spin-up consists of an initialization with idealized temperature–depth profiles, a 100-year geometry relaxation run and a 200 kyr thermo-mechanically coupled run with fixed geometry for thermal equilibration. For those stages, the non-sliding SIA is used on a 16 km horizontal grid. After re-gridding the output (except the geometry) onto the final 8 km grid, the model runs for 30 kyr using full model physics and a freely evolving geometry. The initial ice sheet geometry for the spin-up is based on Bedmap2 and is refined in the Recovery Glacier area with additional ice thickness datasets . The historical simulation from January 2005 until the end of December 2014 employs the NorESM1-M-RCP8.5 atmospheric and oceanic forcing.

C2

DOE_MALI

MPAS-Albany Land Ice (MALI) uses a three-dimensional, first-order “Blatter–Pattyn” momentum balance solver solved using finite-element methods . Ice velocity is solved on a two-dimensional map plane triangulation extruded vertically to form tetrahedra. Mass and tracer transport occur on the Voronoi dual mesh using a mass-conserving finite volume first-order upwinding scheme. Mesh resolution is 2 km along grounding lines, in all marine regions of West Antarctica and in marine regions of East Antarctica where present-day ice thickness is less than 2500 m to ensure that the grounding line remains in the fine-resolution region even under full retreat of West Antarctica and large parts of East Antarctica. Mesh resolution coarsens to 20 km in the ice sheet interior and no greater than 6 km in the large ice shelves. The horizontal mesh has 1.6 million cells. The mesh uses 10 vertical layers that are finest near the bed (4 % of total thickness in deepest layer) and coarsen towards the surface (23 % of total thickness in shallowest layer). Ice temperature is based on results from and held fixed in time. The model uses a linear basal friction law with spatially varying basal friction coefficient. The basal friction of grounded ice and the viscosity of floating ice are inferred to best match observed surface velocity using an adjoint-based optimization method and then kept constant in time. The grounding line position is determined using hydrostatic equilibrium, with sub-element parameterization of the friction. Sub-ice-shelf melt rates come from and are extrapolated across the entire model domain to provide nonzero ice shelf melt rates after grounding line retreat. The surface mass balance is from the RACMO2.1 1979–2010 mean . Maps of surface and basal mass balance forcing are kept constant with time in the ctrl_proj experiment. Time-varying anomalies of surface and basal mass balance relative to the original fields are applied in all other experiments. The ice front position is fixed at the extent of the present-day ice sheet. After initialization, the model is relaxed for 99 years so that the geometry and grounding lines can adjust.

C3

ILTS_PIK_SICOPOLIS

The model SICOPOLIS version 5.1 (Greve and SICOPOLIS Developer Team, 2019; http://www.sicopolis.net/, last access: 6 July 2020) is applied to the Antarctic ice sheet with hybrid shallow-ice–shelfy-stream dynamics for grounded ice and shallow-shelf dynamics for floating ice. Ice thermodynamics is treated with the melting-CTS (cold-temperate transition surface) enthalpy method (ENTM) by . The ice surface is assumed to be traction-free. Basal sliding under grounded ice is described by a Weertman–Budd-type sliding law with sub-melt sliding and subglacial hydrology . The model is initialized by a paleoclimatic spin-up over 140 000 years until 1990, forced by Vostok $\mathit{\delta }D$ converted to $\mathrm{\Delta }T$ , in which the topography is nudged towards the present-day topography to enforce a good agreement . The basal sliding coefficient is determined individually for the 18 IMBIE-2016 basins by minimizing the RMSD between simulated and observed logarithmic surface velocities. The historical run from 1990 until 2015 employs the NorESM1-M-RCP8.5 atmospheric and oceanic forcing. For the last 2000 years of the spin-up, the historical run and the future climate simulations, a regular (structured) grid with 8 km resolution is used. In the vertical, we use terrain-following coordinates with 81 layers in the ice domain and 41 layers in the thermal lithosphere layer below. The present-day surface temperature is parameterized , the present-day precipitation is by and , and runoff is modeled by the positive-degree-day method with the parameters by . The 1960–1989 average SMB correction that results diagnostically from the nudging technique is used as a prescribed SMB correction for the future climate simulations. The bed topography is Bedmap2 , the geothermal heat flux is by , and isostatic adjustment is included using an elastic-lithosphere–relaxing-asthenosphere (ELRA) model (parameters by ). Present-day ice-shelf basal melting is parameterized by the ISMIP6 standard approach (Eq. ). A more detailed description of the setup, which is consistent with the one used for the LARMIP-2 and ABUMIP initiatives, will be given elsewhere .

C4

IMAU_IMAUICE

The finite-difference model uses a combination of SIA and SSA solutions, with velocities added over grounded ice to model basal sliding . The model grid at 32 km horizontal resolution covers the entire Antarctic ice sheet and surrounding ice shelves. The grounded ice margin is freely evolving, while the shelf extends to the grid margin and a calving front is not explicitly determined. We use the Schoof flux boundary condition at the grounding line with a heuristic rule following . For the ISMIP6 projections the sea level equation is not solved or coupled . We run the thermodynamically coupled model with constant present-day boundary conditions to determine a thermodynamic steady state. The model is first initialized for 100 kyr using the average 1979–2014 SMB and surface ice temperature from RACMO 2.3 . Bedrock elevation is fixed in time with data taken from the Bedmap2 dataset , and geothermal heat flux data are from . We then run this for 30 kyr with constant ice temperature from the first run to get to a dynamic steady state, which was our initial condition for initMIP. For IMAUICE1 we assign this steady state to the year 1978 and run the historical period 1979–2014 unforced, keeping the initial SMB constant and sub-shelf basal melting at zero. This model setup is provided for comparison with initMIP. For IMAUICE2 we assign the steady state to the year 1900 and run a 79-year experiment with constant SMB and sub-shelf basal melt rates estimated for the modeled ice draft at 1900 using the shelf melt parameterization of with a thermal forcing derived from the World Ocean Atlas (WOA) at 400 m depth. We continue with the historical period 1979–2014, keeping the initial sub-shelf basal melt rates constant, with transient SMB variations from RACMO 2.3 .

C5

JPL_ISSM

The JPL_ISSM ice sheet model configuration relies on data assimilation of present-day conditions, followed by a short model relaxation as described in . The model domain covers the present-day Antarctic ice sheet, and its geometry is based on an early version of BedMachine Antarctica . The model is based on the 2D Shelfy-Stream Approximation , and the mesh resolution varies between 1 km along the coast to 50 km in the interior and has a resolution of 8 km or finer within the boundary of all initial ice shelves. The model is vertically extruded into 15 layers. To estimate land ice viscosity ($B$), we compute the ice temperature based on a thermal steady state using three-dimensional higher-order stress balance equations, observations of surface velocities , and basal friction inferred from surface elevations . Thermal boundary conditions are geothermal heat flux from and surface temperatures from . Steady-state ice temperatures are then vertically averaged and used to calibrate the ice viscosity, which is held constant over time. To infer the unknown basal friction coefficient over grounded ice and the ice viscosity of the floating ice, we use data assimilation to reproduce observed surface velocities from . Following this, we run the model forward for 2 years, allowing the grounding line position and ice geometry to relax . The grounding line evolves assuming hydrostatic equilibrium and following a sub-element grid scheme SEP2 in. The ice front remains fixed in time during all simulations performed, and we impose a minimum ice thickness of 1 m everywhere in the domain. The surface mass balance and the ice shelf basal melt rates used in the control experiment are, respectively, from the 1979–2010 mean of RACMO2.1 and from the 2004–2013 mean following .

C6

LSCE_GRISLI

The GRISLI model is a three-dimensional thermo-mechanically coupled ice sheet model originating from the coupling of the inland ice model of and and the ice shelf model of , extended to the case of ice streams treated as dragging ice shelves . In the version used here, over the whole domain, the velocity field consists of the superposition of the shallow-ice approximation (SIA) velocities for ice flow due to vertical shearing and the shallow-shelf approximation (SSA) velocities, which are used as a sliding law . For the initMIP-Antarctica experiments, we used the GRISLI version 2.0 , which includes the analytical formulation of to compute the flux at the grounding line. Basal drag is computed with a power law basal friction . For this study, we use an iterative inversion method to infer a spatially variable basal drag coefficient that insures an ice thickness that is as close as possible to observations with a minimal model drift . The basal drag is assumed to be constant for the forward experiments.

The model uses finite differences on a staggered Arakawa C grid in the horizontal plane at 16 km resolution with 21 vertical levels. Atmospheric forcing, namely near-surface air temperature and surface mass balance, is taken from the 1979–2016 climatological annual mean computed by RACMO2.3p2 regional atmospheric model . Sub-shelf basal melting rates are computed with the nonlocal quadratic parametrization suggested in ISMIP. For the inversion step and the control experiments we use the 1995–2017 climatological observed thermal forcing. The initial ice sheet geometry, bedrock, and ice thickness are taken from the Bedmap2 dataset , and the geothermal heat flux is from .

C7

NCAR_CISM

The Community Ice Sheet Model CISM, uses finite-element methods to solve a depth-integrated higher-order approximation over the entire Antarctic ice sheet. The model uses a structured rectangular grid with uniform horizontal resolution of 4 km and five vertical $\mathit{\sigma }$–coordinate levels. The ice sheet is initialized with present-day geometry and an idealized temperature profile, then spun up for 30 000 years using 1979–2016 climatological surface mass balance and surface air temperature from RACMO2.3 . During the spin-up, basal friction parameters (for grounded ice) and sub-shelf melt rates (for floating ice) are adjusted to nudge the ice thickness during present-day observations. This method is a hybrid approach between assimilation and spin-up, similar to that described by . The geothermal heat flux is taken from . The basal sliding is similar to that of , combining power law and Coulomb behavior. The grounding line location is determined using hydrostatic equilibrium and sub-element parameterization . Basal melt is applied in partly floating grid cells in proportion to the floating fraction as determined by the grounding-line parameterization. The calving front is initialized from present-day observations and thereafter is allowed to retreat but not advance. For the historical run (1995–2014), the SMB anomaly was provided by RACMO2.3, and the basal melt rate anomaly was derived from NorESM1-M RCP8.5 thermal forcing. For the open parameterization of basal melting, we weighted the melt from the standard nonlocal parameterization by $\sin \mathit{\theta }$, where $\mathit{\theta }$ is the ice shelf basal slope angle, with ${\mathit{\gamma }}_{\mathrm{0}}$ recalibrated by N. Jourdain. See for more information about the model.

C8

PIK_PISM

With the Parallel Ice Sheet Model (PISM, , https://pism-docs.org/wiki/doku.php, last access: 8 July 2020; version 1.0, pism version available at 10.5281/zenodo.3903343), we perform an equilibrium simulation on a regular rectangular grid with 8 km horizontal resolution. The vertical resolution increases from 100 m at the top of the domain to 13 m at the (ice) base, with a domain height of 6000 m. PISM uses a hybrid of the Shallow-Ice Approximation (SIA) and the two-dimensional Shelfy-Stream Approximation of the stress balance SSA, over the entire Antarctic ice sheet. The grounding line position is determined using hydrostatic equilibrium, with sub-grid interpolation of the friction at the grounding line . The calving front position can freely evolve using the eigen-calving parameterization . PISM is a thermomechanically coupled (polythermal) model based on the Glen–Paterson–Budd–Lliboutry–Duval flow law . The three-dimensional enthalpy field can evolve freely for given boundary conditions.

The model is initialized from Bedmap2 geometry , with surface mass balance and surface temperatures from RACMOv2.3p2 1986–2005 mean remapped from 27 km resolution. Geothermal heat flux is from . We use the Potsdam Ice-shelf Cavity model PICO,, which extends the ocean box model by for application in three-dimensional ice sheet models to calculate basal melt rate patterns underneath the ice shelves. We use a compilation of observed ocean temperature and salinity values (1975–2012, ; 1955–2018, ) to drive PICO. We apply a power law for sliding with a Mohr–Coulomb criterion relating the yield stress to parameterized till material properties and the effective pressure of the overlaying ice on the saturated till . Basal friction and sub-shelf melting are linearly interpolated on a sub-grid scale around the grounding line . We apply eigen-calving in combination with the removal of all ice that is thinner than $\mathrm{50}$ $\mathrm{m}$ or extends beyond present-day ice fronts .

C9

UCIJPL_ISSM

We initialize the model by using data assimilation of present-day conditions, following the method presented in . The mesh horizontal resolution varies from 3 km near the margins to 30 km inland. The mesh is vertically extruded into 10 layers. We use a higher-order stress balance and an enthalpy-based thermal model . The initialization is a two-step process: we first invert for ice shelf viscosity ($B$) and then invert for basal friction under grounded ice assuming thermo-mechanical steady state. Our geometry is based on BedMachine Antarctica . The thermal model is constrained by surface temperatures from and geothermal heat flux from , both of which are included in the SeaRISE dataset . The surface mass balance used in the control experiment is from RACMO 2.3 .

C10

ULB_FETISH

The f.ETISh (Fast Elementary Thermomechanical Ice Sheet) model version 1.3 is a vertically integrated hybrid finite-difference (SSA for basal sliding; SIA for grounded ice deformation) ice sheet–ice shelf model with vertically integrated thermomechanical coupling. The transient englacial temperature field is calculated in a 3D fashion. The marine boundary is represented by a grounding-line flux condition according to , coherent with a power law basal sliding (power law coefficient of 2). Model initialization is based on an adapted iterative procedure based on to fit the model as closely as possible to present-day observed thickness and flow field . The model is forced by present-day surface mass balance and temperature , based on the output of the regional atmospheric climate model RACMO2 for the period 1979–2011. The PICO model was employed to calculate sub-shelf melt rates, based on present-day observed ocean temperature and salinity , onto which the initMIP forcings for the different basins are added. The model is run on a regular grid of 16 km with time steps of 0.05 year.

C11

UTAS_ElmerIce

The Elmer/Ice model domain covers the present-day Antarctic ice sheet, and its geometry is interpolated from the Bedmap2 dataset . An unstructured mesh in the horizontal is refined using the Hessian of the observed surface velocity, as in . Mesh resolution in the horizontal varies from approximately 4 km near the grounding lines of fast-flowing ice streams to approximately 40 km in the interior. The mesh is extruded to 10 layers in the vertical. The forward simulations solve the Stokes equations directly . Initialization was comprised of the following steps:

• short surface relaxation (20 time steps of 0.001 years);

• inversion for sliding coefficient with constant temperature $T=-\mathrm{20}$ ${}^{\circ }$ ;

• steady-state temperature simulation using the flow field from previous step;

• inversion for sliding coefficient using the new temperature field from the previous step;

• thermo-mechanically coupled steady state temperature–velocity calculation, using the basal sliding coefficient distribution from the previous step;

• inversion for sliding coefficient using the latest temperature field from the previous step;

• surface relaxation (10 years with an increasing time step size).

A linear sliding relation is used, and the ice front is not allowed to evolve. Elmer/Ice solves a contact problem at the grounding line, and no further parameterizations are applied. Thermal boundary conditions are the geothermal heat flux from and surface temperatures from . Steady temperature is solved for during the initialization steps and held constant during the transient simulations. We impose a minimum ice thickness of 40 m everywhere in the domain. The surface mass balance used in the surface relaxation and control experiment is the 1995 to 2014 mean from the MAR model . Basal melt rates are computed using the local quadratic parameterization provided by ISMIP as an alternative to the nonlocal parameterization.

C12

VUB_AISMPALEO

The Antarctic ice sheet model from the Vrije Universiteit Brussel is derived from the coarse-resolution version used mainly in simulations of the glacial cycles . It considers thermomechanically coupled flow in both the ice sheet and the ice shelf, using the SIA and SSA coupled across a transition zone one grid cell wide. Basal sliding is calculated using a Weertman relation inversely proportional to the height above buoyancy wherever the ice is at the pressure melting point. The horizontal resolution is 20 km, and there are 31 layers in the vertical. The model is initialized with a freely evolving geometry until a steady state is reached. The precipitation pattern is based on the compilation used in , updated with accumulation rates obtained from shallow ice cores during the EPICA pre-site surveys . Surface melting is calculated over the entire model domain with the Positive Degree Day (PDD) scheme, including meltwater retention by refreezing and capillary forces in the snowpack . The sub-shelf basal melt rate is parameterized as a function of local mid-depth (485–700 m) ocean water temperature above the freezing point . A distinction is made between protected ice shelves (Ross and Filchner-Ronne) with a low melt factor and all other ice shelves with a higher melt factor. Ocean temperatures are derived from the LOVECLIM climate model , and melt is parameterized with a plume model . Heat conduction is calculated in a slab of bedrock 4 km thick underneath the ice sheet. Isostatic compensation is based on an elastic lithosphere floating on a viscous asthenosphere (ELRA model) but is not allowed to evolve further in line with the initMIP-Antarctica experiments.

C13

VUW_PISM

We use an identical approach to the one described in . Starting from initial bedrock and ice thickness conditions from , together with reference climatology from , we run a multistage spinup that guarantees well-evolved thermal and dynamic conditions without loss of accuracy in terms of geometry. This is achieved through an iterative nudging procedure, 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. In summary, 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 1000 years using full model physics and a present-day climate, refine the grid to 10 km for a further 500 years and then refine the grid to 8 km for a GCM-forced historical run from 1950 to 2000. The resultant configuration is then used as the starting point for each of our forward experiments.

Data availability

Model outputs from the simulations described in this paper will be made available in the CMIP6 archive through the Earth System Grid Federation (ESGF; https://esgf-node.llnl.gov/search/cmip6/ , last access: 11 July 2020) for two-dimensional variables. Scalars computed from two-dimensional fields for this study will be available from archive on Zenodo with the following doi: 10.5281/zenodo.3940766. In order to document CMIP6's scientific impact and enable ongoing support of CMIP, users are obligated to acknowledge CMIP6, participating modeling groups and the ESGF centers (see details on the CMIP Panel website: https://www.wcrp-climate.org/wgcm-cmip, last access: 11 July 2020). The forcing datasets are available through the ISMIP6 wiki (http://www.climate-cryosphere.org/wiki/index.php?title=ISMIP6_wiki_page, ) and will also be archived in a publicly available repository (see assets tab).

Code availability

Data processing, analysis and plotting scripts are archived in permanent repositories on Zenodo and will be available via the following digital object identifier: 10.5281/zenodo.3940768.

Author contributions

HS, SN, AJP, HG and WHL designed the experiments. CA, XA-D, AB, RiC, TH, NCJ, CML, ES, RSS, FS and LDT derived the external forcings. HS, HG, WHL, TA, ReC, CD, BKG-F, RG, NG, RG, MJH, AH, PH, TK, GRL, DPL, MM, FP, TP, SFP, AQ, RR, N-JS, AA, JVB, RSWvdW, RW, CZ, ToZ and TaZ ran the ice flow model simulations. HS analyzed the results with inputs from all authors and wrote the first draft of the manuscript. All authors contributed to the writing of the manuscript.

Competing interests

Eric Larour serves as topical editor for the journal. William Lipscomb, Sophie Nowicki, Helene Seroussi, Ayako Abe-Ouchi and Robin Smith are editors of the Special Issue “The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6)”.

Special issue statement

This article is part of the special issue “The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6)”. It is not associated with a conference.

Acknowledgements

We thank the Climate and Cryosphere (CliC) effort, which provided support for ISMIP6 through sponsoring of workshops, hosting the ISMIP6 website and wiki, and promoting ISMIP6. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP5 and CMIP6. We thank the climate modeling groups for producing their model output and making it available; the Earth System Grid Federation (ESGF) for archiving the CMIP data and providing access to it; the University at Buffalo for ISMIP6 data distribution and upload; and the multiple funding agencies who support CMIP5, CMIP6, and ESGF. We thank the ISMIP6 steering committee, the ISMIP6 model selection group and ISMIP6 dataset preparation group for their continuous engagement in defining ISMIP6. This is ISMIP6 contribution no. 11.

Research was carried out at the Jet Propulsion Laboratory, California Institute of Technology. Helene Seroussi and Nicole Schlegel are supported by grants from NASA Cryospheric Science and Modeling, Analysis, and Predictions Programs. AB was supported by the U.S. Department of Energy (DOE) Office of Science Regional and Global Model Analysis (RGMA) component of the Earth and Environmental System Modeling (EESM) program (HiLAT-RASM project), and the DOE Office of Science (Biological and Environmental Research), Early Career Research program. Heiko Goelzer has received funding from the program of the Netherlands Earth System Science Centre (NESSC), financially supported by the Dutch Ministry of Education, Culture and Science (OCW) under grant no. 024.002.001. Rupert Gladstone and Thomas Zwinger were supported by Academy of Finland grant nos. 286587 and 322430. Chen Zhao and Ben Galton-Fenzi were supported under the Australian Research Council's Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001) and received grant funding from the Australian Government for the Australian Antarctic Program Partnership (Project ID ASCI000002). Support for Xylar Asay-Davis, Matthew Hoffman, Stephen Price and Tong Zhang was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the US Department of Energy (DOE), Office of Science, Advanced Scientific Computing Research, and Biological and Environmental Research Programs. MALI Earth System Grid Federation simulations used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science user facility supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC02-05CH11231. Nicolas Jourdain is funded by the French National Research Agency (ANR) through the TROIS-AS project (ANR-15-CE01-0005-01) and the European Commission through the TiPACCs project (grant no. 820575, call H2020-LC-CLA-2018-2). Philippe Huybrechts and Jonas Van Breedam acknowledge support from the iceMOD project funded by the Research Foundation - Flanders (FWO-Vlaanderen). Ralf Greve was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (grant nos. JP16H02224, JP17H06104 and JP17H06323). Support for Nicholas Golledge and Daniel Lowry was provided by the New Zealand Ministry of Business Innovation and Employment contract no. RTVU1705. The work of Thomas Kleiner has been conducted in the framework of the PalMod project (FKZ: 01LP1511B), supported by the German Federal Ministry of Education and Research (BMBF) as part of the Research for Sustainability initiative (FONA). Support for Mathieu Morlighem and Tyler Pelle was provided by the National Science Foundation (NSF, grant no. 1739031). Development of PISM is supported by NASA (grant no. NNX17AG65G) and the NSF (grant nos. PLR-1603799 and PLR-1644277). Luke Trusel was supported under NSF Antarctic Glaciology Program award no. 1643733. The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Federal State of Brandenburg for supporting this project by providing resources on the high-performance computer system at the Potsdam Institute for Climate Impact Research. Computer resources for this project have been also provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre (https://www.lrz.de/, last access: 8 July 2020) under Project ID pr94ga and pn69ru. Ronja Reese was supported by the Deutsche Forschungsgemeinschaft (DFG) under grant no. WI 4556/3-1 and through the TiPACCs project, which receives funding from the European Union’s Horizon 2020 Research and Innovation Programme under grant agreement no. 820575. Torsten Albrecht is supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the priority program “Antarctic Research with comparative investigations in Arctic ice areas” by grant no. WI4556/2-1. Reinhard Calov was funded by the Bundesministerium für Bildung und Forschung (BMBF) grants PalMod-1.1 and PalMod-1.3. Gunter Leguy and William Lipscomb were supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement no. 1852977. Computing and data storage resources for CISM simulations, including the Cheyenne supercomputer (10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. Funding support for Nicholas Golledge and Daniel Lowry was provided by the New Zealand Ministry of Business, Innovation and Employment through Victoria University of Wellington (RTUV1705), the Antarctic Science Platform (ANTA1801) and the Royal Society of New Zealand (grant no. RDF-VUW1501).

Financial support

This research has been supported by the U.S. Department of Energy, Office of Science, the Netherlands Earth System Science Centre (grant no. 024.002.001), the Academy of Finland (grant nos. 286587 and 322430), the Australian Research Council (grant no. SR140300001), the Agence Nationale de la Recherche (grant no. ANR-15-CE01-0005-01), the European Commission (TiPACCs grant no. 820575), the Research Foundation – Flanders, the Japan Society for the Promotion of Science (grant nos. JP16H02224, JP17H06104 and JP17H06323), the New Zealand Ministry of Business Innovation and Employment (grant no. RTVU1705), the German Federal Ministry of Education and Research, the Office of Polar Programs (grant no. 1739031), the National Science Foundation (grant nos. 1603799, 1644277,1852977, and 1916566), the National Aeronautics and Space Administration (grant nos. NNX17AG65G and NNX17AI03G), the Deutsche Forschungsgemeinschaft (grant nos. WI4556/2-1 and WI4556/3-1), and the Norwegian Research Council (grant nos. 280727 and 295075).

Review statement

This paper was edited by Christina Hulbe and reviewed by two anonymous referees.