Evolution of the Antarctic Ice Sheet Over the Next Three Centuries From an ISMIP6 Model Ensemble

Abstract

The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) is the primary effort of CMIP6 (Coupled Model Intercomparison Project–Phase 6) focusing on ice sheets, designed to provide an ensemble of process‐based projections of the ice‐sheet contribution to sea‐level rise over the twenty‐first century. However, the behavior of the Antarctic Ice Sheet beyond 2100 remains largely unknown: several instability mechanisms can develop on longer time scales, potentially destabilizing large parts of Antarctica. Projections of Antarctic Ice Sheet evolution until 2300 are presented here, using an ensemble of 16 ice‐flow models and forcing from global climate models. Under high‐emission scenarios, the Antarctic sea‐level contribution is limited to less than 30 cm sea‐level equivalent (SLE) by 2100, but increases rapidly thereafter to reach up to 4.4 m SLE by 2300. Simulations including ice‐shelf collapse lead to an additional 1.1 m SLE on average by 2300, and can reach 6.9 m SLE. Widespread retreat is observed on that timescale in most West Antarctic basins, leading to a collapse of large sectors of West Antarctica by 2300 in 30%–40% of the ensemble. While the onset date of retreat varies among ice models, the rate of upstream propagation is highly consistent once retreat begins. Calculations of sea‐level contribution including water density corrections lead to an additional ∼10% sea level and up to 50% for contributions accounting for bedrock uplift in response to ice loading. Overall, these results highlight large sea‐level contributions from Antarctica and suggest that the choice of ice sheet model remains the leading source of uncertainty in multi‐century projections.

Introduction

The Antarctic Ice Sheet has undergone rapid acceleration, thinning, and grounding-line retreat over the past few decades, raising global mean sea level by 14 mm during 1979–2017 (Rignot et al., 2019; Shepherd et al., 2019). Estimates of the Antarctic Ice Sheet contribution to sea-level over the twenty-first century vary from a few millimeters to more than 1 m sea-level equivalent (SLE, see Table A1 for the list of all acronyms used in the text, tables and figures) and is the largest source of uncertainty in sea level projections (Cornford et al., 2016; DeConto & Pollard, 2016; DeConto et al., 2021; Edwards et al., 2021; Ritz et al., 2015; Schlegel et al., 2018). Simulations from the Ice Sheet Model Intercomparison Project for CMIP6 (Coupled Model Intercomparison Project–phase 6) (Nowicki et al., 2016, 2020) suggested a sea-level contribution from the Antarctic Ice Sheet between −5 and 43 cm SLE by 2100 in addition to the current background trend in response to past climate warming (Payne et al., 2021; Seroussi et al., 2020). This projection was based on an ensemble of runs using 13 different ice sheet models, all forced with similar oceanic and atmospheric conditions derived from the Coupled Model Intercomparison Project–phase 5 (CMIP5) and CMIP6 simulations (Barthel et al., 2020; Jourdain et al., 2020; Payne et al., 2021). The results suggested that the large uncertainty in Antarctic evolution to 2100 is dominated by the choice of ice-flow models, the parameters they use, and initialization methods, but that the impact of the climate forcing increases steadily until the end of the century (Seroussi et al., 2023). Rapid mass loss associated with potential instability mechanisms, such as the Marine Ice Sheet Instability (Schoof, 2007; Thomas & Bentley, 1978; Weertman, 1974) and the Marine Ice Cliff Instability (DeConto & Pollard, 2016; Pollard et al., 2015), do not play a major role through 2100, but have been suggested to possibly happen after 2100 (DeConto et al., 2021).

Some previous studies of multi-century ice-sheet evolution suggest that the Antarctic Ice Sheet could contribute as much as 10 m SLE by 2300 under high-emission scenarios (DeConto & Pollard, 2016; DeConto et al., 2021), while others suggest a contribution of no more than 4 m SLE (Bindschadler et al., 2013; Chambers et al., 2022; Coulon et al., 2023; Golledge et al., 2015; Greve et al., 2023; Klose et al., 2023; Lipscomb et al., 2021; Lowry et al., 2021). The impact of the climate scenario, as represented by the choice of Representative Concentration Pathways, is limited during the twenty-first century, but starts to emerge around 2150 and quickly becomes an important driver of the Antarctic Ice Sheet sea-level contribution. Simulations from Lowry et al. (2021) point to a multi-meter gap between emission scenarios by 2300, from a mean of 1.0 m for the low-emission scenario RCP2.6 to 3.7 m for the high-emission scenario RCP8.5. However, the climate forcing in these previous experiments was not based on conditions simulated by climate models for the entire period, but rather on conditions simulated for the twenty-first century and then extended beyond 2100. The forcing after 2100 was either held constant at conditions from the end of the twenty-first century (Bindschadler et al., 2013; Chambers et al., 2022; Lipscomb et al., 2021; Lowry et al., 2021), held at constant values after a certain atmospheric temperature threshold was reached (DeConto et al., 2021), or extrapolated based on the twenty-first century climate combined with indices from climate models run until 2300 (Greve et al., 2023). As the atmosphere and Southern Ocean respond differently and on different time scales to climate warming, such extensions likely introduce shortcomings in the applied forcings. A number of global climate models with oceanic and atmospheric components have been run until 2300 under low- and high-emission scenarios as part of CMIP5 and CMIP6 (Tebaldi et al., 2021), creating an opportunity to use consistent forcing for ice-flow models over this longer period.

The extension of the shared socio-economic pathway scenarios beyond 2100 is described in Meinshausen et al. (2020). Under SSP5-8.5, the emissions of greenhouse gases increase until the second half of the twenty-first century, then decrease linearly from 2100 to zero or very low values in 2250. The extended RCP8.5 (Representative Concentration Pathway) scenario is similar in terms of carbon emissions, but keeps higher values of other greenhouse gases until 2300. The extended SSP1-2.6 and RCP2.6 scenarios correspond to emissions that stabilize to very low values after 2100. Under the extended SSP5-8.5 scenario, the Antarctic climate is much warmer than today in 2300, with almost no sea ice remaining and annual surface air temperatures up to 20°C warmer in some models (Mathiot & Jourdain, 2023).

In this study, we investigate the Antarctic Ice Sheet response to warming climate conditions until 2300 from an ensemble of 16 ice-flow models, using atmospheric and oceanic forcings derived from global climate model simulations performed for CMIP5 and CMIP6 (Barthel et al., 2020; Eyring et al., 2016; Knutti et al., 2013) in a framework similar to previous ISMIP6 studies (Goelzer, Nowicki, et al., 2020; Payne et al., 2021; Seroussi et al., 2020). We refer to the previous Antarctic Ice Sheet model ensemble (Seroussi et al., 2020) as “ISMIP6 Antarctica 2100” and the new ensemble as “ISMIP6 Antarctica 2300.” Some experiments in the new ensemble are based on climate forcing derived from simulations run until 2300, while others are based on 2080–2100 conditions maintained beyond 2100, to compare the results with previous studies and investigate the impact of this choice. Most experiments are based on high-emission scenarios to assess vulnerable basins and provide upper bounds on sea level contributions. Several experiments are also designed to investigate the impact of ice-shelf collapse and low-emission scenarios. We describe the climate forcings and experimental protocol in Section 2 and list the ice flow simulations and model characteristics in Section 3. We analyze and discuss the resulting large-scale and local ice sheet changes and SLE contribution in Sections 4 and 5. We finish with general conclusions and implications for future research.

Climate Forcings

The forcing for ice sheet models is derived from selected global climate model simulations following the approach used in ISMIP6 Antarctica 2100 (Barthel et al., 2020; Jourdain et al., 2020; Nowicki et al., 2020; Seroussi et al., 2020). This forcing includes atmospheric and oceanic forcings, as well as prescribed ice-shelf collapse. This section summarizes the generation of these different forcings and the experiments performed for ISMIP6 Antarctica 2300. More details on the forcing and climate model choices for ISMIP6 Antarctica 2100 can be found in Barthel et al. (2020), Jourdain et al. (2020), and Nowicki et al. (2020).

Selection of CMIP5 and CMIP6 Climate Models

Four global climate models are used to generate forcings for the ice sheet models: two CMIP5 models, the Community Climate System Model (CCSM4) and the Hadley Center Global Environment Model (HadGEM2-ES), and two CMIP6 models, the Community Earth System Model (CESM2-WACCM) and the UK Earth System Model (UKESM1-0-LL). These four models were all used for ISMIP6 Antarctica 2100 and were chosen based on the availability of extended simulations until 2300. One difference is that the extended CESM2 simulation uses the Whole Atmosphere Community Climate Model (WACCM) atmospheric core that resolves the upper atmosphere, while the CESM2 version used in ISMIP6 Antarctica 2100 was run with the Community Atmosphere Model (CAM6) atmospheric core that does not resolve the upper atmosphere (Danabasoglu et al., 2020). The non-atmospheric CESM2 components are identical in the two cases. Three of these four models (HadGEM2-ES, CESM2-WACCM, and UKESM1-0-LL) have equilibrium climate sensitivities, the change in the temperature in response to a doubling in the atmospheric carbon dioxide concentration, at the upper end of the 90% confidence interval estimated in the IPCC-AR6 (Meehl et al., 2020).

These four models form the base of most experiments (Table 1). Two additional experiments are based on the Norwegian Earth System Model (NorESM1-M) with repeated forcing until 2300 to allow comparison with previous results, since this model was a reference in the ISMIP6 Antarctica 2100 ensemble.

Table 1 List of ISMIP6 Antarctica 2300 Experiments in Tiers 1 and 2, With Forcing Derived From Global Climate Model Simulations From CMIP5 and CMIP6

Experiment	Global climate model	Scenario	Forcing	Ice shelf collapse	Tier
Historical	None	None	Free	No	Tier 1
ctrlAE	None	None	Free	No	Tier 1
expAE01	NorESM1-M	RCP2.6	Repeat	No	Tier 1
expAE02	CCSM4	RCP8.5	To 2300	No	Tier 1
expAE03	HadGEM2	RCP8.5	To 2300	No	Tier 1
expAE04	CESM2	SSP5-85	To 2300	No	Tier 1
expAE05	UKESM	SSP5-85	To 2300	No	Tier 1
expAE06	UKESM	SSP5-85	Repeat	No	Tier 1
expAE07	NorESM1-M	RCP8.5	Repeat	No	Tier 2
expAE08	HadGEM2	RCP8.5	Repeat	No	Tier 2
expAE09	CESM2	SSP5-85	Repeat	No	Tier 2
expAE10	UKESM	SSP1-26	To 2300	No	Tier 2
expAE11	CCSM4	RCP8.5	To 2300	Yes	Tier 2
expAE12	HadGEM2	RCP8.5	To 2300	Yes	Tier 2
expAE13	CESM2	SSP5-85	To 2300	Yes	Tier 2
expAE14	UKESM	SSP5-85	To 2300	Yes	Tier 2

Atmospheric Forcing

The derivation of the atmospheric forcing is similar to the approach described in Nowicki et al. (2020). It consists of forcings providing the surface mass balance and temperature anomalies at the surface of the Antarctic Ice Sheet compared to the reference period of 1995–2014. These files provide annual anomalies for the period 2015–2300. The surface mass balance anomalies are based on the precipitation, runoff, evaporation, and sublimation calculated by the global climate models and are provided as water-equivalent quantities. These anomalies are added to the reference surface mass balance used during model initialization, which varies among the ice sheet models. All these global climate model simulations were run with fixed ice sheet topography, although the actual Antarctic surface mass balance is influenced by changing surface elevation over time (Weertman, 1961). The inclusion of a surface-elevation feedback to correct the surface mass balance in response to changing ice sheet geometry is left to the discretion of modeling groups.

Oceanic Forcing

The oceanic forcing is derived from global climate model outputs as described in Jourdain et al. (2020) and used in ISMIP6 Antarctica 2100 (Seroussi et al., 2020). The ocean fields, including temperature and salinity, are first extrapolated into ice shelf cavities and other areas outside the ocean domains in the global climate models, where the ocean could advance during the simulations. Annual-mean forcing files of ocean temperature, salinity, and thermal forcing are provided for 2015–2300.

Ice flow models typically rely on sub-shelf melt parameterizations to convert ocean thermal forcing to basal melt. The choice and calibration of melt parameterization are left to the discretion of the modelers: they can be based on the parameterizations proposed for ISMIP6 (Jourdain et al., 2020) or any other parameterizations (DeConto & Pollard, 2016; Lazeroms et al., 2018; Martin et al., 2011; Pelle et al., 2019; Reese et al., 2018). The only constraint is that these parameterizations must use the ocean conditions provided in ice shelf cavities. Unlike Seroussi et al. (2020), we make no distinction between the ISMIP6 parameterizations and other approaches (i.e., no “standard” vs. “open” parameterizations). The choice of sub-shelf melting parameterization is treated as one of many modeling decisions left to the discretion of the individual modeling group, similar to calving and sliding laws.

Ice Shelf Collapse Forcing

Several experiments include ice shelf collapse due to hydrofracture (Pollard et al., 2015; Trusel et al., 2015). In order to parameterize this process, regions with more than 72.5 mm of annual liquid water for 10 years or more are considered prone to hydrofracture and therefore likely to collapse. This threshold is based on observations of recent ice shelf collapse in the Antarctic Peninsula and simulated conditions from regional climate models in this area at the time of their collapse (Trusel et al., 2015). It is identical to the condition applied in ISMIP6 Antarctica 2100 (Nowicki et al., 2020). Annual collapse masks are provided to specify the maximum extent of ice shelves at a given time, based on the amount of liquid precipitation at the ice sheet surface. The amount of liquid water is calculated based on the surface air temperature simulated by global climate models, using the nonlinear relationship between surface melting and summer air temperature derived in Trusel et al. (2015). This parameterization does not take into account the impact of water retention in firn (Donat-Magnin et al., 2021; van Wessem et al., 2023) or the required mechanical preconditioning (Lai et al., 2020) but can be easily included in a large ensemble of ice sheet models (Nowicki et al., 2020).

List of Experiments

Table 1 lists the experiments included in ISMIP6 Antarctica 2300. There is a historical run, a control run, and 14 projection experiments. The historical run, historical, extends from the model initialization date (typically in the late twentieth or early twenty-first century) to the simulation start date of January 2015. Some groups initialize their simulation to the beginning of 2015 and therefore do not submit a historical run. The other experiments, including the control, run from January 2015 until the end of 2300. In the control experiment, ctrlAE, climate conditions are unchanged and remain similar to the 1995–2014 period, with the exact surface mass balance and ocean conditions left to the discretion of modelers.

The 14 projection experiments, expAE01–expAE14, are divided into Tier 1, with six core mandatory experiments, and Tier 2, with eight additional experiments (see Table 1). These experiments are based on different global climate model simulations from CMIP5 and CMIP6 as described above. Most use a high-emission scenario (RCP8.5 or SSP5-8.5), while two experiments use a low-emission scenario (RCP2.6 or SSP1-2.6) for comparison. The majority of experiments are forced with atmospheric and oceanic conditions simulated by global climate models until 2300, but several experiments apply repeated forcing from years 2080–2100 for the final two centuries to compare the results with previous studies more easily. In the latter case, years are selected randomly from 2080 to 2100 to avoid repeating the same 20-year patterns, and each ice-flow model uses the same randomized forcing. We refer to these experiments as “repeat-forcing” experiments, as opposed to the “2300-forcing” experiments with global climate model forcing for the full period. Four experiments include ice shelf collapse as an additional forcing.

In summary, the projection experiments include:

• eight 2300-forcing experiments based on high-emission scenarios: four with ice-shelf collapse and four without collapse;

• four repeat-forcing experiments based on high-emission scenarios;

• two experiments based on low-emission scenarios, one with repeat forcing and one with 2300 forcing.

Ice Flow Models

Model Set-Up

Similar to the philosophy adopted for initMIP-Antarctica (Seroussi et al., 2019) and ISMIP6 Antarctica 2100 (Seroussi et al., 2020), there are no constraints on the methods or data sets used to initialize ice sheet models or the way physical processes are parameterized. The initialization date is also left to the discretion of each group to accommodate different modeling and initialization approaches, and allow the use of data sets registered to different observational periods. The only requirements are the ability to simulate ice shelves and grounding-line evolution (regardless of the numerical scheme used to do so) and to apply the climate forcing provided. The resulting ensemble includes a variety of model resolutions, stress-balance approximations, initialization methods, sliding and calving laws, and sub-ice-shelf melt parameterizations, representing the diversity of current ice sheet models.

Participating Models

Table 2 lists the 16 modeling groups and ice sheet modelers who submitted simulations to ISMIP6 Antarctica 2300. The ensemble includes 12 different ice flow models, a range of model resolutions, and various stress-balance approximations, basal sliding laws, and sub-ice-shelf melt parameterizations.

Table 2 Contributors, Modeling Groups, and Ice Flow Models for ISMIP6 Antarctica 2300 Projections

Contributors	Group	Ice flow model	Institutions
Jake Twarog	DC	ISSM	Dartmouth College, Hanover, NH, USA
Hélène Seroussi
Holly Kyeore Han	DOE	MALI	U.S. Department of Energy, Los Alamos National Laboratory, Los Alamos, NM, USA
Trevor Hillebrand
Matthew Hoffman
Justine Caillet	IGE	Elmer/Ice	Institut des Géosciences de l'Environnement, Grenoble, France
Fabien Gillet-Chaulet
Pierre Mathiot
Benoit Urruty
Nicolas Jourdain
Ralf Greve	ILTS	SICOPOLIS	Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan
Constantijn Berends	IMAU	UFEMISM	Institute for Marine and Atmospheric research Utrecht, Utrecht University, The Netherlands
Jorge Bernales
Roderik van de Wal
Christophe Dumas	LSCE	Grisli	Laboratoire des Sciences du Climat et de l'Environnement, Université Paris-Saclay, France
Aurélien Quiquet
Gunter Leguy	NCAR	CISM	NSF National Center for Atmospheric Research, Boulder, CO, USA
William Lipscomb
Heiko Goelzer	NORCE	CISM	Norwegian Research Center, Bergen, Norway
Petra Langebroek
David Chandler
Gunter Leguy			NSF National Center for Atmospheric Research, Boulder, CO, USA
William Lipscomb
Ann Kristin Klose	PIK	PISM	Potsdam Institute for Climate Impact Research, Germany
Julius Garbe
Torsten Albrecht
Ronja Reese			University of Northumbria, United Kingdom
Javier Blasco	UCM	Yelmo	Universidad Complutense de Madrid, Madrid, Spain
Alexander Robinson
Jorge Alvarez-Solas
Marisa Montoya
Tyler Pelle	UCSD	ISSM	University of California, San Diego, CA, USA
Violaine Coulon	ULB	Kori	Université libre de Bruxelles, Belgium
Frank Pattyn
Sainan Sun			University of Northumbria, United Kingdom
Sainan Sun	UNN	Úa	University of Northumbria, United Kingdom
Hilmar Gudmundsson
Chen Zhao	UTAS	Elmer/Ice	University of Tasmania, Australia
Yu Wang
Rupert Gladstone			Arctic Center, University of Lapland, Finland
Thomas Zwinger			CSC IT Center for Science, Espoo, Finland
Leopekka Saraste
Benjamin Galton-Fenzi			Australian Antarctic Division, Kingston, Australia
Jonas Van Breedam	VUB	AISMPALEO	Vrije Universiteit Brussel, Belgium
Philippe Huybrechts
Daniel Lowry	VUW	PISM	Antarctic Research Center, Victoria University of Wellington,
Nicholas Golledge	and GNS Science, New Zealand

Table 3 summarizes the main characteristics of the ice flow models, and Appendix C describes the models and initialization methods (including the differences among ensemble members for a given model) in more detail. Each group submitted between 1 and 11 sets of simulations, resulting in an ensemble of 43 sets of simulations. The models are initialized using a combination of long spin-ups, steady-state conditions, Data Assimilation (DA), and relaxation, with initialization years ranging from 1850 to 2015. The stress-balance approximations include combinations of the shallow Ice Approximation (SIA) and Shallow Shelf Approximation (SSA) (Bueler & Brown, 2009), as well as higher-order models (Blatter, 1995; Pattyn, 2003), and several depth-integrated models (Goldberg, 2011; Hindmarsh, 2004). The spatial resolution varies from 4 to 32 km for models with uniform grids. For models with spatially varying resolutions, the resolution can be as fine as 0.5 km near grounding lines and in shear margins, and as coarse as 200 km in the ice sheet interior. The sub-ice-shelf melt parameterizations include linear dependence on thermal forcing (Martin et al., 2011), plume models (Lazeroms et al., 2018), the Potsdam Ice-shelf Cavity mOdel (PICO) (Reese et al., 2018), the Potsdam Ice-shelf Cavity mOdel with plume (PICOP) (Pelle et al., 2019), and quadratic local and non-local ISMIP6 parameterizations (Jourdain et al., 2020). Various schemes were used to adjust melt rates near grounding lines (Leguy et al., 2021; Seroussi & Morlighem, 2018). Ice-front evolution was based on a minimum thickness condition, a retreat-only condition, strain rate (Levermann et al., 2012), or von Mises stress (Morlighem et al., 2016). Several models kept their ice front fixed except when prescribed by ice shelf collapse.

Table 3 Characteristics of Ice Flow Models Contributing to ISMIP6 Antarctica 2300

Ice flow model	Numerics	Stress balance	Resolution (km)	Initialization method	Initial year	Melt in partially floating cells	Ice front	Ice shelf melt parameterization	Bedrock adjustment
DC_ISSM	FE	SSA	2–50	DA	2015	Sub-Grid	MH	Quad. Non-local	No
DOE_MALI_4km*	FE/FV	HO	4–20	DA+	2000	Floating condition	RO	Quad. Non-local	No
DOE_MALI_8km_Ant95	FE/FV	HO	8–30	DA+	2000	Floating condition	RO	Quad. Non-local	No
DOE_MALI_8km_AntMean	FE/FV	HO	8–30	DA+	2000	Floating condition	RO	Quad. Non-local	No
IGE_ElmerIce	FE	SSA	1–50	DA	1995	No	Fix	PICO	No
ILTS_SICOPOLIS	FD	Hybrid	8	SP+	1990	Floating condition	MH	Quad. Non-local	ELRA
IMAU_UFEMISM1*	FD/FV	Hybrid	30–200	SP+	2014	Floating condition	Fix	Quad. Local	No
IMAU_UFEMISM2	FD/FV	Hybrid	30–200	SP+	2014	Floating condition	Fix	Quad. Local	No
IMAU_UFEMISM3	FD/FV	Hybrid	30–200	SP+	2014	Floating condition	Fix	Quad. Local	No
IMAU_UFEMISM4	FD/FV	Hybrid	30–200	SP+	2014	Floating condition	Fix	Quad. Local	No
LSCE_GRISLI*	FD	Hybrid	16	SP+	1995	Floating condition	MH	Quad. Non-local	No
LSCE_GRISLI2	FD	Hybrid	16	SP+	1995	No	MH	Quad. Non-local	No
NCAR_CISM1	FE/FV	L1L2	4	SP+	1995	Sub-Grid	RO	Quad. Non-local	No
NCAR_CISM2*	FE/FV	L1L2	4	SP+	2015	Sub-Grid	RO	Quad. Local	No
NORCE_CISM2	FE/FV	L1L2	8	SP+	1995	Floating condition	RO	Quad. Non-local Slope	No
NORCE_CISM3	FE/FV	L1L2	8	SP+	1995	Floating condition	RO	Quad. Non-local Slope	No
NORCE_CISM3_nonlocal*	FE/FV	L1L2	8	SP+	1995	Floating condition	RO	Quad. Non-local	No
NORCE_CISM3_local	FE/FV	L1L2	8	SP+	1995	Floating condition	RO	Quad. Local	No
NORCE_CISM4	FE/FV	L1L2	16	SP+	1995	Floating condition	RO	Quad. Non-local Slope	No
NORCE_CISM4_nonlocal	FE/FV	L1L2	16	SP+	1995	Floating condition	RO	Quad. Non-local	No
NORCE_CISM4_Local	FE/FV	L1L2	16	SP+	1995	Floating condition	RO	Quad. Local	No
NORCE_CISM4_JRA	FE/FV	L1L2	16	SP+	1995	Floating condition	RO	Quad. Non-local Slope	No
NORCE_CISM5	FE/FV	L1L2	32	SP+	1995	Floating condition	RO	Quad. Non-local Slope	No
NORCE_CISM5_nonlocal	FE/FV	L1L2	32	SP+	1995	Floating condition	RO	Quad. Non-local	No
NORCE_CISM5_local	FE/FV	L1L2	32	SP+	1995	Floating condition	RO	Quad. Local	No
PIK_PISM	FD	Hybrid	8	SP	1850	Floating condition	StR	PICO	No
UCM_Yelmo	FD	L1L2	16	SP+	1990	Sub-Grid	VM	Quad. Non-local	ELRA
UCSD_ISSM	FE	SSA	3–50	DA	2007	Sub-Grid	Fix	PICOP	No
ULB_Kori1*	FD	Hybrid	16	DA*	1950	No	Div	PICO	ELRA
ULB_Kori2	FD	Hybrid	16	DA*	1950	No	Div	Quad. Non-local	ELRA
UNN_Úa	FE	SSA	1–40	DA	2000	No	RO	Quad. Local	No
UTAS_ElmerIce	FE	SSA	1–25	DA	1995	Sub-Grid	Fix	Quad. Local	No
VUB_AISMPALEO	FD	SIA + SSA	20	SP	2000	N/A	MH	Quad. Non-local	No
VUW_PISM1*	FD	Hybrid	8	SP	2015	No	StR	Lin	VE
VUW_PISM1_s1	FD	Hybrid	8	SP	2015	No	StR	Lin	VE
VUW_PISM1_s2	FD	Hybrid	8	SP	2015	No	StR	Lin	VE
VUW_PISM1_s3	FD	Hybrid	8	SP	2015	No	StR	Lin	VE
VUW_PISM1_s4	FD	Hybrid	8	SP	2015	No	StR	Lin	VE
VUW_PISM2	FD	Hybrid	8	SP	2015	Yes	StR	Lin	VE
VUW_PISM2_s1	FD	Hybrid	8	SP	2015	Yes	StR	Lin	VE
VUW_PISM2_s2	FD	Hybrid	8	SP	2015	Yes	StR	Lin	VE
VUW_PISM2_s3	FD	Hybrid	8	SP	2015	Yes	StR	Lin	VE
VUW_PISM2_s4	FD	Hybrid	8	SP	2015	Yes	StR	Lin	VE

Unlike ISMIP6 Antarctica 2100, in which all simulations held the bedrock topography and ocean bathymetry fixed, several models adjusted the bathymetry and bedrock in response to the evolving ice load. Bedrock uplift can slow the retreat and mass loss of Antarctic ice streams and is important over a range of timescales (Adhikari et al., 2014; Barletta et al., 2018; Gomez et al., 2010; Han et al., 2022; Larour et al., 2019). It also has implications on how the sea-level contributions are calculated (Adhikari et al., 2020; Goelzer, Coulon, et al., 2020). This adjustment was based on either a viscoelastic deformable Earth model (Bueler et al., 2007; Lingle & Clark, 1985) or an elastic-lithosphere-relaxing-asthenosphere model (Le Meur & Huybrechts, 1996).

Model Outputs and Processing

Simulation results were submitted on a regular grid with a resolution of 4, 8, 16, or 32 km to be close to the models' native resolution. Outputs included annual values of two-dimensional fields and scalar quantities that capture the geometric evolution, ice flow, and forcings, similar to previous ISMIP6 Antarctic efforts (Seroussi et al., 2019, 2020). Scalar quantities (such as ice volume, ice volume above floatation, and ocean-induced melt) capturing the overall ice-sheet evolution were reprocessed based on the two-dimensional outputs for consistency between models. Scalar values were also calculated separately for the West Antarctic Ice Sheet, East Antarctic Ice Sheet, and Antarctic Peninsula using the Ice sheet Mass Balance Inter-comparison Exercise 2 (IMBIE2) basins (Shepherd et al., 2018).

Results

The results presented below are based directly on outputs from each experiment unless noted otherwise. Unlike what was done in previous ISMIP6 publications (Goelzer, Nowicki, et al., 2020; Payne et al., 2021; Seroussi et al., 2020, 2023), we do not subtract results of the control run. Since the forcing and changes until 2300 are large, the trend in the control run contributes only a relatively small fraction of the changes, and we are focusing mostly on substantial evolutions.

Many ice sheet modeling groups submitted several sets of experiments (e.g., 11 sets of experiments from NORCE) with different model settings, in which case they were asked to identify one submission as their “primary submission.” This decision was made by the modeling groups submitting several sets of experiments, based on the differences between these submissions, the ability of the submissions to capture recent changes, and their expertise. To assess the full range of uncertainty from the ensemble without giving too much weight to one model, we present results in two forms: (a) using the primary submission from each group to provide an equal weight to each modeling group and (b) including all the simulations submitted to capture the full range of results. Table 3 higlights the main submision from each group and Table 4 lists all the simulations submitted.

Table 4 Experiments Performed by Each Modeling Group for ISMIP6 Antarctica 2300 Projections

Experiment	DC_ISSM	DOE_MALI_4km	DOE_MALI_8km_Ant95	DOE_MALI_8km_AntMean	IGE_ElmerIce	ILTS_SICOPOLIS	IMAU_UFEMISM	LSCE_GRISLI	NCAR_CISM1	NCAR_CISM2	NORCE_CISM	PIK_PISM	UCM_Yelmo	UCSD_ISSM	ULB_Kori	UNN_Úa	UTAS_ElmerIce	VUB_AISMPALEO	VUW_PISM	Total number of experiments
Historical		X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X		32
ctrlAE	X	X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X	10	43
expAE01	X	X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X	10	43
expAE02	X	X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X	10	43
expAE03	X	X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X	10	43
expAE04	X	X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X	10	43
expAE05	X	X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X	10	43
expAE06	X	X	X	X	X	X	4	2	X	X	11	X	X	X	2	X	X	X	10	43
expAE07	X	X		X		X	4	2		X		X	X	X	2	X	X	X	10	29
expAE08	X	X		X		X	4	2		X		X	X	X	2	X	X	X		19
expAE09	X	X		X		X	4	2		X		X	X	X	2	X	X	X		19
expAE10	X	X		X		X	4	2		X		X	X	X	2	X	X	X	10	29
expAE11		X		X		X	4	2					X	X	2	X				14
expAE12		X		X		X	4	2					X	X	2	X				14
expAE13		X		X		X	4	2					X	X	2	X				14
expAE14		X		X		X	4	2					X	X	2	X				14

Historical Period and Control Run

The historical period covers the time between the model initialization date and the start of the experiments in January 2015. Since models have different initialization dates, the length of the historical period varies from 0 years (for UTAS_ElmerIce and VUW_PISM) to 175 years (for PIK_PISM). During the historical period, the mass change varies between −880 and 303 Gt/year, respectively for UNN_Úa and LSCE_GRISLI, equivalent to a sea-level contribution between 2.4 and −0.8 mm/year (Figure 1a). The IMBIE2 (Otosaka et al., 2023) estimates an Antarctic contribution to sea level increasing from 19 ± 39 Gt/year in 1997–2001 to 115 ± 55 Gt/year in 2017–2020. All models show a relatively linear and monotonic change. The surface mass balance (Figure 1b) is very similar for most simulations, as it is generally prescribed from regional climate models (Agosta et al., 2019; Mottram et al., 2021; van Wessem et al., 2018). It shows a large interannual variability in most cases, except for IMAU_UFEMISM and DOE_MALI, which used a constant surface mass balance over this period. Unlike the surface mass balance, the sub-ice-shelf melt applied to the different models varies widely among the historical simulations, from 10 to 3,080 Gt/year of melt overall (Figure 1c); observational estimates of sub-ice-shelf melt vary between 700 and 1,500 Gt/year (Paolo et al., 2023; Rignot et al., 2013). The VUW_PISM simulation has large interannual variations in both melting and refreezing, since basal melt is adjusted to reproduce the overall Antarctic mass loss.

[IMAGE OMITTED. SEE PDF]

Figure 2 shows how well the models capture observed ice sheet conditions at the start date of January 2015. We compare the initial ice thickness and velocity to observations from the BedMachine v2 and MEaSUREs Antarctica data sets interpolated onto the same regular grid (Morlighem et al., 2020; Rignot et al., 2013). The results show a root mean square error (RMSE) for ice thickness varying from 72 to 376 m. Models initialized either with DA of present-day conditions, or over longer periods but with the present-day thickness as a target, have the smallest errors. The RMSE for the ice velocity varies between 26 and 258 m/year. Models that assimilate present-day conditions have smaller errors, as expected. This ability to capture present-day conditions has to be balanced with the ability to capture recent changes, as representing both accurate ice sheet conditions and trends at a given time demonstrates the reliability of ice sheet simulations.

[IMAGE OMITTED. SEE PDF]

In the control run, ice sheets evolve under constant climate forcing. Figure 1a shows the evolution between 2015 and 2300, for the main submissions. By 2300, the ice mass change varies from −7.7 · 10⁴ to 15.9 · 10⁴ Gt, equivalent to a range from 210 mm of sea level rise to 440 mm of sea-level drop. Most trends are linear, since the applied forcing is constant and there is little interannual variability in surface mass balance or ice-shelf melt (Figures 1b and 1c). There are, however, large differences among models: the surface mass balance varies from 2050 to 3323 Gt/year and the ice-shelf melt from 175 to 2040 Gt/year. The very high total surface mass balance value of UCM_Yelmo comes from the larger ice shelves than present day in the control run; when the observed present-day ice extent is used for this calculation, the total surface mass balance for this model is much closer to the others with a value of 2450 Gt/year. The models initialized with present-day DA have the largest ice mass changes; DC_ISSM has the largest loss and UNN_Úa the largest gain. Within an ensemble from the same model, the mass changes can be either clustered or far apart, depending on the parameters varying between the simulations.

Evolution Until 2300 Under High-Emission Scenarios

Under high-emission scenarios, surface mass balance forcing from the four climate models is applied in a similar way by all the ice sheet models, with relatively limited differences among ice flow models and similar interannual variability patterns captured by all the ice models (Figure 3a). The initial surface mass balance in 2015 varies from ∼2,000 to 3,600 Gt/year, depending on the ice extent and the choice of reference surface mass balance. It increases slightly until 2100 for all climate models, after which it stays relatively constant for CCSM4 but decreases for the three other climate models. CESM2 has the largest surface mass balance decrease, driving a negative overall surface mass balance for some ice models starting around 2200. CESM2 also has the largest surface mass balance spread in 2300, from ∼−600 to 3,400 Gt/year.

[IMAGE OMITTED. SEE PDF]

Figure 3b shows the total ice-shelf basal melt under evolving ocean forcing. In contrast to the surface mass balance, total melt varies significantly among ice models for a given ocean forcing. The total melt increases overall until 2100, after which the trend depends on the ice model and its response to the forcing. For some models the melt continues to increase, while for others the initially higher melt rates decrease as ice shelves melt and thin, and fronts retreat. Starting around 2150, the ice shelves in several simulations have entirely melted. There is no clear pattern for simulations forced with the same climate model; the differences come primarily from ice model differences such as the ice-shelf geometry and the melt parameterization (Table 3). Ice shelf melting removes mass from the ice sheet, but has a very limited direct effect on sea level, since the ice shelves are already floating. However, these changes raise sea level indirectly by reducing buttressing for the ice streams that feed them (Reese et al., 2017).

Figure 4 shows the Antarctic Ice Sheet mass change for the four Tier 1 experiments with high-emission scenarios extended to 2300 (expAE02–expAE05, one for each climate model). Mass loss is limited until 2100 but increases rapidly thereafter. The overall change in ice mass above floatation between 2015 and 2300 ranges from −16 · 10⁵ to 2 · 10⁵ Gt, the equivalent of 4.4 m to −0.6 m of sea level rise. Both the climate and ice-flow models introduce large uncertainties. Forcing with CCSM4 leads to the lowest mass loss, between −0.6 and 1.9 m SLE. The other three climate models lead to relatively similar mass loss, between −0.4 and 4.4 m SLE by 2300. Considering just the main submission for each ice-flow model, IMAU_UFEMISM has the lowest spread (0.9 m) for the four experiments, from −0.44 to 0.45 m SLE, while ILTS_SICOPOLIS has the largest spread (2.5 m), from 1.9 to 4.4 m SLE. When all simulations are included, models with several submissions show a larger spread: the 11 NORCE submissions have a range of 3.3 m SLE, as do the 10 VUW submissions.

[IMAGE OMITTED. SEE PDF]

Figure 5 shows the changes in modeled sea level contributions at years 2100, 2200, and 2300 under high emission climate forcings for the main three Antarctic regions: the West Antarctic Ice Sheet, the East Antarctic Ice Sheet, and the Antarctic Peninsula. The West Antarctic Ice Sheet is the largest sea-level contributor across all climate models and time periods. By 2100, sea level changes are limited, with less than 0.1 m SLE for all regions; HadGEM2 drives the largest mass loss for both West and East Antarctic ice sheets. The West Antarctic Ice Sheet loses mass for all climate models, while the East Antarctic Ice Sheet gains mass or experiences minimal changes. Mass loss increases rapidly by 2200, reaching on average 0.2–0.5 m SLE for West Antarctica, while the average East Antarctic changes vary between −0.2 and 0.2 m SLE, and the Antarctic Peninsula contribution remains small. The spread associated with the different ice flow models is similar for the West and East Antarctic Ice Sheets. The sea-level contribution continues to increase rapidly until 2300, when the West Antarctic Ice Sheet average is 0.6–1.3 m SLE, and much of the ice in this region has been removed. For East Antarctica, CCSM4 continues to show small changes and a large spread, while the other three models show a contribution of 0.3–0.6 m SLE with a lower spread than in West Antarctica. The contribution from the Antarctic Peninsula again is smaller; HadGEM2 and CESM2 have an average mass loss of 0.09 and 0.14 m SLE, respectively. When all the ice simulations are included (Figure B1), the average contributions are similar for all regions and times, but the spread increases significantly, especially for the East Antarctic contribution in 2100.

[IMAGE OMITTED. SEE PDF]

Figure 6 shows spatial patterns of the Antarctic Ice Sheet thinning (red shading) and thickening (blue shading) for the ensemble mean, as well as the standard deviation of thinning, for the main model submissions under high-emission scenarios (expAE02–expAE05). By 2100, ice shelves thin up to ∼200 m, while grounded ice thinning of ∼50 m is confined to coastal sectors of the Thwaites, Cook, and Totten glaciers, with limited changes elsewhere. Thinning increases rapidly by 2200, with large portions of West Antarctica and coastal sectors of East Antarctica thinning by more than 500 m. By 2300, grounded ice has thinned across the entire West Antarctic Ice Sheet, with thinning of more than 1,000 m projected for Thwaites Glacier and for both grounded and floating ice in the Filchner-Ronne and Ross regions. There is also extensive thinning in the Wilkes and Aurora subglacial basins (housing Cook and Totten Glaciers, respectively) that extends ∼400 km into the East Antarctic interior. By 2300, ice thickens by up to 200 m in large parts of the East Antarctic interior. The intermodel spread in the ensemble is high in some regions; the standard deviation is comparable to the total thinning signal for Pine Island and Thwaites glaciers in 2300. There is greater agreement for the Bungenstock Ice Rise (upstream of the Ronne Ice Shelf) and Siple Coast (upstream of the Ross Ice Shelf), where thinning exceeds 1,000 m and the standard deviation of the ensemble is <200 m. The standard deviation over the large ice shelves is low in 2300, when most ice shelves have completely melted. Figure B2 shows results similar to Figure 6 but including all the submissions.

[IMAGE OMITTED. SEE PDF]

Extensive Regional Retreat

Figure 7 shows the percentage of the ensemble that projects retreat of currently grounded ice for the main submissions under high-emission scenario (expAE02–expAE05) in 2100, 2200, and 2300. In 2100, retreat is mostly confined to the Bungenstock Ice Rise and Siple Coast, where up to 60% of the ensemble ungrounds. By 2200, ensemble retreat in these regions is more widespread and extends farther into the West Antarctic interior, with 90%–100% of the ensemble ungrounding near the Bungenstock Ice rise and 80%–90% along the Siple Coast. Nearly all ensemble members simulate ungrounding of the Korff and Henry Ice Rises, which currently stabilize the Ronne Ice Shelf. Less than 40% of the ensemble projects retreat of Thwaites and Pine Island glaciers at this time. By 2300, however, this number increases to 80%, with 20%–30% simulating a collapse of Pine Island and Thwaites glaciers. Furthermore, 90%–100% of the ensemble projects Siple Coast retreat up to ∼200 km upstream from the present-day grounding line. Along coastal sectors of the East Antarctic Ice Sheet, 60%–70% of the ensemble projects extensive retreat of Moscow University Glacier, Totten Glacier, and the glaciers draining into Queen Maud Land.

[IMAGE OMITTED. SEE PDF]

To further illustrate these spatial and temporal retreat patterns, Figure 8 shows the retreat for two high-emission scenario experiments (expAE03 and expAE04) for each model's main submission along three sections (red lines in Figure 7d) taken through Thwaites Glacier, the Bungenstock Ice Rise, and the Siple Coast. In each panel, the x-axis denotes the distance from the present-day grounding line (Morlighem et al., 2020), the y-axis denotes the simulation year (with time increasing upward), and each colored marker shows the position of the grounding line along the flowline for each year of the simulation. That is, for a particular ensemble member, grounding-line retreat appears as colored markers shifting to the right in each panel, while grounding line stabilization appears as colored markers stacked vertically.

[IMAGE OMITTED. SEE PDF]

For Thwaites Glacier (panels a and d), about half of the models remain grounded near the present-day grounding line through 2300. For other models, retreat ensues primarily after 2100 and 2200 in these two high-emission scenario experiments (expAE03 and expAE04, respectively), with a consistent retreat of ∼3.5 km/year across these models. The grounding line stabilizes on topographic highs and sections of prograde bed topography (bed topography that slopes upward toward the ice sheet interior; see the brown shading at the bottom of Figure 8), with retreat sometimes extending ∼430 km upstream of the present-day grounding line, to a topographic high at 700-m depth.

Along the Bungenstock Ice Rise, rapid retreat across smooth bed topography at the start of both experiments is modeled in nearly all simulations. Retreat accelerates down the retrograde bed topography by more than 10 km/year during 2100–2150 and 2050–2100 in expAE03 and expAE04, respectively. The grounding line stabilizes about 300 km upstream on prograde bed topography in all simulations. Siple Coast retreat is more variable; the retreat starts by 2150, progresses at rates of 0.8–1.5 km/year, and extends up to 420 km into the West Antarctic interior. Unlike the retreat of Thwaites and Bungenstock, there is limited correspondence between retreat rates and the bed topography, likely due to retreat and stabilization occurring in areas around the transect. Including all the simulations leads to similar results (see Figure B3). Overall, these results highlight the heterogeneous nature of retreat in different Antarctic basins and the importance of bed topography in controlling rapid and extensive retreat.

Impact of Ice Shelf Collapse

Observations and modeling studies have shown the importance of ice-shelf collapse on the Antarctic Ice Sheet evolution (Hulbe et al., 2008; Khazendar et al., 2015; Pollard et al., 2015; Scambos et al., 2004; Schannwell et al., 2018; Sun et al., 2020). We therefore included four experiments with ice-shelf collapse prescribed based on the presence of liquid water at the ice surface to test the influence of reduced buttressing from ice shelves (expAE11–expAE14, see Table 1). Figure 9 shows the timing and spatial evolution of shelf collapse for each climate model. The timing and spatial patterns vary significantly among the four models, but most regions with floating ice experience some collapse by 2300. The collapse starts in the Antarctic Peninsula, where most ice shelves are affected by 2100, and rapidly spreads to other regions after 2100. For the large Filchner–Ronne, Ross, and Amery ice shelves, collapse starts around 2150 and propagates inland, except for HadGEM2 where the collapse of these shelves is limited. We can observe the contrast between the CMIP5 and CMIP6 models. The two CMIP5 models, CCSM4 and HadGEM2, have vastly different spatial and temporal patterns of collapse for these large ice shelves: limited for HadGEM2 and extensive for CCSM4. For the CMIP6 models, CESM2 and UKESM, the patterns are similar to each other, with collapse starting in the grounding zone and proceeding inland, leading to the removal of large ice shelves once the collapse reaches the entire grounding zone area.

[IMAGE OMITTED. SEE PDF]

Both CESM2 and UKESM have amplified melt leading to ice shelf collapse across the ice sheet escarpment zones, where downslope katabatic and foehn winds are known to drive localized warming and positive wind–melt–albedo feedbacks today (Lenaerts et al., 2017). This is particularly apparent across the Siple Coast ice streams and Transantarctic Mountain glaciers feeding into the Ross Ice Shelf, as well as across the innermost Filchner–Ronne Ice Shelf in Figure 9. Similar locally enhanced melt in grounding zones occurs more broadly in the CMIP6 models, particularly for the Larsen C, Amery, and Queen Maud Land ice shelves, though this is masked in Figure 9 by the wider areas of high melt over these shelves. These patterns mimic observed melt and lake patterns today (Bell et al., 2018; Moussavi et al., 2020; Trusel et al., 2013), suggesting that these areas will be the sites of high future melt and supraglacial lake formation. The strong amplification of melt in grounding zones in the CMIP6 models implicates adiabatic warming of downslope winds (and the more sensitive treatment of these processes) as a leading driver of future melt and ice shelf collapse. In comparison, shelf collapse in the CMIP5 models appears to represent broader, transient atmospheric warming.

Figure 10 shows how ice shelf collapse affects the overall Antarctic mass loss for the main submission of each modeling group and Figure B4 for all the model submissions. Fewer ice flow models performed these experiments (Table 4), as they were in Tier 2 (Table 1), and some ice flow models are not able to prescribe ice-front evolution. When ice-shelf collapse is included, the mass loss can reach up to 6.9 m SLE, an additional 2.5 m compared to experiments without collapse. The additional mass loss is 1.1 m SLE on average, but varies between 0 and 6.3 m SLE.

[IMAGE OMITTED. SEE PDF]

The top row of Figure 11 shows the percentage of models (out of all models that ran the collapse experiments) that simulate ungrounding of currently grounded ice in years 2100, 2200, and 2300. The bottom row shows the percentage difference between corresponding experiments with and without ice-shelf collapse (i.e., expAE11–expAE14 minus expAE02–expAE05). By 2200 and 2300, shelf collapse results in widespread retreat of coastal glaciers farther upstream than in non-collapse experiments. In particular, by 2200, 40%–50% more models in the ensemble project retreat of Pine Island Glacier up to 50 km inland, while 20%–30% project greater inland retreat of Thwaites Glacier and the Siple Coast glaciers. By 2300, shelf collapse enhances the possibility of large-scale West Antarctic collapse by 30%–40% via widespread inland retreat of Pine Island and Thwaites glaciers, as well as the glaciers feeding the western Ronne ice shelf.

[IMAGE OMITTED. SEE PDF]

Repeat Forcing and 2300 Forcing

Most previous studies investigating the evolution of the Antarctic Ice Sheet beyond 2100 relied on repeat forcing from the end of the twenty-first century, as forcing from global climate models beyond 2100 is not widely available (Lipscomb et al., 2021; Lowry et al., 2021). In order to assess the impact of repeat forcing, we compare the previous high-emission simulations (expAE02–expAE05) with simulations that use repeated, randomly sampled late-twenty-first century forcing after 2100 (expAE06, expAE08, and expAE09; see Section 2.5 and Table 1 for details). Figure 12 shows the evolution of ice sheet mass loss for these three experiments (main model submissions only). Results are identical to the previous experiments until 2100 and diverge thereafter. Compared to the previous experiments, mass change is smaller and varies between −7 · 10⁵ Gt and 2 · 10⁵ Gt by 2300, or between 3 and −0.4 m SLE. The average sea-level contribution is 1.6 m smaller on average, with differences ranging from 0.2 m SLE to 4.4 m SLE, depending on the ice and climate models. The results (Figure B5) have a larger spread between 3 m and −0.4 m when all simulations are included.

[IMAGE OMITTED. SEE PDF]

Low-Emission Scenarios

Two experiments use forcing from low-emission scenarios. Experiment expAE01 uses climate forcing from NorESM under RCP2.6, with repeat forcing from the late twenty-first century after 2100, while expAE10 uses UKESM forcing under SSP1-2.6 extended to 2300. Figure 13 shows the relatively small mass changes in these two low-emission scenario experiments, with a mix of positive and negative sea-level contributions. Mass change in expAE01 varies between −1.7 · 10⁵ and 1.3 · 10⁵ Gt (or 460 and −355 mm SLE). Mass change in expAE10 is overall similar, varying between −2.3 · 10⁵ and 1.4 · 10⁵ Gt (or 622 and −374 mm SLE). The intermodel differences are small until 2100 but increase through 2300.

[IMAGE OMITTED. SEE PDF]

Discussion

Sea Level Contribution of the Antarctic Ice Sheet

The results presented above show the evolution of the Antarctic ice sheet under high-emission scenarios, in order to provide upper bounds to sea level and assess regions most susceptible to large changes. While the sea-level contribution from the Antarctic Ice Sheet is relatively limited during the twenty-first century (less than 30 cm SLE by 2100), it increases rapidly afterward under these high-emission scenarios, reaching up to 1.7 and 4.4 m SLE by 2200 and 2300, respectively (Figure 4). These upper bounds are comparable to previous studies investigating Antarctic evolution over several centuries (Bulthuis et al., 2019; Golledge et al., 2015; Lowry et al., 2021), but much lower than the values up to 8 and 14 m SLE by 2200 and 2300, respectively, reported in DeConto et al. (2021) when the Marine Ice Cliff Instability is included following ice shelf collapse. The lower bounds for the high-emission scenarios differ widely among the different studies, with a small mass loss in our study and in Bulthuis et al. (2019), a medium lower bound of ∼1.5 m SLE by 2300 in Golledge et al. (2015) and Lowry et al. (2021), and a much larger lower bound of 7 m by 2300 in DeConto et al. (2021). Projections from this ISMIP6 ensemble are comparable with the values reported in the IPCC AR6 for the assessed Antarctic ice-sheet contribution in 2300 of −0.28 to +3.13 m SLE for the likely range (IPCC, 2021) and with the ISMIP6 Antarctica 2100 results until 2100 (Seroussi et al., 2020).

The mass loss and retreat come mostly from the West Antarctic Ice Sheet (Figure 5). There is substantial thinning and grounding-line retreat in most West Antarctic basins by 2200, and an extended collapse of West Antarctica in some simulations by 2300. In East Antarctica, thinning is widespread and affects many glaciers in Wilkes Land and Queen Maud Land (Figure 6), but large grounding line retreat is mostly limited to Totten, Moscow University, Cook, and Ninnis glaciers (Figure 7). The timing of retreat onset varies widely among the different models, but the retreat rate, once initiated, is similar for most simulations, as shown for Thwaites Glacier, the Siple Coast, and the Bungenstock Ice Rise in Figure 8. For this large ensemble of simulations with a variety of ice flow models, the speed of grounding-line retreat in regions of retrograde bed is controlled mostly by the bedrock topography, as was suggested by previous studies based on a single ice flow model (Jones et al., 2021; Seroussi et al., 2017).

Adding ice-shelf collapse leads to faster retreat after 2100 in the West Antarctic Ice Sheet and a greater number of simulations showing its collapse by 2300, leading to an additional 1.1 m SLE on average. The difference between the low (Figure 13) and high (Figure 4) emission scenarios is relatively small in the twenty-first century but grows rapidly afterward. The experiments with low-emission forcing have less than 0.5 m SLE by 2300. These results from a large ensemble of models confirm those of Lowry et al. (2021), which suggested that the impact of emission scenario emerges around 2150 due to current limitations in our understanding of ice flow processes and their representation in ice flow models.

Sea-level rise is a main consequence of climate change, and sea-level projection is critical for risk assessment, adaptation planning, and policy analysis. Using results from the high-emission scenario experiments, including each group's main submission with and without ice shelf collapse (expAE02–expAE05 and expAE11–expAE14, respectively), Figure 14 shows the timing of 0.5 and 1 m SLE from the Antarctic Ice Sheet. The first simulations with 0.5 and 1 m SLE appear in 2120 and 2145, respectively, after which the number of simulations reaching these levels grows rapidly. Half of the simulations project 0.5 m SLE by 2160 when ice-shelf collapse is included; this happens 100 years later without shelf collapse. Half of the ensemble reaches 1 m SLE by 2225 with collapse, while only about 50% reach this threshold by 2300 without collapse in these high-emission scenarios.

[IMAGE OMITTED. SEE PDF]

Calculation of Sea-Level Contribution

We have calculated the Antarctic Ice Sheet contribution to global sea level based on the change in mass above floatation and assuming 1 mm SLE for 362.5 Gt of mass loss, similar to most studies during the past decade (Bindschadler et al., 2013; Levermann et al., 2020; Sun et al., 2020). Recent studies, however, have investigated several correction terms that are not included in the mass above floatation (Adhikari et al., 2020; Goelzer, Coulon, et al., 2020). The first correction is based on the density difference between fresh meltwater and saline ocean water. This density correction amounts to about 3% of the ice mass grounded below sea level that is lost (Goelzer, Coulon, et al., 2020). This additional contribution to sea level depends on the ice and ocean water density assumed in each model, as well as the ratio of total ice mass to ice mass above floatation that is changed. A second term, the bedrock correction, is related to changes in subglacial bedrock and ocean bathymetry due to the isostatic response of bedrock (i.e., uplift or subsidence) to variations in ice loading (Gomez et al., 2015; Larour et al., 2019; Pan et al., 2021). Its common function in both Adhikari et al. (2020) and Goelzer, Coulon, et al. (2020) is to counter changes in mass above floatation that occur at grounded portions of the ice sheet due to isostatic changes, but without sea-level contribution from actual ice mass changes.

How to best include these corrections in standalone ice sheet models remains a topic of active research, and the solutions proposed by Adhikari et al. (2020) and Goelzer, Coulon, et al. (2020) differ. Goelzer, Coulon, et al. (2020) include changes in bedrock elevation over the entire grid domain, to take into account the so-called water-expulsion effect (Pan et al., 2021), while Adhikari et al. (2020) include changes only within the evolving grounded ice sheet extent. Table 5 shows the effect of including density and bedrock corrections on the Antarctic Ice Sheet sea-level contribution in 2300 for one high-emission scenario experiment (expAE03). For models with a constant bedrock elevation, only the density correction is included. This correction is similar for both methods and increases the sea level contribution by ∼6%–29% overall, depending on the ratio of ice mass and ice mass above floatation lost during the simulations. The density correction captures the same process in both methods and should be identical for models that do not include bedrock adjustment. The small differences in sea level correction between these two methods in Table 5 come from the different integration of results over grid cells partially filled with ice between the two methods, as well as the integration intervals chosen for Adhikari et al. (2020).

Table 5 Sea Level Contribution of the Antarctic Ice Sheet Calculated Based on Volume Above Floatation (VAF), Sea Level Correction From Goelzer, Coulon, et al. (2020) (G2020) and Sea Level Correction From Adhikari et al. (2020) (A2020) for the Main Submissions of Experiment expAE03

Ice flow model	Bedrock adjustment	Volume above float. (VAF)	Sea level G2020	Sea level A2020
m SLE	m SLE (% diff. with VAF)	m SLE (% diff. with VAF)
DC_ISSM	No	0.98	1.07 (8.9%)	1.06 (8.4%)
DOE_MALI	No	2.90	3.07 (5.7%)	3.07 (5.7%)
IGE_ElmerIce	No	0.46	0.54 (18.6%)	0.54 (19.2%)
ILTS_SICOPOLIS	ELRA	4.25	4.56 (7.4%)	4.51 (6.2%)
IMAU_UFEMISM1	No	0.43	0.50 (16.7%)	0.48 (11.4%)
LSCE_GRISLI	No	2.60	2.75 (6.1%)	2.71 (4.4%)
NCAR_CISM	No	3.33	3.52 (5.9%)	3.53 (6.2%)
NORCE_CISM	No	1.00	1.09 (9.6%)	1.10 (10.2%)
PIK_PISM	No	1.98	2.10 (5.9%)	2.04 (3.3%)
UCM_Yelmo	ELRA	1.31	1.54 (17.6%)	1.45 (11.1%)
UCSD_ISSM	No	0.63	0.68 (8.1%)	0.68 (8.1%)
ULB_Kori	ELRA	1.68	1.80 (7.0%)	1.79 (6.6%)
UNN_Úa	No	1.01	1.14 (12.5%)	1.14 (12.5%)
UTAS_ElmerIce	No	0.31	0.40 (28.7%)	0.49 (26.8%)
VUB_AISMPALEO	No	0.60	0.67 (11.2%)	0.71 (17.8%)
VUW_PISM	VE	1.83	2.88 (57.6%)	2.35 (28.2%)

For models that include bedrock adjustment from an Elastic Lithosphere–Relaxed Asthenosphere (ELRA) or viscoelastic deformable Earth model, the adjustment depends on the correction method. Bedrock corrections for 2300 vary from 7% to 17% based on Goelzer, Coulon, et al. (2020) and from 6% to 12% based on Adhikari et al. (2020) for the three models that include an ELRA bedrock adjustment, while the correction reaches 58% and 28% for the bedrock adjustment based on the visco-elastic deformable Earth, as this method creates a very strong response to ice unloading that extends over the entire Southern Ocean. The Goelzer, Coulon, et al. (2020) correction is consistently larger since it integrates changes in bedrock uplift and therefore considers the water expulsion effect over the entire grid, including under floating ice shelves and in the open ocean, while the Adhikari et al. (2020) only considers the changes in the ice domain at a given time. The difference between the two methods continues increasing further over time and is larger for significant bedrock changes. While it is beyond the scope of this study to decide which method is the most appropriate conversion of mass loss into sea-level contribution, these differences highlight the importance of establishing a unified framework within the community for accurate sea-level corrections in multi-century simulations.

Sources of Uncertainty

Uncertainties in the projections presented here come from several sources, including the choice of ice model, climate model forcing, emission scenario, and the inclusion of processes like ice-shelf collapse. Analysis of the ISMIP6 Antarctica 2100 ensemble showed that uncertainty in the dynamic mass loss is dominated by the ice flow model, but the role of climate forcing increases over the twenty-first century (Seroussi et al., 2023). Here we perform a similar analysis to assess the role of different sources of uncertainty on mass loss.

We use analysis-of-variance theory (ANOVA, Girden 1992; von Storch & Zwiers, 1999) to decompose the variance of the ensemble into the contribution of the different components and the interactions among these components. This statistical technique allows one to compare and analyze the means of two or more groups to determine if statistically significant differences exist between these groups. The total variability in the data is broken down into two components: the variance within each group (within–group variance) and the variance between the groups (between–group variance). The significance of differences in group means is evaluated based on the ratio of these two variance components. This approach allows one to decompose the variance of the Antarctic sea level contribution into the contribution of different components and the interactions between these different components. We include here three sources of uncertainty: ice flow models, climate models, and the inclusion or exclusion of ice shelf collapse. The total uncertainty is therefore partitioned between contributions from climate models, ice flow models, the inclusion or exclusion of ice shelf collapse, the two-way interactions between these terms (referred to as ice–climate, ice-collapse and climate-collapse), and the three-way interactions between these terms (referred to as 3-way interactions). We consider only the high-emission experiments since very few experiments were performed with low and medium emission scenarios. Unlike previous studies where missing experiments were emulated using Gaussian process or neural networks (Seroussi et al., 2023; Van Katwyk et al., 2023), we include in our ensemble only the sets of simulations that include experiments with and without ice-shelf collapse (i.e., both expAE02–expAE05 and expAE11–expAE14). It is impossible to emulate ice shelf collapse for models without collapse, since its impact varies strongly from one model to the next (Section 4.4). The reduced ensemble includes eight ice flow models and four climate models running experiments with and without ice-shelf collapse, for a total of 64 experiments.

Figure 15 shows the evolution of the overall uncertainty and the different sources of uncertainty through 2300. The overall uncertainty increases over the entire period, reaching ∼1.6 m SLE by 2300, and is dominated by the choice of ice flow model over the entire period. The role of climate forcing increases over the twenty-first century to reach ∼15% of the variance by 2100 and remains relatively constant afterward. The role of shelf collapse is negligible until 2100 but increases rapidly thereafter, reaching ∼10% of the variance by 2150 and then remaining roughly constant. The interactions among the different terms represent a smaller part of the variance and sum to less than ∼25% over the entire period. These results show that uncertainties linked to differences in ice flow models dominate the projections even on timescales of several centuries.

[IMAGE OMITTED. SEE PDF]

Seven groups submitted ensembles of simulations, with 2–11 ensemble members per group. The spread of results within these small ensembles depends strongly on the number of simulations and the parameters varied within the ensemble. For example, the spread in the DOE and IMAU ensembles is low compared to the overall spread, but encompasses most of the overall spread in the VUW ensemble (Figure 4). The range of projections is overall similar whether we use the full ensembles or only the main submission for each group, since the upper and lower bounds are similar in each case (Figure 4). However, the inclusion of ensembles affects the uncertainty, especially when the mass loss is low, as for the East Antarctic Ice Sheet in 2100 (see Figures B1 and 5). Both NORCE and VUW simulate large mass loss compared to the other ice flow models but did not submit simulations with ice shelf collapse, so the large number of simulations in these ensembles affects the timing to reach 0.5 and 1 m SLE for simulations without ice shelf collapse. These thresholds are reached about 50 years earlier when including all simulations instead of just the main submissions (see Figures B6 and 14). We need new methods to take advantage of the entire ensemble of simulations without biasing the overall results, for example, based on emulators (Edwards et al., 2021; Van Katwyk et al., 2023).

These Antarctic Ice Sheet projections have limitations associated with both the ice flow models and the external forcing, and are mostly based on high-emission scenario conditions. The limitations of ice flow models, even after a decade of rapid progress, have been described in previous intercomparisons (Bindschadler et al., 2011; Nowicki et al., 2013a, 2013b) and reviews (Goelzer et al., 2018; Pattyn et al., 2018). The use of forcing from global climate models leads to further limitations (Nowicki & Seroussi, 2018); in particular, these models do not capture ocean circulation in ice-shelf cavities (Jourdain et al., 2020) and some contain significant bias in ocean temperatures (Barthel et al., 2020). The derivation of ice-shelf collapse based on surface liquid water is oversimplified in order to be included in as many models as possible (Nowicki et al., 2020). Finally, as these experiments use standalone ice sheet models, they do not capture feedbacks between ice sheets, the rest of the climate system, and the solid-Earth and sea-level. Important missing mechanisms include surface mass balance-elevation feedback (Coulon et al., 2023; Tewari et al., 2021), the melt-stratification feedback, whereby ice melting can warm the ocean near ice shelves by increasing upper ocean stratification and suppressing interaction between warm subsurface waters and the cold atmosphere (Golledge et al., 2019; Silvano et al., 2018), and the rotational and gravitational feedbacks that create negative feedbacks on marine based ice sheet (Coulon et al., 2021; Han et al., 2022). Since these processes are key for the long-term Antarctic evolution, efforts should continue to understand these processes and include them in ice flow and global climate models.

Conclusions

ISMIP6 Antarctica 2300 provides the first multi-century, multi-model projections of the Antarctic Ice Sheet evolution using an ensemble of 16 ice sheet models forced with climate conditions derived from global climate models. Our analysis highlights the large increase in Antarctic mass loss after 2100 under high-emission scenarios, with simulations suggesting sea level contributions of up to 0.28 m of SLE by 2100, but up to 1.7 and 4.4 m SLE by 2200 and 2300, respectively. Antarctic mass loss is dominated by a few key regions in West Antarctica, including the Bungenstock Ice Rise (upstream of the Ronne Ice Shelf), the Siple Coast (upstream of the Ross Ice Shelf), and the Amundsen Sea sector. Simulations of the Amundsen Sea sector show strong variations across ice flow models, with some models simulating collapse of this region between 2200 and 2300, while others suggest modest retreat. The intermodel variation is less for the Bungenstock and Siple Coast regions. There are consistent retreat patterns across models, with periods of slow retreat alternating with rapid retreat in a given region. During fast retreat phases, retreat rates are consistent across ice flow models, suggesting control by the bed topography rather than the details of model physics.

These simulations also highlight the critical role of the emission scenario. While mass losses from low- and high-emission scenarios are relatively similar in the twenty-first century, the difference between scenarios grows rapidly after 2100, underscoring the importance of emission reductions for the long-term Antarctic Ice Sheet stability and the need for CMIP scenarios extending beyond the twenty-first century. This study furthermore sheds lights on the importance of ice-shelf collapse, which is limited in the twenty-first century but can significantly increase the mass loss thereafter. Finally, the results from this ensemble underscore that the choice of ice sheet models continues to be a leading source of uncertainty, even on multi-century timescales, emphasizing the need for continued improvements in numerical ice sheet models and interactive ice sheet coupling in climate models.

Appendix

Appendix A - Table of Acronyms

This appendix provides a table of all the acronyms used in the text, tables, and figures.

Table A1 List of Acronyms Used in the Text, Tables and Figures

Acronym	Name
CAM	Community Atmosphere Model
CCSM	Community Climate System Model
CESM	Community Earth System Model
CISM	Community Ice Sheet Model
CMIP5	Coupled Model Intercomparison Project phase 5
CMIP6	Coupled Model Intercomparison Project phase 6
ELRA	Elastic Lithosphere–Relaxed Asthenosphere
HadGEM	Hadley center Global Environment Model
HO	Higher-Order approximation
IMBIE2	Ice sheet Mass Balance Inter-comparison Exercise 2
ISMIP6	Ice Sheet Model Intercomparison for CMIP6
ISSM	Ice-sheet and Sea-level System Model
L1L2	Depth-integrated stress balance approximation
MALI	Model for Prediction Across Scales
NorESM	Norwegian Earth System Model
PICO	Potsdam Ice-shelf Cavity mOdel
PICOP	Potsdam Ice-shelf Cavity mOdel + Plume
PISM	Parallel Ice Sheet Model
RCP	Representative Concentration Pathway
RMSE	Root Mean Square Error
SIA	Shallow-Ice Approximation
SLE	Sea Level Equivalent
SSA	Shallow-Shelf Approximation
SSP	Shared Socioeconomic Pathways
UFEMISM	Utrecht Finite Element Ice-Sheet Model
UKESM	UK Earth System Model
VAF	Volume Above Floatation
VE	Visco Elastic
WACCM	Whole Atmosphere Community Climate Model

Appendix

Appendix B - Additional Figures Including All Ensemble Members

This appendix replicates figures from the main text including all ensemble members.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Appendix

Appendix C - Description of Ice Flow Model Configuration and Main Characteristics

The descriptions below summarize the main characteristics, model parameters, and initialization procedure used by the different ice sheet modeling groups.

C1

DC_ISSM

The DC_ISSM ice sheet model configuration is similar to the JPL_ISSM model used in initMIP-Antarctica and ISMIP6 Antarctica 2100 (Seroussi et al., 2019, 2020), and uses the Ice-sheet and Sea-level System Model (Larour et al., 2012). The initial ice sheet geometry is based on BedMachine Antarctica v2 data set (Morlighem et al., 2020) to simulate an Antarctic ice sheet close to present-day conditions. The mesh resolution varies between 2 km along the grounding lines and 50 km inland, including a resolution of 8 km or less on all current ice shelves, and remains fixed over time. The stress balance is simulated using the Shelfy-Stream Approximation (MacAyeal, 1989) over the entire model domain and basal sliding follows a Budd sliding law (Budd et al., 1979), with the effective pressure equal to the ice overburden pressure, assuming the subglacial water is connected to the ocean. The model is initialized using DA to reproduce surface velocity (Rignot et al., 2011) and infer for the basal friction coefficient on grounded ice and the ice rheology on the ice shelves (MacAyeal, 1993; Morlighem et al., 2010). The rheology on the grounded ice is estimated using a thermal model. The two-dimensional model is vertically extruded into 15 layers with thinner layers close to the bed. We calculate the ice temperature assuming a thermal steady-state (Dawson et al., 2022; Seroussi et al., 2013), and using a three dimensional higher-order stress balance approximation (Blatter, 1995; Pattyn, 2003). The thermal boundary conditions are geothermal heat flux from Maule et al. (2005) at the ice base and imposed surface temperatures at the ice surface (Lenaerts et al., 2012). Steady-state temperatures are vertically averaged to calculate the ice viscosity, which is held constant over time. After this initialization, we relax the model geometry, grounding lines and velocity for 2 years (Gillet-Chaulet et al., 2012; Seroussi et al., 2011). The grounding lines evolve freely assuming hydrostatic equilibrium. The friction at the grounding line follows a sub-element scheme (SEP2 in Seroussi, Morlighem, Rignot, et al. (2014)) and the ocean-induced melt at the grounding line also follows a sub-element scheme (Seroussi & Morlighem, 2018). The ice front position is based on the levelset method (Bondzio et al., 2016); it retreats when ice becomes thinner than 10 m in any part of the domain and cannot readvance. The ice shelf melt follows the non-local quadratic scheme developed for ISMIP6 (Jourdain et al., 2020). The surface mass balance and ocean conditions in the ice shelf cavities for the control experiment are respectively from the 1979–2010 mean of RACMO2.1 (Lenaerts et al., 2012) and from the 1995–2018 ocean climatology (Jourdain et al., 2020).

C2

DOE_MALI

MPAS-Albany Land Ice (MALI; MPAS: Model for Prediction Across Scales) (Hoffman et al., 2018) is a three-dimensional, thermomechanically coupled, higher-order ice sheet model. It solves Blatter-Pattyn approximation to the Stokes equations (Blatter, 1995; Pattyn, 2003) for the momentum balance on an unstructured triangular Delaunay mesh using the finite element method. Mass and tracer transport take place on the dual Voronoi mesh using the Finite Volume (FV) method. We solve the energy balance using a temperature-based approach, which is more stable than the enthalpy-based method and yields very similar results. Mass and tracer advection use a first-order upwind scheme. We use a forward Euler time integration scheme for transport, and backward Euler for vertical heat diffusion. We use Nye's generalization of Glen's flow law (Glen, 1955; Nye, 1957) for the constitutive equation, with the effect of temperature on ice viscosity following Paterson and Budd (1982).

We used two different variable resolution meshes: an 8–30 km mesh (98,341 cells) and a 4–20 km mesh (385,379 cells), which both used five terrain-following vertical layers. Resolution is dependent both on observed ice flow speed and distance from the modern grounding line. For the speed-based cell spacing function, minimum cell spacing is used where log₁₀(speed [m/year]) ≥ 2.5, linearly decreasing to the maximum cell spacing where log₁₀(speed [m/year]) ≤ 0.5. For the distance-to-grounding-line cell spacing function, minimum cell spacing is applied within 40 km of the grounding line, linearly increasing to the maximum spacing ≥250 km from the grounding line. The final cell spacing is taken as the element-wise minimum of these two functions. This ensures high resolution over modern ice shelves and ice streams, with lower resolution toward the ice sheet interior.

To initialize our model, we simultaneously solve for the spatially-varying basal friction law coefficient and a spatially-varying depth-averaged ice stiffness parameter that minimized the misfit to observed ice velocities using the adjoint optimization method of Perego et al. (2014). The temperature field for both the 4–20 km and 8–30 km resolution meshes is taken from the initialization on the 8–30 km mesh. The basal shear stress boundary condition is a power-law sliding relationship (Weertman, 1957). We used an exponent of 1/3 for the optimization, but we recalculated the basal friction coefficient for a range of exponents from 1 to 1/10 for testing in forward simulations. We found an exponent of 1/5 to give the best trade-off between matching overall ice-sheet mass balance and preventing unrealistic rapid mass loss from Thwaites Glacier, so we use this value in all forward simulations. We used the following data sets in the inversion: BedMachine Antarctica v2 (Morlighem et al., 2020) ice thickness and bed topography (bed topography under Lake Vostok was raised to make the ice grounded); MEaSUREs 1996–2016 composite InSAR-Based Antarctica Ice Velocity Map, Version 2 (Mouginot et al., 2012, 2017; Rignot, 2017; Rignot et al., 2011); geothermal flux from Martos et al. (2017); 1979–2010 mean surface air temperature surface air temperature from RACMO 2.1 (Lenaerts et al., 2012).

We impose a fixed calving front, meaning that ice that advances beyond the initial ice extent is immediately calved away, but no other iceberg calving is imposed and the ice terminus can retreat if local mass balance thins ice thickness to zero. We use the non-local quadratic sub-shelf melt parameterization from Jourdain et al. (2020). In our 4–20 km simulations, we used the median parameter value from the MeanAnt calibration, while for the 8–30 km simulations we ran separate ensembles with both the median and 95th percentile values. We use the melt parameterization of Rignot et al. (2016) to impose melt at grounded marine margins, using the thermal forcing value from the seafloor and making the conservative assumption of zero subglacial discharge. Our historical and control simulations use mean 1995–2017 RACMO 2.3p2 surface mass balance (van Wessem et al., 2018) and the 1995–2017 observational ocean thermal forcing provided by ISMIP6 (Jourdain et al., 2020).

In preliminary simulations, we found that the grounding line advanced rapidly in many areas where it is known to be relatively stable in observations. We therefore lowered the modern seafloor by 20 m everywhere except the ASE to prevent spurious advance. While this is an ad-hoc solution, 20 m is well within the reported uncertainty in seafloor elevation at most locations near the grounding line. We also removed some ridge artifacts from several submarine troughs in the Weddell Sea Sector. Interpolation to our meshes also removed the pinning point from the Thwaites Ice Shelf, so we raised the bed topography in order to create a pinning point in the observed location with 40 m thickness above floatation. We found that the Thwaites Glacier grounding line retreated rapidly in all preliminary simulations on the 4–20 km mesh. In order to alleviate this, we created a thickness field for the ASE using the BedMap2 (Fretwell et al., 2013) surface elevation product and the BedMachine v2 bed topography. Because the BedMap2 surface corresponds to an early 2000s ice geometry, this results in much slower grounding line retreat at Thwaites Glacier during the historical period, which falls within the reported range of mass loss from Rignot et al. (2008). This change did not alleviate retreat of Thwaites Glacier in the 8–30 km simulations were complete, so this change was not included in those simulations.

Finally, we performed a 10-year simulation with melting in the ASE reduced to zero to remove fast transients from the model. The model state at the end of this simulation is considered our initial condition, with a nominal year of 2000.

C3

IGE_ElmerIce

The ice-sheet model configuration from IGE uses version v9.0 of the Elmer/Ice finite element model (Gagliardini et al., 2013). The ice dynamics is computed by solving the SSA of the Stokes equations (MacAyeal, 1989), assuming an isotropic rheology following Glen's flow law (Glen, 1955) and a linear friction law. The grounding line position is determined using a floatation criterion and the SEP3 sub-grid scheme is applied for the friction in partially floating elements (Seroussi, Morlighem, Rignot, et al., 2014). The mesh resolution varies from 1 km both close to the grounding line and in areas where observed surface velocities and thickness show high curvatures, to 50 km in the very interior of the ice sheet. The model domain does not change over time as we assume a steady calving front. The ice thickness is subject to a lower limit of 1 m and elements that reach this limit are considered deglaciated. For stability reasons, the thickness of isolated icebergs (ice-covered area with less than seven elements disconnected from the ice sheet) is set to the critical thickness of 1 m if they appear during the simulation. The model is initialized using DA. The initial present-day ice-sheet geometry is based on BedMachine Antarctica v2 (Morlighem et al., 2020), slightly modified as Lake Vostok is not considered. Initial viscosity and friction parameter are inferred using inverse methods (Brondex et al., 2019; Gillet-Chaulet et al., 2012) to match the 2015–2016 surface ice velocities described in Mouginot, Rignot, Bjørk, et al. (2019). If the ice gets grounded in areas where it was floating in the initial geometry, the friction parameter is set to 1 Pa m⁻¹ year. To match the rate of ice-sheet mass change estimated by the IMBIE team (Shepherd et al., 2018), the inverted friction parameter is reduced by 10% all over Antarctica. The ice viscosity is constant, that is, not affected by potential changes in temperature or damage. To reduce the inconsistencies between input data and inverted data when we switch from a diagnostic to a prognostic simulation (Gillet-Chaulet et al., 2012), we conducted a relaxation phase over the 1995–2014 period under the present-day forcing that we considered as our historical run. The reference present-day Surface Mass Balance is based on the 1995–2014 climatology of the RACMO-2.3.p2 regional climate model (van Wessem et al., 2018). The PICO box model (Reese et al., 2018) is used to parameterize ice-shelf basal melting under fully floating elements, with the same parameters as in Reese et al. (2023). The present-day sea floor temperature and salinity for each of the 19 regions defined in Reese et al. (2018) are extracted from the ISMIP6 ocean climatology (Jourdain et al., 2020) and averaged within 50 km of the ice-shelf front as described in Burgard et al. (2022). The resulting ocean temperature is corrected in individual regions to match the range of the 1994–2018 melt estimates from Adusumilli et al. (2020), and the correction is kept constant in all simulations.

C4

ILTS_SICOPOLIS

The model SICOPOLIS version 5.3 (SICOPOLIS Authors, 2022) is applied to the Antarctic ice sheet with hybrid shallow-ice–shelfy-stream dynamics for grounded ice (Bernales et al., 2017) and shallow-shelf dynamics for floating ice. Ice thermodynamics is treated with the melting-CTS (cold-temperate transition surface) enthalpy method by Greve and Blatter (2016). 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 (Sato & Greve, 2012) and subglacial hydrology (Calov et al., 2018; Kleiner & Humbert, 2014). The model is initialized by a paleoclimatic spin-up over 140,000 years until 1990, forced by Vostok δD converted to ΔT (Petit et al., 1999), in which the topography is nudged toward the present-day topography to enforce a good agreement (Rückamp et al., 2019). The basal sliding coefficient is determined individually for the 18 IMBIE-2016 basins (Rignot & Mouginot, 2016) 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 (Fortuin & Oerlemans, 1990), the present-day precipitation is by Arthern et al. (2006) and Le Brocq et al. (2010), and present-day and past runoff is modeled by the positive-degree-day method with the parameters by Sato and Greve (2012). The 1960–1989 average surface mass balance correction that results diagnostically from the nudging technique is used as a prescribed surface mass balance correction for the future climate simulations. The bed topography is Bedmap2 (Fretwell et al., 2013), the geothermal heat flux is by Martos et al. (2017), and isostatic adjustment is included using an elastic-lithosphere–relaxing-asthenosphere (ELRA) model (parameters by Sato and Greve (2012)). Ice-shelf basal melting is parameterized by the quadratic non-local ISMIP6 approach (Jourdain et al., 2020). The set-up is essentially the same as that used for ISMIP6 Antarctica 2100 (Seroussi et al., 2020) and the follow-up studies by Chambers et al. (2022) and Greve et al. (2023). For a more detailed description, see Greve et al. (2020).

C5

IMAU_UFEMISM

The Utrecht Finite Volume Ice-Sheet Model (UFEMISM; Berends et al., 2021) is a vertically integrated ice-sheet model. The stress balance is described by the hybrid shallow ice/SSA (Bueler & Brown, 2009). The equations are solved on an unstructured triangular grid, with a medium-high resolution near the grounding line (either 30 km or 16 km depending on the model version), and a low resolution (approximately 200 km) over the ice-sheet interior and the open ocean. Basal sliding is described by a Budd-type sliding law, and the thermomechanical coupling is realized by a 3-D thermodynamical module, affecting the ice viscosity by means of an Arrhenius-type relation. The basal friction coefficient is scaled with the sub-grid grounded fraction of each velocity cell, allowing the model to resolve the grounding line to within a single grid cell. The model is initialized with a hybrid DA/spin-up approach, based loosely on Pollard and DeConto (2012b). Here, the subglacial bed friction coefficients underneath grounded ice, and the sub-shelf ocean temperatures underneath floating ice, are progressively nudged until the modeled ice sheet reaches a steady state that matches the BedMachine Antarctica v1 geometry (Morlighem et al., 2019). Both the ice geometry and ice temperature are allowed to evolve freely throughout this procedure, which lasts for 50,000 model years. During this time, the atmospheric forcing is kept constant at the ensemble mean climate for the 1950–1979 reference period, of the global climate models that were used for the ISMIP6 Antarctica 2100 projections (Seroussi et al., 2020). In the last 10,000 years of the initialization, the modeled ice volume does not change by more than 1 cm SLE, and the root-mean-square error of the modeled ice thickness with respect to BedMachine v1 is approximately 100 m. During the initialization, no basal melt parameterization is applied; instead, basal melt rates are determined entirely by the nudging procedure. During the historical period and the projections, basal melt is calculated as a quadratic function of the local thermal forcing (ocean temperature at the ice base relative to the pressure- and salinity-dependent melting point) following Jourdain et al. (2020). The model does not include a calving law, except for ice beyond the present-day ice front (which is not allowed). Within that front, grid cells can only become ice-free through basal or surface melt.

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 Ritz (1992) and Ritz et al. (1997) and the ice shelf model of Rommelaere and MacAyeal (1997), extended to the case of ice streams treated as dragging ice shelves (Ritz et al., 2001). Over the whole domain, the velocity field consists of the superposition of the SIA velocities for ice flow due to vertical shearing and the SSA velocities, which are used as a sliding law (Bueler & Brown, 2009). Here we used the GRISLI version 2.0 (Quiquet et al., 2018), which includes the analytical formulation of Schoof (2007) to compute the flux at the grounding line. Basal drag is computed with a power law basal friction (Weertman, 1957). For this study, we use an iterative inversion method to infer a spatially variable basal drag coefficient that ensures an ice thickness that is as close as possible to observations with a minimal model drift (Le clec'h et al., 2019). 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 (van Wessem et al., 2018). Sub-shelf basal melting rates are computed with the non-local quadratic parametrization suggested in ISMIP6. 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 data set (Fretwell et al., 2013), and the geothermal heat flux is from Shapiro and Ritzwoller (2004). We use an identical setup as for ISMIP6 Antarctica (Seroussi et al., 2020) which is fully described in Quiquet and Dumas (2021). The submission LSCE_GRISLI2 corresponds exactly to the version used for (Seroussi et al., 2020) while in LSCE_GRISLI we account for oceanic basal melt in partially floating grid cell.

C7

NCAR_CISM

The Community Ice Sheet Model (CISM, Lipscomb et al., 2019, 2021) uses finite-element methods to solve a depth-integrated higher-order approximation (Goldberg, 2011) over the entire ice sheet. The model uses a structured rectangular grid with uniform 4-km horizontal resolution and five vertical σ-coordinate levels. The grounding line location is determined using hydrostatic equilibrium and a sub-element parameterization, with basal melt applied in partly floating grid cells in proportion to the floating fraction (Leguy et al., 2021). Submission 1 is configured to be consistent with the NCAR_CISM submission to Seroussi et al. (2020), whereas submission 2 uses some different settings. Basal friction combines power-law and Coulomb behavior, following either Schoof (2005) (submission 1) or Zoet and Iverson (2020) (submission 2). Sub-shelf melting is computed using the ISMIP6 quadratic nonlocal scheme (submission 1) or quadratic local scheme (submission 2) (Jourdain et al., 2020). The geothermal heat flux is from Shapiro and Ritzwoller (2004) (submission 1) or Martos et al. (2017) (submission 2). The ice sheet is initialized as in Lipscomb et al. (2021) with present-day geometry (Morlighem et al., 2020) and an idealized temperature profile, then spun up for 20,000 years using 1979–2016 climatological surface mass balance and surface air temperature from RACMO2 (van Wessem et al., 2018). During the spin-up, basal friction parameters (for grounded ice) and thermal forcing correction parameters δT (for floating ice) are adjusted to nudge the ice thickness toward present-day observations. For submission 2 we modified the initialization to assimilate recent observations of ice-shelf mass loss (Smith et al., 2020). During the spin-up, we apply a positive mass balance correction in basins (including Thwaites and Pine Island) where floating ice is observed to be thinning on average. This correction is removed at the start of each forward run, such that floating ice will thin in agreement with observations.

C8

NORCE_CISM

The CISM version run at NORCE is identical to NCAR_CISM, but some other modeling choices have been taken. The basal sliding relation is similar to that of Zoet and Iverson (2020) in all experiments. The model has been applied with a range of horizontal resolutions of 4 km (CISM2), 8 km (CISM3), 16 km (CISM4) or 32 km (CISM5) and uses three different sub-shelf melt parameterizations (quadratic nonlocal-slope, Lipscomb et al. (2021) as the default and quadratic nonlocal, quadratic local as indicated in Table 3). The ice sheet is initialized similar to NCAR_CISM with present-day geometry (BedMachine V2) and an idealized temperature profile, but then spun up for 10,000 years using 1979–2017 climatological surface mass balance and surface air temperature from MAR v3.6.4 (Agosta et al., 2019). In one case (NORCE_CISM4_JRA), we use MAR downscaling of JRA-55 (Harada et al., 2016; Kobayashi et al., 2015) and in all other cases ERA-Interim (Dee et al., 2011). Each model version of the NORCE_CISM ensemble uses a different combination of resolution, input surface mass balance and melt parameterization and receives an individual spin-up. After the spin-up, the model is relaxed for 1,000 years, the result is assigned to the year 1995, and the historical period is run under NorESM1-M anomalies until end 2014.

C9

PIK_PISM

Projections with the thermomechanically coupled Parallel Ice Sheet Model (PISM; Bueler & Brown, 2009; Winkelmann et al., 2011; ) use the model setup of Reese et al. (2023), which is based on PISM release v1.2.2. All simulations are performed on a regular rectangular horizontal grid of 8 km resolution and with 121 vertical layers, spaced quadratically from 13 m at the ice base to 100 m at the surface. In PISM, the SIA and the two-dimensional shallow-shelf approximation (Bueler & Brown, 2009; MacAyeal, 1989) of the stress balance are combined over the entire ice sheet. A generalized power law (Schoof & Hindmarsh, 2010) is applied to parameterize basal sliding. The basal friction coefficient depends on the effective pressure and till friction angle, that is parameterized using a heuristic, piecewise linear function of the bed elevation (Martin et al., 2011). The movement of the grounding line is parameterized at a subgrid scale and its position results from hydrostatic equilibrium. No subgrid interpolation of melt is used, that is, melt is applied in all cells that are floating according to the floatation criterion. The ice front (also simulated at subgrid scale) evolves freely, with the calving rate following the “eigencalving” law (Levermann et al., 2012). This is combined with a minimum thickness criterion of 50 m at the calving front, and the removal of ice extending beyond present-day extents of the Antarctic Ice Sheet and ice shelves as given in Morlighem et al. (2020). Sub-shelf melt rates are computed by the PICO (Reese et al., 2018). We use the coefficients from Reese et al. (2023), which are the result of a new parameter optimization that uses melt sensitivity estimates as a target. Glacial isostatic adjustment is not included in the projections presented here. Surface mass balance and atmospheric temperatures do not adapt to changes in ice-sheet geometry. A quasi-equilibrium initial ice-sheet state under constant historical climatic boundary conditions is obtained in a spin-up approach. Following a thermal equilibrium spin-up under constant geometry (with ice thickness and bed topography from Morlighem et al. (2020)) on a coarser 16 km grid, a full-physics spin-up ensemble with varying model parameters on 8 km horizontal resolution is run for 25,000 years. Present-day climatologies for the ocean and atmosphere are modified with NorESM1-M anomalies (Bentsen et al., 2013; with respect to the period 1995–2014) to derive climate conditions around 1850 (Reese et al., 2023). Surface mass balance and atmospheric temperatures are from RACMO2.3p2 (averaged over the period 1995–2014; van Wessem et al., 2018). For the ocean, we use observed present-day temperatures and salinities (derived from Schmidtko et al. (2014)), combined with basin-wide temperature corrections that match aggregated melt rates close to present-day (Reese et al., 2023). A historic simulation is run from 1850 to 2015 using the NorESM1-M forcing in the ocean and atmosphere. Members of the full-physics spin-up ensemble are compared to observations in terms of present-day ice thickness, ice-stream velocities, as well as deviations in grounded and floating area, and the average distance to the observed grounding-line position (Morlighem et al., 2020; Mouginot, Rignot, & Scheuchl, 2019), using the scoring methods described in Albrecht et al. (2020) and Reese et al. (2020). The ensemble member performing best in the aggregated score was chosen as the initial state for projections.

C10

UCM_Yelmo

UCM_Yelmo is a finite difference thermomechanical ice-sheet-shelf model (Robinson et al., 2020, 2021). The model configuration covers the whole Antarctic domain with 381 × 381 grid cells of 16 × 16 km resolution and 20 vertical layers in sigma-coordinates with finer spacing at the ice base. The stress balance is simulated via the higher-order depth-integrated viscosity approximation (DIVA; Goldberg & Sergienko, 2011; Lipscomb et al., 2019). The model is initialized with the present-day ice-sheet geometry (BedMachine v2, Morlighem et al., 2019) and spun up for 20,000 years with 1979–2016 climatological surface accumulation and surface air temperature from RACMO2.3 (van Wessem et al., 2018) and ocean temperature and salinity from the ISMIP6 1995–2018 ocean climatology (Jourdain et al., 2020). For the first 15,000 years of the spin-up, basal friction coefficients and the thermal forcing under the ice shelves are optimized to fit present-day ice thickness observations (similar to the procedure defined by Lipscomb et al. (2021)). The remaining 5,000 years of the spinup allow the ice sheet to evolve freely and ensure that the ice sheet has come sufficiently close to thermodynamic equilibrium. The grounding-line position is derived via the floatation criterion. Basal friction at the ice base is parameterized via a regularized-Coulomb law scaled by effective pressure and with a transition speed of 100 m/year as suggested by Zoet and Iverson (2020). Effective pressure is determined as a function of till water saturation (Bueler & van Pelt, 2015) and thus is coupled to the thermodynamics of the model. The friction at the grounding line as well as the ocean-induced melt is weighted by the floating fraction of the grid cell (PMP following the notation of Leguy et al. (2021)). Calving at the ice front is parameterized using a von-Mises-like calving law (following Lipscomb et al. (2021)) and the ice front is free to evolve. This results in slightly larger ice shelves than observations. If the temperature at the ice base reaches the pressure melting point, then the temperature is set to the pressure melting point, and the basal mass balance is diagnosed as Cuffey and Paterson (2010), where the geothermal heat flux field is obtained by Shapiro and Ritzwoller (2004). The glacial isostatic adjustment is computed with the ELRA method (Le Meur & Huybrechts, 1996), where the relaxation time of the asthenosphere is set to 3,000 years.

C11

UCSD_ISSM

The UCSD_ISSM ice sheet model configuration is similar to the configuration of the UCI_ISSM ice sheet model used in the previous ISMIP6 Antarctica 2100 ensemble (Seroussi et al., 2020). The model domain covers the present-day Antarctic Ice Sheet as defined in BedMachine Antarctica v2 (Morlighem et al., 2020) and is discretized into an unstructured triangular mesh with edge lengths varying from 1 km in dynamic ice streams up to 50 km in regions of stagnant flow. The governing ice flow equations are solved through implementation of the two-dimensional Shelfy-Stream stress balance approximation to the Stokes flow equations (MacAyeal, 1989). We initialize our ice sheet model with static inversions of interferometric synthetic aperture radar (InSAR) derived surface ice velocities from 2011 to obtain the basal friction coefficient for grounded elements and ice stiffness for floating elements (Morlighem et al., 2013; Rignot et al., 2011). These fields remain unchanged throughout all projections. Basal friction is solved with a Budd sliding law (Budd et al., 1979). Ice viscosity is computed assuming an ice temperature of −10°C using the table provided in Cuffey and Paterson (2010). We apply water pressure at the ice-ocean interface and a stress-free boundary condition at the ice-air interface. Bed topography and initial ice geometry is taken from BedMachine Antarctica v2 (Morlighem et al., 2020) and the present day surface mass balance field is from the Regional Atmospheric Climate Model-2.3 (van Wessem et al., 2018). We use a floatation criterion to determine the position of the grounding line in our model simulations. The precise location of the grounding line within individual elements is tracked using the “SEP2” sub-element grounding line migration parameterization, in which basal friction is integrated only over the grounded portion of a mesh element (Seroussi, Morlighem, Larour, et al., 2014). Ice shelf basal melting is computed with the PICOP ice shelf melt rate parameterization (Pelle et al., 2019), which is a combination of an ocean box model and a buoyant plume parameterization that uses basin-averaged ocean condition on the seafloor adjacent to the ice sheet edge as inputs. Ice shelf melt is applied only to elements that are completely floating, as to follow the recommendations of Seroussi and Morlighem (2018). Lastly, the ice front remains fixed in all projections, as we do not model iceberg calving.

C12

ULB_Kori

The Kori-ULB ice flow model (Coulon et al., 2023; Klose et al., 2023), previously called f.ETISh (Pattyn, 2017), is a vertically integrated hybrid finite-difference ice sheet/ice shelf model with vertically integrated thermomechanical coupling. Grounded ice flow is represented as a combination of the SSA for basal sliding and SIA for ice deformation, while SSA is applied for floating ice shelves (Bueler & Brown, 2009; Winkelmann et al., 2011). Transient thermodynamics are solved in 3d, using shape functions for the vertical profile of horizontal and vertical velocities (Lliboutry, 1979). A flux condition (related to the ice thickness at the grounding line; Schoof, 2007) is imposed at the grounding line following the implementation by Pollard and DeConto (2012a, 2020), as described in Pattyn (2017). All simulations are performed on a regular rectangular horizontal grid of 16 km resolution. Basal sliding is introduced as a Weertman sliding law, with sliding exponent m = 3. Basal melting underneath the floating ice shelves is determined by different sub-shelf melt schemes: the main model submission (ULB_Kori1) uses the PICO model (Reese et al., 2018) with $${\gamma }_{T}^{\star }=3\times 1{0}^{-5}\,\mathrm{m}\,{\mathrm{s}}^{-\mathrm{1}}$$ and C = 2 × 10⁶ m⁶ s⁻¹ kg⁻¹ while the second model submission (ULB_Kori2) uses the ISMIP6 quadratic non-local parameterization with median values (Jourdain et al., 2020). Calving at the ice front depends on the combined penetration depths of surface and basal crevasses, relative to total ice thickness, following Pollard et al. (2015). Note that the influence of hydrofracturing is not considered here. Prescribed input data include the present-day ice-sheet geometry and bedrock topography from the BedMachine v2 data set (Morlighem et al., 2020) and the geothermal heat flux by Shapiro and Ritzwoller (2004). Present-day mean surface mass balance and temperature are obtained from van Wessem et al. (2018), based on the regional atmospheric climate model RACMO2.3p2. Changes in bedrock elevation due to changes in ice load are modeled by the commonly used ELRA model, using an asthenosphere relaxation time τ of 3,000 years and a flexural rigidity of the lithosphere D of 10²⁵ N m (Le Meur & Huybrechts, 1996). Ice-sheet initial conditions and spatially-varying basal sliding coefficients are obtained by an inverse method following the nudging scheme by Pollard and DeConto (2012b), using surface mass balance forcing for the year 1950 (anomalies from 1945 to 1955 relative to the period 1995–2014 derived from the CMIP5 climate model NorESM1-M are added to the present-day climatology for 1995–2014). In the inverse procedure, basal sliding coefficients under grounded ice, and sub-shelf melt rates under floating ice (Bernales et al., 2017) are adjusted iteratively to reduce the misfit with observed ice thickness (the calving front is kept to its observed present-day position). The resulting sub-shelf melt rates, treated as balance melt rates, remain independent of ocean boundary conditions. To limit an initial shock caused by the transition from the balance sub-shelf melt rates to the imposed sub-shelf melt parameterization scheme, a short 10-year relaxation is run after the model initialization, before the historical simulation, using constant atmospheric and oceanic forcings for the year 1950 (Coulon et al., 2023). Hindcasts from 1950 to 2014 are produced using changes in oceanic and atmospheric boundary conditions derived from NorESM1-M (Bentsen et al., 2013).

C13

UNN_Úa

The Úa ice flow model is a vertically integrated finite element model (Gudmundsson, 2020). The SSA flow approximation was used, although the model allows for the use of other flow approximations as well, such as SIA and hybrid formulations. The momentum and mass conservation equations are solved simultaneously using an implicit approach. Ice thickness was set to a minimum value of 0.1 m and this thickness constraint enforced using the active set method. Vertically integrated density is allowed to vary horizontally, and the resulting impact on the force term is included. Various basal sliding laws are implemented, and here Weertman sliding law was used. For given basal stress (m) and ice flow exponents (n), inverse methods are used to solve simultaneously for basal drag (C), and the ice rheological factor (A) in Glen's flow law. Inversion was performed once at the starting point of the runs, and no spin-up period was used, and no adjustments made to surface mass balance. Surface velocity data was used for the inversion and the geometry based on BedMachine v2 (Morlighem et al., 2020). Resolution was spatially variable and ranged from 1 to 40 km depending on location. Linear triangular elements were used. Around grounding lines resolution was generally about 1 km and the coarsest resolution was used in the interior of East Antarctica. Calving front positions were held fixed.

C14

UTAS_ElmerIce

The Elmer/Ice (Gagliardini et al., 2013) model domain covers the present-day Antarctic ice sheet, and its geometry is interpolated from the MEaSUREs BedMachine Antarctica, Version 3 data set (Morlighem et al., 2020). An unstructured finite-element mesh is created following coastline positions provided by BedMachine V3, and refined using MMG (Dobrzynski & Frey, 2008) based on geometry gradients and the Hessian matrix of the observed velocity fields provided by MEaSUREs Phase-Based Antarctica Ice Velocity Map, Version 1 (Mouginot, Rignot, Bjørk, et al., 2019). Mesh resolution in the horizontal varies from approximately 1 km near the grounding lines of fast-flowing ice streams to approximately 25 km in the interior. For the stress balance calculation, we use the SSA (MacAyeal, 1989) to solve the vertically integrated ice dynamics. An inverse method is used to simultaneously infer both basal drag coefficient and enhancement factor through reducing the mismatch between the simulated and observed MEaSUREs velocity. A linear Weertman sliding law (Weertman, 1957) is applied in the inversion and further used in forward simulations, during which the ice front is not allowed to evolve. For partially floating cells, no melting is applied and the basal drag is calculated on a sub-grid scale using 20 integration points inside a partly grounded element. We impose a minimum ice thickness of 40 m everywhere in the domain. The surface mass balance used in the surface relaxation and the control experiment is forced by the 1995–2014 mean from the MAR model (Agosta et al., 2019). Basal melt rates are computed using the local quadratic parameterization provided by ISMIP6 as an alternative to the nonlocal parameterization. A historical run from 1995 to 2014 is conducted with forcing from NorESM1-M RCP8.5.

C15

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 (Huybrechts, 1990, 2002) and paleoclimatic studies (Van Breedam et al., 2021, 2023). It considers thermomechanically coupled flow in both the ice sheet and the ice shelf, using the SIA and SSA coupled across a one grid cell wide transition zone where all the stress component contribute in the effective stress in the flow law. 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 Giovinetto and Zwally (2000) compilation used in Huybrechts et al. (2000), updated with accumulation rates obtained from shallow ice cores during the EPICA pre-site surveys (Huybrechts, 2007). 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 (Janssens & Huybrechts, 2000). 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 (Beckmann & Goosse, 2003). 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 (Goelzer et al., 2016). 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 ISMIP6 protocol.

C16

VUW_PISM

Using PISM v2.0.3, we followed the approach described in Golledge et al. (2019) and Seroussi et al. (2020). Starting from initial bedrock and ice thickness conditions from Morlighem et al. (2020), together with reference climatology from van Wessem et al. (2014), we ran a multistage spinup that guarantees a 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. We started with an initial 32 km resolution 20-year smoothing run in which only the SIA is used. Then, holding the ice geometry fixed, we ran a 250,000 years, 32 km resolution, thermal evolution simulation in which temperatures were allowed to equilibrate. Refining the grid to 16 km and resetting bed elevations and ice thicknesses we ran a further 1,000 years using full model physics and a present-day climate, refined the grid to 10 km for a further 500 years and then refined the grid again to 8 km for a historical run forced by global climate model from 1950 to 2000. The resultant configuration provided the starting point for each of our simulations to 2300. In contrast to the fixed-bed approach of Seroussi et al. (2020), we used a visco-elastic Earth deformation model (Lingle & Clark, 1985) to track vertical loading and unloading as ice thickness changed. Using the supplied climate anomaly fields, we ran an ensemble of simulations that explored both parameter uncertainty as well as forcing uncertainty. The two main runs, PISM1 and PISM2, were run without and with subgrid basal melt at the grounding line, respectively. These runs used a pseudo-plastic q of 0.4 and mantle viscosity of 1e20 Pa s. Sensitivity experiments were also conducted to gauge model response to ice geometry response to the linearity of the basal traction parameter q (PISM1-s1: q = 0.3; PISM1-s2: q = 0.5) and to the mantle viscosity value used in the Earth deformation scheme (PISM1-s3: mv = 5e18; PISM1-s4: mv = 5e21). These runs were then repeated with subgrid melting applied at the grounding line (PISM2-s1 to s4). Compared to our previous simulations (Seroussi et al., 2020), our new runs employ a surface PDD standard deviation value of 10, chosen to minimize the mismatch between modeled and observed ice geometries. We also used a higher SSA enhancement factor (0.8 instead of 0.6), and explored a range of values for the pseudo-plastic q parameter (from 0.3 to 0.5). We set our eigencalving K value to 1e16, and a minimum thickness calving threshold of 50.

Acknowledgments

We thank the Climate and Cryosphere (CliC) effort, which provided support for ISMIP6 through sponsoring of workshops, hosting the ISMIP6 website and wiki, and promoted ISMIP6. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modeling, coordinated and promoted CMIP5 and CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the CMIP data and providing access, the University at Buffalo for ISMIP6 data distribution and upload, and the multiple funding agencies who support CMIP5 and CMIP6 and ESGF. We thank the ISMIP6 steering committee, the ISMIP6 model selection group and ISMIP6 data set preparation group for their continuous engagement in defining ISMIP6. This is ISMIP6 contribution No 33. Hélène Seroussi was supported by grants from NASA Cryospheric Science Program (#80NSSC21K1939 and #80NSSC22K0383) and the Novo Nordisk Foundation under the Challenge Programme 2023—Grant number NNF23OC00807040. Jake Twarog was supported by a URAD grant from Dartmouth College's Undergraduate Advising & Research program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. Tyler Pelle was supported by grants from the NASA Cryospheric Science Program (#80NSSC22K0387), the NSF (#OPP-2114454), and the Cecil H. and the Ida M. Green Foundation for Earth Sciences at the Institute of Geophysics and Planetary Physics at Scripps Institution of Oceanography. WHL was supported by the NSF National Center for Atmospheric Research, which is a major facility sponsored by the U.S. National Science Foundation under Cooperative Agreement No. 1852977. Support for Xylar Asay-Davis, Holly Kyeore Han, Trevor Hillebrand, and Matthew Hoffman was provided through the Scientific Discovery through Advanced Computing (SciDAC) and Early Career Research programs, funded by the US Department of Energy (DOE), Office of Science, Advanced Scientific Computing Research and Biological and Environmental Research programs. MALI simulations were performed on machines at the National Energy Research Scientific Computing Center, a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award using NERSC awards ERCAP0024081 and ERCAP0023782. Heiko Goelzer has received funding from the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement No. 869304, PROTECT and from the Research Council of Norway under project 343397 (CLIM2Ant). Rupert Gladstone was supported by Academy of Finland Grants 322430 and 355572, and by the Finnish Ministry of Education and Culture and CSC—IT Center for Science (Decision diary number OKM/10/524/2022). Chen Zhao, Yu Wang, and Ben Galton-Fenzi received grant funding from the Australian Government as part of the Antarctic Science Collaboration Initiative program (ASCI000002). Chen Zhao is the recipient of an Australian Research Council Discovery Early Career Researcher Award (project number DE240100267) funded by the Australian Government. Thomas Zwinger was supported by the Academy of Finland Grant 286587. UTAS_ElmerIce simulations were enabled by computational resources provided by CSC-IT Center for Science Ltd. Ralf Greve was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant JP16H02224, JP17H06104, and JP17H06323. Alexander Robinson received funding from the European Union (ERC, FORCLIMA, 101044247). Marisa Montoya received funding from the Spanish State Research Agency (MARINE, PID2020-117768RB-I00). Jorge Alvarez-Solas received funding from the Spanish State Research Agency (ICEAGE, PID2019-110714RA-I00). Javier Blasco received funding from the Dutch Research Council (HiRISE, OCENW.GROOT.2019.091). Torsten Albrecht was funded by the German climate modeling project PalMod (FKZ: 01LP1925D and 01LP2305B) supported by the German Federal Ministry of Education and Research (BMBF) as a Research for Sustainability initiative (FONA). Torsten Albrecht acknowledges support by OCEAN:ICE, which is co-funded by the European Union, Horizon Europe Funding Programme for research and innovation under grant agreement Nr. 101060452 and by UK Research and Innovation. This is O:I contribution number 8. Ann Kristin Klose, Ronja Reese, Julius Garbe and Torsten Albrecht gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for supporting this project by providing resources on the high performance computer system at the Potsdam Institute for Climate Impact Research (PIK). Development of PISM is supported by NASA Grants 20-CRYO2020-0052 and 80NSSC22K0274 and NSF Grant OAC-2118285. Ann Kristin Klose acknowledges support by the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 820575 (TiPACCs) and No. 869304 (PROTECT). Julius Garbe, Ronja Reese, David Chandler, Petra Langebroek, Benoît Urruty, Fabien Gillet-Chaulet, Pierre Mathiot, Sainan Sun, G. Hilmar Gudmundsson and Nicolas Jourdain were supported by the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement No. 820575 (TiPACCs). Petra Langebroek was also supported by the Research Council of Norway through the Centre of Excellence iC3 (project number 332635). Justine Caillet, Fabien Gillet-Chaulet, Pierre Mathiot and Nicolas Jourdain were supported by the French National Research Agency (ANR) under grant ANR-19-CE01-0015 (EIS). IGE team work was granted access to the high-performance computing (HPC) resources of TGCC under allocation A0140106066 attributed by GENCI. Constantijn J. Berends and Violaine Coulon were supported by PROTECT, receiving funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 869304. This is PROTECT contribution number 113. Jorge Bernales was funded by NWO, Grant OCENW.KLEIN.515. UFEMISM simulations were carried out on the Dutch National supercomputer Snellius, hosted at SURF. Nick Golledge and Dan Lowry are supported by Antarctic Science Platform (ANTA1801) and Our Changing Coast (RTVU2206) grants from the New Zealand Ministry of Business, Innovation, and Employment. Computational resources to perform Kori simulations have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant 2.5020.11 and by the Walloon Region. Philippe Huybrechts and Jonas Van Breedam acknowledge support from the Research Foundation Flanders (FWO Vlaanderen) with Grant G091820N. Roderik van de Wal received funding from the Dutch Polar Program (Grant: ALWPP.2019.003).

Data Availability Statement

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 at ). Scalars computed from two dimensional fields for this study, routines used to compute them, and scripts for analysis and figures are permanently available on Zenodo: (Seroussi & Pelle, 2024). Original forcings data sets and simulations results for two-dimensional fields are available on Ghub () under “ISMIP6 23rd Century Forcing Data sets”: and (Nowicki & ISMIP6 Team, 2024a). Model outputs from the simulations described in this paper are available on Ghub () under “ISMIP6 23rd Century Projections”: and (Nowicki & ISMIP6 Team, 2024b).

References

Adhikari, S., Ivins, E., Larour, E., Seroussi, H., Morlighem, M., & Nowicki, S. (2014). Future Antarctic bed topography and its implications for ice sheet dynamic. Solid Earth, 5(1), 569–584. [DOI: https://dx.doi.org/10.5194/se-5-569-2014]