Experimental protocol for sea level projections from ISMIP6 stand-alone ice sheet models

Projection of the contribution of ice sheets to sea level change as part of the Coupled Model Intercomparison Project Phase 6 (CMIP6) takes the form of simulations from coupled ice sheet–climate models and stand-alone ice sheet models, overseen by the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6). This paper describes the experimental setup for process-based sea level change projections to be performed with stand-alone Greenland and Antarctic ice sheet models in the context of ISMIP6. The ISMIP6 protocol relies on a suite of polar atmospheric and oceanic CMIP-based forcing for ice sheet models, in order to explore the uncertainty in projected sea level change due to future emissions scenarios, CMIP models, ice sheet models, and parameterizations for ice–ocean interactions. We describe here the approach taken for defining the suite of ISMIP6 stand-alone ice sheet simulations, document the experimental framework and implementation, and present an overview of the ISMIP6 forcing to be used by participating ice sheet modeling groups.

1 Introduction

The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) is a targeted activity of the Climate and Cryosphere (CliC) project of the World Climate Research Project (WCRP) and has been formally endorsed by the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016). Its aim is to provide process-based projections of the sea level contributions from the Greenland and Antarctic ice sheets that are tightly linked to the wider suite of CMIP6 climate projections and employ forcing from the CMIP atmosphere–ocean general circulation model (AOGCM) ensemble. Nowicki et al. (2016) describe the overall design of ISMIP6, which includes an assessment of the impact of initial conditions on projections (ISMIP6-initMIP; Goelzer et al., 2018; Seroussi et al., 2019), experiments in which ice sheet models are fully coupled within Earth system models (ISM–ESMs), and experiments with stand-alone ice sheet models (ISM) forced by output from the CMIP AOGCM.

ISMs numerically simulate the dynamic flow of glacial ice – from the continental interior, across geographical distances, and potentially transitioning to floating ice shelves, prior to terminating into the global ocean. An ISM receives input fields in the form of accumulating snowfall, surface and ocean temperature, and other time-varying conditions describing the ice sheet surface and its lateral boundaries, and it provides output fields of ice velocities and the distribution of ice mass. In a coupled ISM–ESM framework, the input fields become functions of how the ice sheets vary with time. Changes in the ice sheet topography can influence the atmospheric circulation, while the selective discharge of ice and meltwater may alter ocean circulation. The coupling of dynamical ice sheet models with ESMs is highly complex, as a mismatch between the relatively high spatial resolution and long integration time step of the ISM and the relatively coarse spatial resolution and short integration time step of the ESM atmosphere and ocean fields must be negotiated (Vizcaino, 2014; Fyke et al., 2018). To date, only a limited number of ESMs have been coupled with ISMs, which motivates the need for simulations with state-of-the-art stand-alone ISMs in order to explore the uncertainty in projected sea level change.

Nowicki et al. (2016) described a number of possible avenues for the stand-alone component of ISMIP6 but were limited in their final protocol design because, at the time, the CMIP6 simulations had not started. Subsequent to the publication of Nowicki et al. (2016), ISMIP6 formed focus groups to evaluate the polar climate in the CMIP AOGCMs and to finalize how the output of CMIP AOGCMs would be translated into forcing for ice sheet models. Here we present a revised version of the ISMIP6 protocol, based on the improved understanding of ISM needs and CMIP6 AOGCM outputs. The complex issues of providing the offline AOGCM output at a high spatial resolution suitable for ISM modeling needed to be addressed in a uniform, standardized manner that would allow broad participation from the current generation of ice sheet models. Specific challenges included (i) the translation of the various AOGCM resolutions and grids to the various ISM grid resolutions; (ii) the poor representation of steep gradients in the surface topography of the ice sheet margins within AOGCMs, which underestimates large gradients in atmospheric forcing; (iii) the quality of AOGCM polar climate; and (iv) the mismatch in the spatial extents of ice sheets within the AOGCMs and initialized ice sheet extent within certain ISMs. Additionally, oceanic variables from AOGCMs needed to be extrapolated from continental shelves to provide boundary conditions underneath ice shelves and at the calving front, as AOGCMs typically do not resolve ice shelf cavities and proglacial fjords. The ability to provide boundary conditions to ISM modeling groups in a timely manner is another factor that influenced the final ISMIP6 protocol design. The particular implementation of ice–atmosphere and ice–ocean interactions within each participating ISM had to be considered, and these ISM-specific constraints on the protocol were guided by lessons learned from the ISMIP6 initMIP efforts (Goelzer et al., 2018; Seroussi et al., 2019). Finally, the protocol needed to explore the uncertainty in sea level projections due to the choice of emissions scenario, the choice of CMIP AOGCM, the ice sheet model physics (structural uncertainty), and how poorly constrained parameters within the model are (parameter uncertainty), and the uncertainty due to formulation of ice–ocean interactions.

This paper describes the detailed experimental protocol used for stand-alone experiments with ice sheet models of Greenland and Antarctica using forcing from CMIP AOGCMs, and it presents the novel atmospheric and oceanic forcing datasets prepared by ISMIP6. We begin by providing an overview of the projection framework and the purpose of the simulations in Sect. 2. We next present the protocol for initializing the projections, including schematic experiments needed to explore the uncertainty in ice sheet evolution due to initial state and historical simulations in Sect. 3. The atmospheric forcing and implementations are then described in Sect. 4, the oceanic forcing and implementation in Sect. 5, and Antarctic ice shelf fracture strategy in Sect. 6. We summarize the protocol and discuss the expected outcomes and impacts of ISMIP6 in Sect. 7.

2 Overview of the projection setup and their purpose

Following the CMIP6 protocol (Eyring et al., 2016), all ISMIP6 projections start in January 2015 and end in December 2100. Although extensions beyond 2100 are available for some climate models in the CMIP5 archive and possible in the CMIP6 ScenarioMIP protocol (O'Neill et al., 2016), ISMIP6 focus remains on the end of the 21st century, as it is constrained by the availability and quality of the polar climate forcing. Beyond 2100, the ice sheets' surface elevation will likely have deviated significantly from the fixed ice sheet elevation configuration used by the CMIP models, which may affect projected polar climate.

The projection setup strategy for the ISMIP6 stand-alone ISMs is illustrated in Fig. 1. For a given CMIP AOGCM and future climate scenario, an ISM atmospheric forcing dataset ultimately takes the form of surface mass balance (SMB) and surface temperature, while oceanic forcing data include oceanic temperature and salinity in order to infer Antarctic ice shelf basal melt and Greenland calving and frontal melt. For the Antarctic ice shelves, atmospheric properties may result in surface melting, which in turn can trigger ice shelf collapse (Trusel et al., 2015). An issue faced by ISMIP6 is the mismatch in spatial resolution and spatial extent between available AOGCM fields and ISM needs: AOGCMs do not generally resolve oceanic flow within the Greenland fjords or beneath the Antarctic ice shelves, and SMB varies rapidly over the steep topography at the ice sheet margins, but these SMB gradients are not captured by AOGCMs. The implication is that extrapolation and downscaling of atmospheric and oceanic AOGCM fields may be required to produce realistic ice sheet projections. For the Greenland ice sheet atmospheric fields, downscaling was done via the use of a regional climate model. For both Greenland and Antarctica, far-field ocean temperature and salinity were extrapolated through Greenland fjords and beneath Antarctic ice shelves using rules that account for the blocking effects of bathymetric sills.

Figure 1 Overview of ISMIP6 (a) Antarctic and (b) Greenland projection framework, illustrating the strategy for translating the CMIP atmospheric and oceanic properties into climatic forcing for ice sheet models.

Table 1 CMIP model used to obtain atmospheric and oceanic forcing for the Greenland and Antarctic ice sheets.

Ice–ocean interactions remain an active area of research, and the current generation of ice sheet models use a variety of representations. For example, the Antarctic ice–ocean interaction representations range from simple linear relationships between oceanic temperature change and melt rate (e.g., Rignot and Jacobs, 2002) to more complex parameterizations (see Favier et al., 2019, for a review). Similarly, the Greenland ice–ocean interactions in large-scale ice sheet models range from ad hoc methods (e.g., Price et al., 2011; Bindschadler et al., 2013; Goelzer et al., 2013; Nick et al., 2013; Fürst et al., 2015) to estimated submarine melt rates (Aschwanden et al., 2019). The ISMIP6 protocol allows for sampling the uncertainty in the representation of ice–ocean interaction via the use of “open” and “standard” experiments (Sect. 5). In the standard experiments, ice sheet models implement parameterizations and forcings designed by the ocean focus groups. The standard experiments further test the parameter uncertainty in the basal melt formulations by sampling “low”, “mid”, and “high” values of the parameters in the melt rate parameterizations for both ice sheets. In addition, for the Antarctic ice sheets, two calibrations are investigated (see Jourdain et al., 2020): one based on observed mean sub-shelf basal melt over Antarctica (MeanAnt) and one based sub-shelf basal melt near the grounding line of Pine Island glacier (Pine Island grounding line calibration, PIGL). In the open experimental setup, ice sheet models can use their existing parameterization for ice–ocean interaction, driven by the provided extrapolated AOGCM oceanic forcing datasets. The open experiments further allow for inclusion of additional physical processes that are not taken into account in the standard framework. An example is exploring the Greenland ice sheet response to surface meltwater reaching the base of the ice sheet and affecting basal sliding. As there is disagreement on the implication of this process for ice sheet evolution (Shannon et al., 2013; Koziol and Arnold, 2018), this forcing is not part of the standard experiments, but could be explored in the open approach.

Table 2 Antarctic core (Tier 1) experiments.

Table 3 Greenland core (Tier 1) experiments.

Future climate scenarios are defined in CMIP5 as Representative Concentration Pathways (RCPs; van Vuuren et al., 2011) and in CMIP6 as Shared Socioeconomic Pathways (SSPs; Eyring et al., 2016; O'Neill et al., 2016). The ISMIP6 protocol samples the potential sea level from ice sheets under two different climate scenarios: the high-emission CMIP5 RCP8.5 and CMIP6 SSP5-8.5 scenarios are our primary focus, but the lower-emission scenarios CMIP5 RCP2.6 and CMIP6 SSP1-26 are included for a few CMIP models. The focus on CMIP5 scenarios in this revised ISMIP6 protocol is due to the delay in CMIP6 model simulations, which prevents a full analysis of the CMIP6 models. The use of the CMIP multi-model mean for forcing ISMs is not feasible in our experimental protocol, as it is not possible to downscale atmospheric fields via regional climate models (RCMs) with a multi-model mean. It is also not feasible to use all of the CMIP models, due to both the time and computational effort needed to prepare the forcing dataset, as well as the time and computational effort required for running the ISMs. The ISMIP6 strategy is to sample the CMIP ensemble in order to select a manageable number of CMIP models for our projections that is representative of the spread in the full CMIP ensemble. As described in Barthel et al. (2020), six CMIP5 AOGCMs were selected per ice sheet based on the following criteria: (i) present-day polar climate in agreement with observations (evaluated by model biases over the historical period; for example Agosta et al., 2015), (ii) sampling a diversity of future climates (evaluated by difference in projections and code similarities), and (iii) a focus on models with RCP8.5 and RCP2.6 which also have the fields required for RCM downscaling. The CMIP5 models were selected independently for the Greenland and Antarctic ice sheets, using atmospheric and oceanic metrics appropriate for each ice sheet. Four of the CMIP5 models were chosen to be the same for both ice sheets (CSIRO-MK3.6.0, HadGEM2-ES, IPSL-CM5A-MR, NorESM1-M), the fifth choices were closely related (MIROC5 for Greenland and MIROC-ESM-CHEM for Antarctica), and the sixth choices were unrelated (ACCESS1.3 for Greenland and CCSM4 for Antarctica). For more information, see Table 1. The key characteristics of models taking part in CMIP5 are summarized in Table 9.1 and 9.A.1 of Flato et al. (2013). As CMIP6 models started to become available in late spring–summer 2019, ISMIP6 selected four CMIP6 AOGCMs (CESM2, CNRM-CM6-1, CNRM-ESM2-1, UKESM1-0-LL) based solely on their availability and the fact that two of these models would be taking part in the coupled climate–ice sheet component of ISMIP6, which allows for future scientific analysis of the difference in projection from stand-alone ISMs versus ISMs that are fully coupled to climate models.

Table 4 Initialization experiments and examples of duration of experiment for different choices for initial state and projection start state. Modeling groups are free to choose initial state dates that are not indicated in this table, but the projection start state should always correspond to January 2015.

^∗ Denotes experiments that are only needed if the initial state used for the projections is different than that submitted for the initMIP effort. The ctrl may need to be extended from an original initMIP submission.

The projection protocol consists of “core”, or Tier 1, experiments, which modeling groups are required to perform, and “targeted experiments”. The targeted experiments are optional and further divided into higher-priority Tier 2 and lower-priority Tier 3 experiments. Core experiments are designed to explore the range of CMIP5 model uncertainty with three AOGCMs under two future emissions scenarios, the impact of Antarctic ice shelf fracture, and the uncertainty in ocean parameters for groups that participate in the standard oceanic protocol. Groups that can participate in both the open and standard implementations of ocean forcing are encouraged to do so. This setup results in five open and eight standard projections for Antarctica (Table 2) and four open and six standard experiments for Greenland (Table 3). The Tier 2 experiments consist of the remainder of the three CMIP5 AOGCMs and the three CMIP6 models. Tier 3 experiments further explore the uncertainty in the standard ocean parameterization, the impact of ice shelf fracturing, and simulations driven by atmosphere-only forcing (no change in ocean) and ocean-only forcing (no change in atmosphere) to help understand the source of mass loss from the corresponding full simulations, in which changes are due to both atmosphere and ocean. The complete list of Tier 2 and Tier 3 experiments is presented in Appendix A. This mix of core and targeted experiments follows the approach taken in Shannon et al. (2013): it provides a flexible framework that allows less computationally expensive models to explore the full set of experiments, allowing modelers to choose to focus on certain aspects of the protocol that fit their research interests, while ensuring that all groups perform a subset of identical experiments.

Participating groups may decide to investigate the impact of ice sheet model uncertainty in projections via different model setup choices. These include exploring mesh resolution as well as model parameterizations such as the basal sliding law, parameters in ice sheet flow approximation, and ice shelf basal melt parameterization. Unlike the original initMIP effort, SMB and bedrock adjustment in response to a changing ice sheet is allowed. In some cases, it may be necessary to treat these modeling setup choices as different models and repeat the initialization method, the initMIP experiments, and the historical and control runs described below.

Ice sheet modeling groups are requested to submit the variables listed in Appendix B, as long as these variables are applicable to their models. These consist of state variables (such as ice thickness), flux variables (such as SMB), and integrated scalar values (such as total ice-sheet-wide SMB flux). To facilitate the analysis and intercomparison, groups should save their model output on one of the ISMIP6 grids with spatial resolution closest to the model native grid (see Appendix B for details).

3 Initial state, control runs, and historical run

Ice sheet model initial states are typically obtained via different methods: long interglacial spinup or data assimilation of present-day observations, as well as hybrid combinations of these two methods (Nowicki et al., 2013a, b; Goelzer et al., 2017; Pattyn et al., 2017). Interglacial spinups have the advantage of obtaining an ice sheet that can capture transients due to past climatic conditions but the disadvantage of producing an ice sheet geometry that may differ from the present day (Seroussi et al., 2013). Assimilation methods, on the other hand, capture the present-day geometry, but projections often suffer from unrealistic drifts due to the model responding to inconsistencies in the input datasets. Time-dependent data assimilation methods allow for more realistic transients, but to date they have been limited to regional studies or synthetic ice sheet setup (e.g., Goldberg et al., 2015; Gillet-Chaulet, 2020). Other methods include combinations of these techniques, as demonstrated in the initMIP Greenland and Antarctica efforts (Goelzer et al., 2018; Seroussi et al., 2019). These multiple approaches for initialization, and the use of different observations for assimilation methods, create a challenge in the design protocol, as the initial state date becomes model specific and ranged from the 1990s to the 2010s in the initMIP efforts, for example.

Modeling groups taking part in the ISMIP6 projections are free to reuse their “initial state” submitted as part of initMIP or create a new initial state. In the latter case, groups are asked to rerun the 100-year-long initMIP schematic experiments – anomalies in surface mass balance (“asmb”) for Greenland and Antarctica and anomalies in ice shelf basal melt (“bsmb”) for Antarctica – as well as a control run (“ctrl”). This control run is needed to capture any drift present in the projections as a result of the initialization method. The control run is implemented as a forward run without any anomaly forcing, such that for example any surface mass balance used in the initialization would continue unchanged. The control run begins from the same initial state as the initMIP schematic experiments and the historical run and lasts as long as the initMIP experiments and the projections. The control run may need to be extended from an original initMIP submission (where the control was set to 100 years) in order to cover the time period from the initial state to the end of the projections. Table 4 illustrates typical setup and time span for the ISMIP6 protocol.

The “historical run” bridges the gap between the time of the ice sheet initial state and the “projection start state” in January 2015. The projections then branch from a single historical run for each ice sheet model. Because the time of each ice sheet model's initial state varies, the duration of the historical run will therefore also vary between models. Ice sheet modeling groups are left to decide how to perform the historical run and bring their models to the projection start state, a choice motivated by (i) the distinct initialization procedures used in the ice sheet modeling community; (ii) the lack of a known set of historical atmospheric and oceanic forcing that can reproduce observed changes, due in part to the limited observational record; and (iii) the challenges associated with our revised strategy of using multiple CMIP models and scenarios to sample the uncertainty in future climate. The latter would require multiple historical runs from each ice sheet model, which may then result in a distinct projection start state for a given ice sheet model, complicating the projection forcing strategy as well as interpretation of the simulations. Nonetheless, AOGCM-derived historical datasets are provided for each AOGCM projection dataset. Modeling groups are free to use one of the AOGCM historical dataset, or a reanalysis, or a combination of multiple datasets. To test the impact of the choice of historical dataset on the projections and associated model drift, groups are required to submit a “projection control”. This simulation is also an unforced ice sheet model run (implemented with zero anomalies), which starts from the ice sheet projection start state (January 2015) and runs until December 2100.

4 Atmospheric forcing and implementation

Atmospheric forcing for stand-alone ice sheet model simulations consists of surface mass balance (SMB) and surface temperature derived from CMIP AOGCMs. SMB provides mass gain from accumulation (snow and rain) and mass loss from ablation or surface melting to the ice sheet model. Current AOGCM outputs can be directly used to compute SMB, but this approach does not capture the narrow peripheral region with steep SMB gradients, which are key for ice sheet simulations. In this case, AOGCM climate can be downscaled using high-resolution regional climate models (RCMs), which is the technique used for ISMIP6 Greenland projections. RCMs currently provide more realistic surface climate than the direct output of CMIP5 AOGCMs for both ice sheets (Fettweis et al., 2013; Noël et al., 2018; van Wessem et al., 2018; Agosta et al., 2019). In the future, it may be possible to bypass the use of RCMs for downscaling, as a new generation of climate models have implemented multiple elevation classes (CESM, Lipscomb et al., 2013; UKESM1, Sellar et al., 2019; ModelE; Fischer et al., 2014), allowing SMB to be computed at multiple elevations within a horizontal grid cell in order to capture SMB gradients, as well as improvements in the parameterization of polar surface processes in AOGCMs (e.g., Cullather et al., 2014; van Kampenhout et al., 2017; Alexander et al., 2019). The use of CMIP5 models for this protocol therefore requires AOGCMs to be selected for their skills in simulating forcing fields for RCMs for both Greenland and Antarctica, instead of their skills in simulating SMB. Due to time constraints and computational demands with the use of RCMs, the atmospheric forcing for the ice sheet simulations was only obtained with an RCM for the Greenland ice sheet and derived directly from the AOGCM output for the Antarctic ice sheet simulations.

Regardless of the methodology chosen to obtain surface forcing for the ice sheet models from CMIP AOGCM, the strategy of investigating the uncertainty due to CMIP climate models and scenarios (which requires the use of multiple CMIP AOGCMs), as well as the distinct initialization methods used by ice sheet models (which uses diverse SMB and temperature sources), prevents the direct application of SMB and surface temperature as a boundary condition for ice sheet models. Instead, surface forcings are implemented via annual SMB and temperature anomalies. This choice assumes that the interannual variability in SMB and temperature anomalies is greater than any differences in SMB and temperature climatologies from both the CMIP AOGCMs and the climatologies used in the initialization of ice sheet models.

Figure 2 Surface temperature anomalies (K) over the Antarctic ice sheet under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios. (a) Time series of mean surface temperature anomaly from 1950 to 2100 for all AOGCMs and scenarios selected and (b–e) surface temperature anomaly over the time period 2081–2100 for (b) NorESM1-M under RCP2.6, (c) NorESM1-M under RCP8.5, (d) MIROC-ESM-CHEM under RCP8.5, and (e) CCSM4 under RCP8.5.

Figure 3 Surface mass balance anomalies (mm w.e. yr⁻¹) over the Antarctic ice sheet under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios. (a) Time series of mean surface mass balance anomaly from 1950 to 2100 for all AOGCMs and scenarios selected and (b–e) mean surface mass balance anomaly over the time period 2081–2100 for (b) NorESM1-M under RCP2.6, (c) NorESM1-M under RCP8.5, (d) MIROC-ESM-CHEM under RCP8.5, and (e) CCSM4 under RCP8.5.

4.1 Antarctic atmospheric forcing and implementation

For the Antarctic ice sheet, ISMIP6 provides yearly averaged anomalies in SMB (computed as precipitation minus evaporation minus runoff) and surface temperature, along with the respective climatologies for SMB and temperature used to compute the anomalies. The data were prepared by ISMIP6 using direct input from the CMIP models listed in Table 1 and interpolated onto the ISMIP6 input grids. More information on how the datasets were prepared is available in Appendix C. The anomalies were computed as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{aSMB}}_{\mathrm{AOGCM}}\left(x,y,t\right)& ={\mathrm{SMB}}_{\mathrm{AOGCM}}\left(x,y,t\right)\\ & -{\mathrm{SMB}}_{\mathrm{CLIM},\mathrm{AOGCM}}\left(x,y\right)\end{array}\end{array}$$

and

$$\begin{array}{}\text{(2)}& \mathrm{a}{T}_{\mathrm{AOGCM}}\left(x,y,t\right)=T_{\mathrm{AOGCM}}\left(x,y,t\right)-T_{\mathrm{CLIM},\mathrm{AOGCM}}\left(x,y\right),\end{array}$$

where ${\mathrm{aSMB}}_{\mathrm{AOGCM}}\left(x,y,t\right)$ and a$T_{\mathrm{AOGCM}}\left(x,y,t\right)$ are the anomalies in SMB and temperature, ${\mathrm{SMB}}_{\mathrm{AOGCM}}\left(x,y,t\right)$ and $T_{\mathrm{AOGCM}}\left(x,y,t\right)$ are the annual SMB and temperature for a given AOGCM, while ${\mathrm{SMB}}_{\mathrm{CLIM},\mathrm{AOGCM}}\left(x,y\right)$ and $T_{\mathrm{CLIM},\mathrm{AOGCM}}\left(x,y\right)$ are the climatologies. The climatologies were computed by taking the mean values of SMB_AOGCM and T_AOGCM over the Antarctic reference period (January 1995 to December 2014). The anomaly datasets cover the time period of 1950 to 2100.

During the simulations, ice sheet models need to reintroduce the climatology that best fits their simulation and compute surface input (${\mathrm{SMB}}_{\mathrm{ISM},\mathrm{AOGCM}}\left(x,y,t\right)$ and $T_{\mathrm{ISM},\mathrm{AOGCM}}\left(x,y,t\right))$ as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathrm{SMB}}_{\mathrm{ISM},\mathrm{AOGCM}}\left(x,y,t\right)& ={\mathrm{SMB}}_{\mathrm{REF}}\left(x,y\right)\\ & +{\mathrm{aSMB}}_{\mathrm{AOGCM}}\left(x,y,t\right)\end{array}\end{array}$$

and

$$\begin{array}{}\text{(4)}& T_{\mathrm{ISM},\mathrm{AOGCM}}\left(x,y,t\right)=T_{\mathrm{REF}}\left(x,y\right)+\mathrm{a}{T}_{\mathrm{AOGCM}}\left(x,y,t\right),\end{array}$$

where SMB_REF(x,y) and T_REF(x,y) are the SMB and temperature that the ice sheet model would have used over the reference period, which is the same for all core and targeted experiments. SMB_REF(x,y) and T_REF(x,y) should be computed as the average over the reference period for time-dependent SMB and temperatures, or simply set to the climatology used in the case of time-independent input. Note that the anomalies are constant over the entire year and change step-wise at the beginning of every year.

Figure 4 Surface temperature anomalies (K) over the Greenland ice sheet under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios. (a) Time series of mean surface temperature anomaly from 1950 to 2100 for all AOGCMs and scenarios selected and (b–e) surface temperature anomaly over the time period 2081–2100 for (b) MIROC5 under RCP2.6, (c) MIROC5 under RCP8.5, (d) NorESM1-M under RCP8.5, and (e) HadGEM2-ES under RCP8.5.

Figure 5 Surface mass balance anomalies (mm w.e. yr⁻¹) over the Greenland ice sheet under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios. (a) Time series of mean surface mass balance anomaly from 1950 to 2100 for all AOGCMs and scenarios selected and (b–e) mean surface mass balance anomaly over the time period 2081–2100 for (b) MIROC5 under RCP2.6, (c) MIROC5 under RCP8.5, (d) NorESM1-M under RCP8.5, and (e) HadGEM2-ES under RCP8.5.

The surface temperatures and SMB fields to be used in forcing the ISM simulations are illustrated in Figs. 2 and 3 respectively. Time series of values averaged over the fixed Antarctic ice sheet mask under RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios for all Tier experiments are shown, while the spatial patterns of mean anomalies from 2081 to 2100 are shown for Tier 1 experiments only (Table 2). These datasets were generated by ISMIP6 and it is their first presentation in the literature. We note that CCSM4 has a finer native model resolution ($\mathrm{0.94}{}^{\circ }\times \mathrm{1.25}{}^{\circ }$ grid; see Table C1 in Appendix C) compared to NorESM1-M ($\mathrm{1.91}{}^{\circ }\times \mathrm{2.50}{}^{\circ }$) and MIROC-ESM-CHEM ($\mathrm{2.81}{}^{\circ }\times \mathrm{2.81}{}^{\circ }$). As expected, the SMB and temperature anomalies are correlated: CMIP models that project warmer surface conditions (Fig. 2a) also project an increased SMB (Fig. 3a). From the Tier 1 simulations, the largest temperature increases are found in CCSM4 and MIROC-ESM-CHEM fields under the RCP8.5 scenario. The temperature anomalies from these models are about 3.6 to 6.4 K warmer in the 2081–2100 period than at the start of the projection. This contrasts with the temperature fields of the selected NorESM1-M CMIP5 model, which indicates negligible temperature increases by the end of the century for the RCP2.6 scenario and a 2 to 3.8 K increase for the RCP8.5 scenario. As with the continent-averaged time series, the spatial patterns of temperature anomalies for the NorESM1-M RCP2.6 output for the 2081 to 2100 period are similarly muted (Fig. 2b), with the greatest warming over the peninsula, Ronne and Amery ice shelves, and the east Antarctic plateau under RCP8.5 (Fig. 2c). The MIROC-ESM-CHEM and CCSM4 models (Fig. 2d, e) under the RCP8.5 scenario project similar spatial patterns of temperature anomalies as with the NorESM1-M, but with greater magnitude. The MIROC-ESM-CHEM output indicates greater warming than for the CCSM4 over the peninsula, the ice shelves adjacent to Dronning Maud Land, and Enderby and Kemp Land glaciers in East Antarctica and over the Getz ice shelf and the ice shelves fed by Marie Byrd Land glaciers in West Antarctica. CCSM4 shows less warming over the steep margins of the ice sheet.

The coarser grid of MIROC-ESM-CHEM impacts the spatial distribution of projected SMB anomalies (Fig. 3d), with negative values over Marie Byrd Land, the Getz and Abbot ice shelves, glaciers feeding the Ronne and Filchner ice shelves in East Antarctica, the peninsula, the Brunt ice shelf, Enderby Land, and on both flanks of the Amery ice shelf, due to large surface runoff in these regions. The projected SMB anomalies in MIROC-ESM-CHEM are negligible over the East Antarctic plateau and increase towards the coast. This pattern of minimal SMB change over the interior plateau and increased SMB for coastal regions is also captured in CCSM4 and the NorESM1-M models under RCP8.5 (Fig. 3c, e), with the particular exception of the northern Antarctic Peninsula, where the anomalies are negative due to sufficiently warm conditions that promote ablation. The relatively fine grid used by CCSM4 allows for more spatial variations in SMB anomalies, particularly over the steep ice sheet margins. Aside from negative SMB anomalies over the tip of the peninsula, the NorESM1-M model indicates little change with time in conditions under RCP2.6 (Fig. 3b), which is also reflected in the mean ice-sheet-wide SMB anomalies (Fig. 3a). The large areas with negative SMB anomalies projected by the MIROC-ESM-CHEM simulation balance with increases in other regions such that the ice-sheet-averaged SMB anomalies are comparable to those of the NorESM1-M model under RCP8.5. CCSM4 is the Tier 1 model with largest ice-sheet-wide projected increase in SMB, which is comparable to the Tier 2 CMIP6 models under SSP5-8.5 (Fig. 3a).

Figure 6 (a) Bathymetry modified from Fretwell et al. (2013) and IMBIE2 basins (Shepherd et al., 2018) used in the sub-ice-shelf extrapolation of oceanic conditions. (b) Extrapolated ocean temperature from the observational climatology.

4.2 Greenland atmospheric forcing and implementation

For the Greenland ice sheet, ISMIP6 generated surface forcing from CMIP AOGCMs that have been reinterpreted through the Modèle Atmosphérique Régionale (MAR) regional climate model (version 3.9.6; Delhasse et al., 2020). Although MAR uses a fixed topography for the Greenland ice sheet, it also allows for SMB–height feedback to be included in the ice sheet model simulations by providing vertical SMB and temperature gradients on each horizontal grid cell (Franco et al., 2012). The vertical gradients are also used to downscale the original MAR results computed on a 15 km grid to a finer 1 km grid, allowing for the resolution of steep topography. In addition, MAR calculates potential SMB and temperature in regions that are outside the MAR ice sheet mask, allowing for surface forcing to be computed for ice sheet spatial extent that differs slightly from the MAR ice sheet mask.

For each MAR downscaled CMIP model, ISMIP6 provides annual anomalies in SMB and surface temperature from 1950 to 2100, the vertical SMB and temperature gradients over the same time period, and the respective climatologies computed over the reference period of January 1960 to December 1989. The reference period is distinct from the one for Antarctica, because the Greenland ice sheet is considered to have been in steady state with the climate during this period (e.g., Fettweis et al., 2017; Mouginot et al., 2019). The surface forcing anomalies were computed in a manner similar to Antarctica (Eqs. 1, 2) and need to be added to the ice sheet model climatology during the simulation. The implementation of surface forcing during the ice sheet simulation differs slightly from the setup for Antarctica, as the MAR-computed vertical SMB and temperature gradient fields allow for implementation of SMB and temperature feedback in the Greenland framework. During the ice sheet projection, surface forcings are implemented as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{SMB}}_{\mathrm{ISM},\mathrm{RCM}}\left(x,y,t\right)& ={\mathrm{SMB}}_{\mathrm{REF}}\left(x,y\right)\\ & +{\mathrm{aSMB}}_{\mathrm{RCM}}\left(x,y,t\right)\\ & +{\displaystyle \frac{{\mathrm{dSMB}}_{\mathrm{RCM}}\left(x,y,t\right)}{\mathrm{d}z}}\\ & \left(h_{\mathrm{ISM}}\left(x,y,t\right)-h_{\mathrm{REF}}\left(x,y\right)\right)\end{array}\end{array}$$

and

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{T}_{\mathrm{ISM},\mathrm{RCM}}\left(x,y,t\right)& =T_{\mathrm{REF}}\left(x,y\right)+\mathrm{a}{T}_{\mathrm{RCM}}\left(x,y,t\right)\\ & +{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{RCM}}\left(x,y,t\right)}{\mathrm{d}z}}\left(h_{\mathrm{ISM}}\left(x,y,t\right)\right.\\ & \left.-h_{\mathrm{REF}}\left(x,y\right)\right),\end{array}\end{array}$$

where $h_{\mathrm{ISM}}\left(x,y,t\right)$ is the ice sheet model time-dependent surface elevation, h_REF(x,y) is the ice sheet model surface elevation at the start of the projection, $\frac{{\mathrm{dSMB}}_{\mathrm{RCM}}\left(x,y,t\right)}{\mathrm{d}z}$ the time-dependent SMB vertical gradient, and $\frac{\mathrm{d}{T}_{\mathrm{RCM}}\left(x,y,t\right)}{\mathrm{d}z}$ the time-dependent temperature vertical gradient.

The MAR SMB forcing can only be directly applied when the ice sheet model projection start state is close to the present-day geometry used in the MAR simulations. However, the initMIP-Greenland experiments (Goelzer et al., 2018) show that for some ice sheet models the present-day ice sheet can differ substantially from the observed ice sheet configuration, especially for models that initialize their ice sheet to the present day via interglacial spinup. This can also be the case for models that use assimilation techniques and a long relaxation scheme that results in large geometric changes. In these cases, the surface forcing anomalies and vertical gradients should be corrected and remapped to the modeled ice sheet, using the technique described in Goelzer et al. (2020). The method uses the strong dependence of SMB and temperature on elevation to remap the MAR field and reduce unphysical biases while preserving the overall surface forcing patterns. Once the surface forcings have been remapped by ISMIP6 to an individual ice sheet configuration, modelers should implement surface forcings in the same manner as described previously.

The MAR-derived surface temperature and SMB fields to be used in forcing the ISMs are shown in Figs. 4 and 5, respectively. Time series of ice sheet area-averaged values under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios for all Tier experiments over the 21st century were computed using the fixed Greenland ice sheet present-day area mask. The spatial patterns of anomalies averaged for the 2081–2100 period are illustrated using the Tier 1 models (Table 3). These datasets are specific to the ISMIP6 project and shown here for the first time in the literature. Although the time series of both surface temperature and SMB anomalies show considerable interannual variability, the experiments project an overall warming and an associated negative SMB trend over the margins due to increased surface runoff. At the top of the range of Tier 1 models under the RCP8.5 scenario, the HadGEM2-ES projects an increase in surface temperature of about 7 K by 2100 (Fig. 4a), with the greatest warming over the northeastern Greenland ice sheet and smaller values concentrated over the southern ice sheet margins that are characterized by steep topography (Fig. 4e). The patterns of mean SMB change for HadGEM2-ES (Fig. 5e) project increases over the interior of the ice sheet due to enhanced precipitation and large negative SMB anomaly values over the periphery on the ice sheet. This general pattern of SMB change is also found with the MIROC5 and NorESM1-M models under RCP8.5 (Fig. 5c and e), but unlike the HadGEM2-ES, positive SMB anomalies are also located over the south plateau. The region is projected to warm less in the MIROC5 and NorESM1-M projections compared to the HadGEM2-ES (Fig. 4c–e). For the entire Greenland ice sheet, NorESM1-M projects the lowest area-averaged increase in surface temperature – and the smallest decrease in SMB – compared to the MIROC5 and HadGEM2-ES integrations under RCP8.5 (Figs. 4a and 5a). This is also reflected in the spatial patterns of surface temperature (Fig. 4d), in which the NorESM1-M output indicates smaller temperature increases over the periphery of the ice sheet but temperature increases comparable to MIROC5 over the central western and southern regions. For the selected simulations of the RCP2.6 emission scenario, surface temperatures are generally projected as slightly increased compared to present-day conditions (Fig. 4b) but with regional variations: for the MIROC5 output, the northern Greenland ice sheet is projected to warm more than the southern region. The output indicates an increase in SMB for the interior of the ice sheet and decrease around the periphery. Areas of negative SMB anomalies are found to extend further inland compared to the RCP8.5 integrations, despite the smaller magnitude (Fig. 5b, c).

Figure 7 Mean cavity basal melt rates (m w.e. yr⁻¹) under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios over the 21st century for MeanAnt and PIGL calibrations and nonlocal parameterizations in the Amundsen Sea sector. (a) Time series for Pine Island and Thwaites ice shelves obtained with MeanAnt (thin lines) and PIGL (thick lines) calibrations (adapted from Jourdain et al., 2020). Spatial patterns of mean sub-shelf basal melt rate from 2081 to 2100 for the Tier 1 models assuming MeanAnt calibration for (b) NorESM1-M under RCP2.6, (c) MIROC-ESM-CHEM under RCP8.5, (d) CCSM4 under RCP8.5, and (e) NorESM1-M under RCP8.5 and PIGL calibration for (f) NorESM1-M under RCP8.5. The projections from the MeanAnt and PIGL calibrations assume that the ice shelf cavities have not evolved and are set to the present day.

Figure 8 Mean cavity basal melt rates (m w.e. yr⁻¹) under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios over the 21st century for MeanAnt and PIGL calibrations and nonlocal parameterizations in the Weddell Sea sector. (a) Time series for the Ronne–Filchner ice shelf obtained with MeanAnt (thin lines) and PIGL (thick lines) calibrations (adapted from Jourdain et al., 2020). Spatial patterns of mean sub-shelf basal melt rate from 2081 to 2100 for the Tier 1 models assuming MeanAnt calibration for (b) NorESM1-M under RCP2.6, (c) MIROC-ESM-CHEM under RCP8.5, (d) CCSM4 under RCP8.5, and (e) NorESM1-M under RCP8.5 and PIGL calibration for (f) NorESM1-M under RCP8.5. The projections from the MeanAnt and PIGL calibrations assume that the ice shelf cavities have not evolved and are set to the present day.

5 Oceanic forcing and implementation

Oceanic forcing for stand-alone ice sheet model simulations consists of temperatures or melt rates at the ice–ocean interfaces (at grounded ice fronts and beneath floating ice shelves), or retreat rate for Greenland models that are not explicitly resolving calving. AOGCM output requires extrapolation under the Antarctic ice shelves and into Greenland fjords, as these regions are not resolved in CMIP models. ISMIP6 protocol allows for ice sheet models to implement their own methods for simulating ice–ocean interactions (open approach) using the ISMIP6 extrapolated datasets for oceanic conditions. ISMIP6 also proposes a standard approach, where the representation of ice–ocean interactions is specified and datasets that allow for exploring the uncertainty in the standard method are provided. This dual approach of open and standard approaches is designed to explore the impact of the range of uncertainty in ice–ocean interactions.

5.1 Antarctic oceanic forcing and implementation

For the Antarctic ice sheet, ISMIP6 provides three-dimensional anomalies of ocean ambient temperature, salinity, and thermal forcing (temperature minus freezing temperature), averaged yearly from 1850 to 2100, for the CMIP simulations listed in Table 1 and Appendix A, which were added to an observational climatology. The observational climatology was produced from a combination of the Marine Mammals Exploring the Oceans from Pole to Pole (MEOP; Roquet et al., 2013, 2014; Treasure et al., 2017), a prerelease of the World Ocean Atlas 2018 (WOA18; Locarnini et al., 2019; Zweng et al., 2019) and the Met Office ENA4 datasets (Good et al., 2013). The CMIP model climatologies, computed over the reference period from January 1995 to December 2014, are provided for legacy purposes and should not be used in ice sheet parameterizations. Instead, modelers are advised to use anomalies added to the observational climatology. All datasets were extrapolated under ice shelves using an algorithm (described in Jourdain et al., 2020) that accounts for sills and troughs, as these bathymetric features affect the flow of oceanic currents. The bathymetry (based on Bedmap2; Fretwell et al., 2013) used for the data preparation, and an example of the resulting sub-ice-shelf extrapolation of ocean temperature, is shown in Fig. 6. The extrapolated datasets are available on a 60 m vertical and 8 km horizontal ISMIP6 Antarctic grid for use by any ice sheet model implementing its own method for prescribing oceanic forcing (open approach). The extrapolation allows for oceanic fields to vary spatially beneath the ice shelves and results, for example, in ambient temperatures that vary with depth.

For the standard approach, the Antarctic ISMIP6 ocean focus group (Jourdain et al., 2020) developed datasets for two sub-ice-shelf melt rate parameterizations: nonlocal quadratic and local quadratic functions of thermal forcing. These parameterizations were evaluated in Favier et al. (2019) for an idealized Pine Island glacier geometry against other commonly used parameterizations, such as the plume parameterization of Lazeroms et al. (2018), the box parameterization of Reese et al. (2018), and a three-dimensional ocean–ice sheet coupled model. The nonlocal parameterization was found to be in closer agreement with the coupled simulations than the local parameterization and is therefore the preferred approach for ISMIP6 standard simulations. However, the nonlocal parameterization may be more complex to implement for ice sheet modeling groups as some quantities have to be averaged over regions, which might not be straightforward in large-scale parallel models. Ice sheet modeling groups are therefore free to choose between either the nonlocal or local quadratic parameterizations when participating in the standard experiments.

The sub-ice-shelf parameterizations take a slightly different form from that proposed by Favier et al. (2019), to allow for regional temperature corrections. The regional sectors are based on the IMBIE2 basins (Shepherd et al., 2018; Mouginot et al., 2017) and extrapolated to the shelf break (Fig. 6 and Jourdain et al., 2020) to allow for modeled ice sheet extent to be distinct from that of the IMBIE2 observations. The nonlocal quadratic sub-ice-shelf parameterization, m(xy), takes the form of

$$\begin{array}{}\text{(7)}& \begin{array}{rl}m\left(x,y\right)& ={\mathit{\gamma }}_{\mathrm{0}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{sw}}{c}_{\mathrm{pw}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}}}\right)}^{\mathrm{2}}\left(\mathrm{TF}(x,y,z_{\mathrm{draft}})\right.\\ & \left.+\mathit{\delta }{T}_{\mathrm{sector}}\right)\left|\langle \mathrm{TF}{\rangle }_{\mathrm{draft}\in \mathrm{sector}}+\mathit{\delta }{T}_{\mathrm{sector}}\right|,\end{array}\end{array}$$

where γ₀ is a calibration coefficient, ρ_sw and ρ_i are the seawater and ice densities, c_pw is the specific heat of seawater, L_f is the fusion latent heat of ice, TF(xyz_draft) is the thermal forcing at the ice–ocean interface, 〈TF〉_{draft∈sector} is the thermal forcing averaged over all the ice shelves in a sector, and δT_sector is a sector temperature correction. The latter is needed in order to reproduce observation-based melt rate at the sector scale and accounts for biases in observational products. The local quadratic sub-ice-shelf melt parameterization takes the form of

$$\begin{array}{}\text{(8)}& \begin{array}{rl}m\left(x,y\right)& ={\mathit{\gamma }}_{\mathrm{0}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{sw}}{c}_{\mathrm{pw}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}}}\right)}^{\mathrm{2}}\left(\max \left[\mathrm{TF}\left(x,y,z_{\mathrm{draft}}\right)\right.\right.\\ & {\left.\left.+\mathit{\delta }{T}_{\mathrm{sector}},\mathrm{0}\right]\right)}^{\mathrm{2}}.\end{array}\end{array}$$

As described in Jourdain et al. (2020), both quadratic parameterizations are calibrated against observational estimates using two methods: the “MeanAnt” and the “PIGL” methods. The MeanAnt method calibrates γ₀ and δT_sector so that the sub-ice-shelf melt parameterizations reproduce the mean Antarctic melt rates of Rignot et al. (2013) and Depoorter et al. (2013). The PIGL method calibrates γ₀ and δT_sector so that the sub-ice-shelf melt parameterizations reproduce the spatial patterns of melt rates observed close to the grounding line of Pine Island ice shelf (Rignot et al., 2013). To explore the sensitivity of ice sheet simulations to uncertainties in the basal melt rate arising from the calibration, observational melt rates were randomly sampled 10⁵ times to obtain a distribution of possible low (5th percentile), median (50th percentile), and high (95th percentile) values for γ₀ for both calibration method (MeanAnt and PIGL) and parameterization type (local and nonlocal). Once γ₀ had been determined, the median value for δT_sector was obtained from random sampling (10⁵ times) of temperature.

To highlight the differences between sub-shelf melt rates, we compare projections under RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 obtained from the nonlocal sub-shelf parameterizations with the MeanAnt and PIGL calibrations using the median values for γ₀ and δT_sector. These projections assume that the ice shelf cavities are fixed in time and set to the present day. As described in Jourdain et al. (2020), the average melt rates over 1995–2014 correspond to the observational climatology. The time series of mean cavity sub-shelf melt rate for Pine Island, Thwaites, and Ronne–Filchner ice shelves (Figs. 7a and 8a) project a gradual increase for MeanAnt calibration, with values generally less than twice the present-day conditions by 2100. In contrast, the PIGL calibration projects greater melt rates than its MeanAnt counterpart, as well as larger interannual variability. Both calibrations display regional differences due to the choice of CMIP model. For example, MIROC-ESM-CHEM projects negligible change in melt rate over the Ronne–Filchner ice shelves under RCP8.5 but is one of the highest CMIP models for Pine Island and Thwaites ice shelves with the PIGL calibration. The timing of increasing melt rate differs between the CMIP models: NorESM1-M and CCSM4 project similar melt rates with the PIGL calibration in the last 3 decades of the 21st century for both regions under RCP8.5, but the increase in melt rates from CCSM4 only begins to be significant in the 2060s for Ronne–Filchner and in the 2040s for Pine Island and Thwaites. HadGEM2-ES projects the largest increase in melt rate for both calibrations overall under RCP8.5, despite projecting (like most models) negligible change over Ronne–Filchner until the middle of the century. The only model with large Ronne–Filchner melt rate increase prior to the 2050s is UKESM1-0-LL under SSP5-8.5 and PIGL calibration. For the Pine Island and Thwaites ice shelves, the three Tier 1 models MIROC-ESM-CHEM, NorESM1-M, and CCSM4 project melt rates at the end of the 21st century under RCP8.5 that are comparable to those from the SSP5-8.5 models CNRM-ESM2-1 and UKESM1-0-LL and larger than CNRM-CM6-1 under SSP5-8.5. Over the Ronne–Filchner ice shelf, however, the Tier 1 models NorESM1-M and CCSM4 are closer to CNRM-CM6-1 under SSP5-8.5 and project a third of the melt rate from CNRM-ESM2-1 and UKESM1-0-LL.

The spatial pattern of mean sub-shelf melt rate over 2081–2100 is illustrated with the Tier 1 simulations listed in Table 2. In the Amundsen Sea sector, the largest melt rates are located over the Pine Island, Thwaites, Crosson, and Dotson ice shelves for all scenarios, CMIP models, and calibrations (Fig. 7b–f). Despite similar patterns, the amplitude of the PIGL calibration (Fig. 7f) is an order of magnitude larger than with the MeanAnt calibration (Fig. 7e). The spatial variation in projected melt rate is more apparent over the larger Ronne–Filchner ice shelf for all scenarios, with the melt rate increasing towards the ice shelf–ice sheet junctions, due to the deeper ice–ocean interface (Fig. 8b–f). As illustrated with NorESM1-M, the magnitude of melt rate from the PIGL calibration is again much larger than with the MeanAnt calibration (Fig 8e, f). As expected from the time series, the melt rates from CCSM4 (Fig. 8d) and NorESM1-M (Fig. 8e) are similar under RCP8.5 and larger than MIROC-ESM-CHEM (Fig. 8c). The latter are also smaller than those from NorESM1-M under RCP2.6 (Fig. 8a).

5.2 Greenland oceanic forcing and implementation

The ISMIP6 Greenland ocean focus group (Slater et al., 2019, 2020) proposed two methods for implementing oceanic forcing: the “retreat implementation” and the “submarine melt implementation.” The two approaches are needed in order to maximize participation from current state-of-the-art ice sheet models. The retreat implementation prescribes the temporal evolution of ice extent using masks that specify annual ice sheet extent. Both implementations rely on yearly average datasets of subglacial discharge per glacier and two-dimensional ocean thermal forcing (temperature minus freezing temperature), varying horizontally but not with depth, from 1950 to 2100 for the CMIP simulations listed in Table 3 and Appendix A. In the retreat implementation, these datasets are used to create projections of marine-terminating glacier retreat from 2015 to 2100. In the submarine melt implementation, these datasets allow ice sheet models to calculate submarine melt rate and combine with a calving rate to obtain total frontal ablation. Here we provide an overview of the ocean forcing strategy for Greenland ice sheet models taking part in ISMIP6. For details of the retreat parameterization, readers are referred to Slater et al. (2019), and full details of the two implementations can be found in Slater et al. (2020).

Subglacial discharge is freshwater that emerges from beneath the ice into proglacial fjords at marine-terminating glacier termini. The discharge is a result of surface meltwater runoff that reaches the base of the ice sheet and is routed through the subglacial hydrological system prior to reaching the ice margin. Because many ice sheet models do not include a physical representation of the surface processes that generate runoff, nor do they simulate evolving subglacial hydrology, the ISMIP6 Greenland ocean focus group (Slater et al., 2020) recommends approximating subglacial discharge as a spatial aggregation of surface runoff, calculated by a regional climate model, over a static delineation of each glacier's subglacial hydrological catchment. In other words, it is assumed that surface meltwater runoff generated within a particular catchment is instantaneously transported to the ice sheet bed and routed to the terminus, according to water routing that remains constant over the projection time period. As described in Slater et al. (2020), the annual subglacial discharge dataset was produced from downscaled MAR surface runoff (Sect. 4.2), with subglacial hydrological catchments delineated using ice sheet geometry corresponding to present-day surface elevation (Howat et al., 2014) and basal topography (BedMachine3; Morlighem et al., 2017). In addition, because the CMIP AOGCMs used as forcing in MAR may differ from present-day climate, the MAR surface runoff was bias corrected so that it better matches present-day surface runoff for each glacier from the time period between 1995 to 2014. This was done using surface runoff estimates from the RACMO2.3p2 regional climate model, which is forced by ERA-Interim atmospheric reanalysis (Noël et al., 2018), as a measure of the present-day runoff values (Slater et al., 2020). Ice sheet models with a configuration that differs substantially from these observational datasets should consider computing their own subglacial discharge dataset using the appropriate surface runoff and water-routing scheme. This guidance applies to groups that use the SMB remapping technique described in Sect. 4.2.

As described in Slater et al. (2020), the CMIP oceanic datasets have been extrapolated into fjords for use by groups choosing the submarine melt implementation or in a model's own scheme. The extrapolation takes into account the ocean bathymetry and subglacial topography following the method of Morlighem et al. (2019). The technique identifies an “effective depth”, which is the deepest point within a fjord connected with the open ocean. The method then assumes that ocean mass shallower than the effective depth is in contact with the open ocean, while water deeper than the effective depth is sheltered from the open ocean. For depths shallower than the effective depth the temperature and salinity are set to the closest ocean conditions at that particular depth, whereas for deeper regions the temperature and salinity are set to the values corresponding to the effective depth in that region. In addition, both salinity and temperatures have been bias corrected with present-day observations from the Hadley Centre EN4.2.1 dataset (Good et al., 2013). Note that in the preparation of the retreat implementation datasets (Slater et al., 2019), bathymetry was not taken into account; instead oceanic properties were averaged over the seven sectors between 200 and 500 m.

Figure 9 Retreat for selected Greenland ice sheet glaciers under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios. (a) Time series of retreat over the 21st century for the northwest sector (adapted from Slater et al., 2019) and retreat from 2015 to 2100 for (b) MIROC5 under RCP2.6, (c) MIROC5 under RCP8.5, (d) NorESM1-M under RCP8.5, and (e) HadGEM2-ES under RCP8.5. The bedrock contour from Morlighem et al. (2017) is shown as black lines; gray indicates ice sheet mask. Red numbers indicate King Oscar Glacier (51) and the unnamed glacier (53) discussed in the text.

Figure 10 Melt rate anomalies (m d⁻¹) for selected Greenland ice sheet glaciers under the RCP2.6, RCP8.5, SSP1-2.6, and SSP5-8.5 scenarios. (a) Time series of mean melt rate anomaly over the 21st century and mean melt rate anomaly over the time period 2081–2100 for (b) MIROC5 under RCP2.6, (c) MIROC5 under RCP8.5, (d) NorESM1-M under RCP8.5, and (e) HadGEM2-ES under RCP8.5. The bedrock contour from Morlighem et al. (2017) is shown as black lines; gray indicates ice sheet mask. Red numbers indicate King Oscar Glacier (51) and the unnamed glacier (53) discussed in the text. Note that for HadGEM2-ES missing pixels are due to bias correction making ocean thermal forcing in some of the cooler models slightly negative, resulting in no melt rate anomaly. The calculation of melt rate anomaly assumes that the ice sheet thickness and terminus position do not evolve with time.

In the standard approach, the retreat implementation is designed to be simple enough to be implemented by most ice sheet models (Slater et al., 2019, 2020). The approach follows the parameterization of Cowton et al. (2018) and combines anomalies in subglacial runoff and ocean thermal forcing to produce terminus retreat rates that are implemented via the use of a time-variable ice mask. As ice sheet models may not capture small outlet glaciers and may have distinct locations for individual outlet glaciers compared to the observations, it is not practical to provide a dataset for each observed outlet glacier. Instead, separate retreats are provided for seven ice–ocean sectors:

$$\begin{array}{}\text{(9)}& \mathrm{d}L=\mathit{\kappa }\mathrm{d}\left(Q^{\mathrm{0.4}}\mathrm{TF}\right),\end{array}$$

where dL is the retreat distance of each glacier in a particular sector, κ is a calibration constant described in Slater et al. (2019), and Q is the summer subglacial runoff generated from the mean of June, July, and August. The uncertainty associated with the retreat implementation is investigated via low-, median-, and high-retreat scenarios for each CMIP model. These were obtained from an ensemble of 10⁴ ice-flux-weighted trajectories for each ice–ocean sector and set to the 25th percentile (low), 50th percentile (median), and 75th percentile (high) for each CMIP model. Implementation for a specific ice sheet model requires the model's initial ice mask and the observed basal topography, in order to identify ice prone to outlet glacier retreat. For each model participating in ISMIP6 and experiment, a specific time-variable retreat mask was produced using the method described in Appendix D. These annual “land_ice_area_fraction” masks cover the period from 2015 to 2100. The mask values are set to 0.0 over ice-free regions, to 1.0 over fully ice-covered regions, and values in between for grid cells that are partially covered. Ice sheet models should apply full retreat for mask values set to 0.0, no retreat for masks set to 1.0, and partial retreat for mask values in between 0.0 and 1.0. The retreat datasets are not optimum for glaciers that have a floating ice shelf, as the datasets were calibrated using the ice front position of glaciers that do not have ice shelves. Nonetheless, for the few Greenland glaciers with floating ice shelves, it is suggested that the retreat be imposed at the ice front, rather than the grounding line. Groups that are able to compute sub-ice-shelf melt may use the basin thermal forcing datasets provided for each basin.

In the submarine melt implementation, the melt rate and calving rate along each marine-terminating glacier's terminus must be calculated by the ice sheet models and combined to determine the frontal ablation rate. Previous work has shown that submarine melt rate can be parameterized as a function of ocean thermal forcing and subglacial discharge (Xu et al., 2013; Rignot et al., 2016) and we utilize this same parameterization:

$$\begin{array}{}\text{(10)}& \dot{m}=\left(\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}h{q}^{\mathrm{0.39}}+\mathrm{0.15}\right){\mathrm{TF}}^{\mathrm{1.18}},\end{array}$$

where $\dot{m}$ is the submarine melt rate, h is the depth at the grounding line, q is the annual subglacial runoff, and TF is the ocean thermal forcing. The melt rate is the average rate of ice melt across the entire terminus. Thus, at a given ice sheet model time step, q and TF can be sampled from the ISMIP6 forcing files at the location of each marine-terminating glacier terminus and the calculated melt rate applied across each terminus face. Each group is free to simulate calving using any approach they see fit; the summation of melt and calving then provides the ice sheet model with total frontal ice ablation. This method may require substantial model development and a high model resolution of the outlet glaciers, but, whereas the retreat implementation imposes the same retreat amount for all glaciers within each sector, the submarine melt implementation allows for adjacent glaciers to retreat at a different rate, since the melt rate is computed for each individual marine-terminating outlet glacier. Ice sheet models that already have a different submarine melt parameterization are encouraged to implement the Rignot et al. (2016) parameterization as well. The difference between their existing melt parameterization (open experiment) and the proposed standard submarine melt implementation provides insight into the uncertainty in frontal melt parameterizations and the impact of the resulting ocean forcing on ice sheet evolution. At the same time, the possible wide range in submarine melt arising from different parameterizations obfuscates our understanding of how climate forcing uncertainty propagates into submarine melt uncertainty. Thus, the proposed implementation of a standard submarine melt parameterization allows us to examine how the spread in climate projections drives the spread in submarine melt rates and glacier retreat.

To illustrate the difference in the two ocean forcing approaches, we compare retreat against melt rates for two neighboring glaciers chosen at random and focus on the Tier 1 simulations under RCP2.6 and RCP8.5. In the retreat implementation, these two glaciers will retreat over the same distance in any given AOGCM and RCP scenario because all glaciers in this sector are forced with the same rate of terminus retreat (Fig. 9). However, these two glaciers exhibit very different projected melt rate anomalies, defined as the difference in mean melt rates between 2081–2100 and 1995–2014 (Fig. 10). King Oscar Glacier (glacier 51) will undergo a substantial but varied increase in melt rates in most RCP8.5 projections (Fig. 10a), whereas glacier 53 shows modest increases in melt rates in all projections (Fig. 10a). Thus, the ocean forcing in the melt implementation will be stronger for King Oscar Glacier than glacier 53. The melt implementation allows for greater spatial variability compared to the retreat implementation. The melt implementation also allows for greater variability for a given glacier as its terminus retreats and feedback with bedrock topography, e.g., higher melt rates where the bedrock is deeper (Fig. 10b–e). For a given glacier, the spatial patterns of projected melt rate anomalies are similar between the different CMIP models, but the magnitude differs as a result of different thermal ocean forcing and subglacial discharge. For the retreat implementation, the choice of a CMIP model impacts the timing and location of ice retreat (Fig. 9b–e), with NorESM1-M resulting in earlier retreat for a given location, and the retreat reaches further inland by 2100 compared to the other CMIP models. Nonetheless, both the retreat implementation and melt implementation show consistencies with each other (Figs. 9a and 10a): higher ocean forcing is projected for NorESM1-M, followed by MIROC5 under RCP8.5. HadGEM2-ES projects half of the retreat rate and melt rate at the end of the 21st century compared to NorESM1-M for this region and climate scenario.

Figure 11 Time series of ice shelf collapse area under the RCP8.5 and SSP5-8.5 scenarios over the 21st century (a) and corresponding ice shelf collapse masks' spatial evolution (b–g). The CMIP5 models are CCSM4 (blue, b), HadGEM2-ES (green, c), CSIRO-Mk3.6.0 (green, d), IPSL-CM5A-MR (cyan, e), MIROC-ESM-CHEM (red, f), and NorESM1-M (pink, g). The masks' spatial extent shown originates from the 4 km datasets and corresponds to the years 2015, 2030, 2050, 2070, and 2090.

6 Antarctic ice shelf fracture

The 2002 collapse of the Larsen B ice shelf was attributed to enhanced surface melting (Sergienko and MacAyeal, 2005; van den Broeke, 2005). Enhanced surface melt and water ponding can then trigger hydrofracturing and ice shelf collapse (Vaughan and Doake, 1996; Scambos et al., 2000, 2009). Other mechanisms associated with ice shelf collapse include rheological weakening, ocean waves, and surface load changes (MacAyeal et al., 2003; Braun and Humbert, 2009; Borstad et al., 2012). It has been hypothesized that ice cliffs resulting from ice shelf collapse could themselves become inherently unstable, a process called marine ice cliff instability (Bassis and Walker, 2011; DeConto and Pollard, 2016). As the processes for ice shelf collapse remain poorly understood and are rarely implemented in continental ice sheet models, ISMIP6 focuses on ice shelf collapse due to enhanced surface melt and provides time-varying masks for ice shelf fracture. The objective is to investigate the impact of ice shelf collapse, and the resulting loss of buttressing, by performing similar experiments without and with ice shelf collapse. The datasets were prepared following the method described in Trusel et al. (2015), which derives an annual surface melt from CMIP near-surface temperatures:

$$\begin{array}{}\text{(11)}& M_{\mathrm{AOGCM}}\left(x,y,t\right)=\mathrm{1183}\times e^{\left(\mathrm{0.4557}\times \mathrm{T}\mathrm{2}{\mathrm{m}}_{\mathrm{AOGCM},\mathrm{adjusted}}\left(x,y,t\right)\right)},\end{array}$$

where $M_{\mathrm{AOGCM}}\left(x,y,t\right)$ is the annual CMIP-derived surface melt flux (mm w.e. yr⁻¹) and T2m${}_{\mathrm{AOGCM},\mathrm{adjusted}}\left(x,y,t\right)$ is the downscaled, bias-corrected, CMIP near-surface temperature. Ice shelves are assumed to collapse following a 10-consecutive-year period with annual melt above 775 mm w.e. yr⁻¹, a threshold suggested by Trusel et al. (2015). ISMIP6 provides annual masks of ice shelf collapse covering 1995 to 2100 for CCSM4, CSIRO-Mk3.6.0, HadGEM2-ES, IPSL-CM5A-MR, MIROC-ESM-CHEM, and NorESM1-M under RCP8.5 (Fig. 11). Ice shelf collapse masks were created on the 4 km ISMIP6 Antarctica grid, and these were conservatively interpolated to generate masks at 2, 8, 16, and 32 km resolutions. Modelers should use the grid that is the most appropriate for their models. In the 4 km dataset, the mask values are set to 1.0 when the ice shelf is prone to collapse and to 0.0 for no collapse. The interpolation to coarser dataset may result in fractional mask values, with a value between 0.0 and 1.0 indicating partial collapse. Ice sheet models should remove the ice shelf when the mask is set to values greater than 0.0. If flagged, collapse is assumed to occur on 1 January of each year. A modeled ice shelf extent may not always exactly correspond to an ice shelf in the collapse mask dataset. In the event where the ice sheet model considers the ice to be grounded, but the ice shelf collapse mask indicates ice shelf removal, the collapse mask should not be imposed. Application of the mask may also result in ice shelf regions that are now detached from the ice shelf. In this case, these “icebergs” should be removed as well. Finally, a collapsed ice shelf should not regrow. The ice flow response to collapsed ice shelf is left at the discretion of modeling groups.

As indicated in Table 2, the Tier 1 model for this experiment is CCSM4. This model was prioritized over the two other models selected as part of the core experiments (MIROC-ESM-CHEM and NorESM1-M) because of its largest projected ice shelf collapse area under RCP8.5 (Fig. 11a). In CCSM4, ice shelf collapse is initially sporadic, with a collapse in the mid-2020s that lasts until the early 2030s. The 2030s and 2040s do not experience collapse. The second period of ice shelf collapse begins in the mid-2050s, lasts approximately 3 decades, and impacts ice shelves on both sides of the peninsula. The last, more rapid, phase of collapse begins in the early 2080s and impacts ice shelves in all regions around Antarctica. The spatial map of ice shelf collapse (Fig. 11b) indicates that the early collapse occurs predominantly at the fringes of the ice shelves in the peninsula, and in particular over the Wilkins ice shelf. The Larsen ice shelves and the George VI ice shelf collapse during the second period and are gone in 2070. The final phase of collapse, depicted by the conditions in 2090, includes the Abbot ice shelf and fringes of the Getz ice shelf, in West Antarctica, but also occurs at the terminus of small ice shelves throughout the East Antarctic. The large Ronne–Filchner and Ross ice shelves are not projected to experience ice shelf collapse, nor are the ice shelves fed by Pine Island and Thwaites glaciers. This “step-wise” increase over the 21st century in the projected ice shelf area is seen in all the datasets; however the timing, duration, and area affected differ between the models. The projected ice shelf collapse from HadGEM2-ES from 2080 to 2100 closely follows that from CCSM4 in terms of total area affected, but HadGEM2-ES projects a smaller area and delayed collapse for ice shelves on both sides of the peninsula and the Abbot ice shelf (Fig. 11a, c), comparable timing and extent for the Getz ice shelf, and a larger number of East Antarctic ice shelves at the end of the century. The projected behaviors for MIROC-ESM-CHEM and CSIRO-Mk3.6.0 converge in the period 2080 to 2100, despite MIROC-ESM-CHEM trailing behind initially. MIROC-ESM-CHEM projects some collapse over the fringes of the East Antarctica ice shelves towards the end of the century, while CSIRO-Mk3.6.0 does not, explaining the greater total ice shelf area collapse for MIROC-ESM-CHEM at the end of the century. In contrast, both NorESM1-M and IPSL-CM5A-MR only project ice shelf collapse in the peninsula (Fig. 11d, f), with IPSL-CM5A-MR being the model with the smallest area affected (Fig. 11a). As none of the selected CMIP5 models project ice shelf collapse over the Pine Island and Thwaites ice shelves, we do not expect this suite of experiments to trigger a collapse of the West Antarctic ice sheet. The ice shelf loss is concentrated in the peninsula for all models, but the timing and areas affected differ between the them. Therefore, the Antarctic-wide projected sea level change from this implementation of ice-shelf-triggered hydrofracture is anticipated to be small.

7 Discussion and conclusion

The protocols presented in this paper differ in detail between Greenland and Antarctica because the key processes by which ice is lost are different for each ice sheet. Certain aspects of the protocol are also influenced by the time limitations faced by ISMIP6 (due to the delay in CMIP6 simulations) and computational constraints limiting the use of RCM to Greenland only. Both atmospheric and oceanic forcings are required for both ice sheets, presented to the ice sheet model as an evolving flux of ice entering or leaving the ice sheet, or as fraction of ice that should be removed. Examples include the balance between snow accumulation and meltwater runoff on the ice sheet's surface (SMB), melt and/or refreezing on the underside of floating ice shelves, melt from the vertical faces on marine-terminating outlet glaciers, and the calving of icebergs (both from marine-terminating outlet glaciers and in the collapse of floating ice shelves). Several factors determine the choice of how this mass flux is calculated. These include the scale over which the mass flux varies spatially compared to the native resolution of AOGCMs. For instance, steep topographic gradients at the edge of the ice sheet introduce sharp spatial gradients in SMB, or the need to resolve runoff over the narrow Antarctic ice shelves, that are crucial in determining ice sheet response but are typically not resolved by AOGCMs. Other factors influence the quality of the SMB derived from AOGCMs besides resolution, such as the choice of physical parameterizations (e.g., Palerme et al., 2017) or limitation in the processes included in snowpack models, for example. In some cases, the key mass fluxes may not be determined by the AOGCMs at all. This is true, for instance, of melt from the underside of ice shelves and from Greenland marine-terminating outlet glaciers, or the feedback of the ice sheet's evolving geometry on the mass fluxes themselves.

In general, one of three approaches is used to determine each mass flux. Firstly, in some cases it is possible to employ the mass flux determined within the AOGCM directly. Second, an RCM can be used to simulate the necessary mass fluxes. This approach is used to provide SMB forcing for Greenland whereby the MAR RCM is forced using lateral and surface boundary conditions from the relevant AOGCMs (see Fettweis et al., 2017). The SMB forcing from Antarctica was not chosen to be obtained from an RCM due in part to the considerable time and computational resources required with this approach. Finally, a parametrization can be used to create the relevant mass flux based on variations in the relevant climate variables supplied by AOGCMs. This approach is used in the case of iceberg calving and melt from the underside of Antarctic ice shelves and Greenland marine-terminating outlet glaciers. Note that in addition to the primary forcing obtained from the AOGCM (for instance, an index of subsurface ocean temperature for ice shelf melt) these parameterizations may also require ancillary information; for example, the impact of surface runoff on the frontal melt of marine-terminating outlet glaciers in Greenland. In some instances, climatic forcing cannot be imposed directly as a mass flux and instead is implemented as a removal of ice. This is the case for the retreat rate implementation of Greenland outlet glaciers and for the fracture of Antarctic ice shelves due to increased surface melt.

The ISMIP6 experimental protocol attempts to balance the need for a full exploitation of the various sources of uncertainty on sea level projections against the high computational costs of some of the models likely to particulate in the intercomparison. Four types of uncertainty are considered. The first is the choice of emission scenarios – we chose to focus primarily on the high-end RCP8.5 scenario in addition to a limited amount of work with the RCP2.6 to bracket that uncertainty associated with emissions. The second form of uncertainty is the choice of the AOGCM within the overall CMIP ensemble for a given emission scenario. Here, a great deal of care was taken in selecting AOGCMs that represent the range of climate projections within the ensemble, as well as simulating the present-day climate of the ice sheets adequately. Issues related to delays within the CMIP6 led to a focus on CMIP5 supplemented by available CMIP6 results. This has the advantage that existing work on understanding variations within the CMIP5 ensemble (Agosta et al., 2015) could be used to aid model selection. Please see Barthel et al. (2020) for a detailed description of the procedures employed. The third source of uncertainty is that associated with the ice sheet models themselves. This includes the physics used by individual models (structural uncertainty) and the values chosen for poorly constrained parameters within the model (parametric uncertainty). These types of model uncertainty can be difficult to explore, but there are several ways in which ISMIP6 aims to do this. Firstly, the standard experiments, described by the protocols presented here and required from each participant, allow for direct comparison between ISMs. These models can vary widely in ice flow physics and numerical techniques – and therefore comparing the outcome of the standard experiments from different models, when all other experimental choices are equal, helps quantify the structural uncertainty. Secondly, modeling groups are also encouraged to submit open experiments in which they can employ their own parameterizations and include physical processes which may not be included in the standard experiments (which were designed to be compatible with all the participating models). One form of uncertainty that ISMIP6 does not sample in a systematic manner is ice sheet model parameter uncertainty. This requires quantifying the spread in projections by constructing an ensemble of simulations with model parameters varied in a methodical way as it would be computationally prohibitive for many of the participating ISMs. Modeling groups are however encouraged to repeat the experiments with different versions of their models. The final source of uncertainty is linked to our current lack of understanding of ice–ocean interactions: some of the ISMIP6 forcings that are used to parameterize ice–ocean processes in the ice sheet models have been provided with a range of values, to reflect uncertainty in the associated parameterizations.

Experiments are separated into core (or Tier 1) experiments that every participant is expected to complete and optional targeted experiments (Tiers 2 and Tier 3). This approach allows participants with models with faster run times to explore particular sources of uncertainty in more detail, while allowing all models to be compared across a smaller set of core experiments. The experiments are also separated to allow comparison between projections in which only one type of uncertainty changed: for instance, holding parameter and CMIP model choice fixed, to allow the impact of emission scenarios to be isolated. CMIP model selection results in a range of future climate scenarios, with some CMIP models projecting a warmer atmosphere but colder ocean conditions, and vice versa. As illustrated in this paper, there are strong regional differences in the atmospheric and oceanic forcings, which will transfer to differences in ice flow response and projected sea level. The datasets presented in this paper were specifically prepared by ISMIP6 and are available for community use. The preparation of this ISMIP6 protocol for a stand-alone ice sheet model brought together multiple communities towards the goal of understanding the uncertainty in future sea level change. ISMIP6 simulations and scientific achievements are expected to support the WCRP Grand Science Challenges on “Melting Ice and Global Consequences” and “Regional Sea-level Change and Coastal Impacts” and contribute to the projected sea level change in the IPCC Sixth Assessment Report. We anticipate that our protocol will benefit sea level research through the design of a consistent set of ISMIP6 simulations which will provide a basis for projections, comparative analysis, and further targeted ice sheet modeling activities. Progress towards resolving the pressing scientific and practical problem of projecting the contribution to global sea level change from the ice sheets, and reducing its uncertainty, requires continued collaboration between multiple disciplines in the Earth system sciences. Towards this end, we hope that ISMIP6 will succeed in promoting the scientific study of ice sheets and climate as a coupled system.

Appendix A: ISMIP6 Tier 2 and Tier 3 experiments

The optional Tier 2 and Tier 3 ISMIP6 experiments are described in Tables A1–A4 for Antarctica and Tables A5–A6 for Greenland. These experiments complement the mandatory Tier 1 experiments. Groups are highly encouraged to perform Tier 2 experiments, which complete the CMIP5 models selected in Barthel et al. (2020) for ISMIP6. In contrast, the CMIP6 models were a selection of opportunity. Tier 3 experiments consist of experiments with atmosphere- or ocean-only forcing, or the remainder of the ice shelf hydrofracture experiments. Groups can choose to focus on one aspect of these experiments.

Table A1 Description of ISMIP6 Antarctica Tier 2 simulations.

Table A2 Description of ISMIP6 Antarctica Tier 3 simulations, which are performed with ocean-only (OO) or atmosphere-only (AO) forcing, as indicated in the CMIP AOGCM column. N/A indicates not applicable.

Table A3 Description of ISMIP6 Antarctica Tier 3 simulations to investigate impact of ocean forcing.

Table A4 Description of ISMIP6 Antarctica Tier 3 simulations to investigate impact of ice shelf fracture.

Table A5 Description of ISMIP6 Greenland Tier 2 simulations.

Table A6 Description of ISMIP6 Greenland Tier 3 simulations, which are performed with ocean-only (OO) or atmosphere-only (AO) forcing, as indicated in the CMIP AOGCM column.

Table A7 Description of ISMIP6 Greenland Tier 3 simulations to investigate impact of ocean forcing.

Appendix B: ISMIP6 grids and variable request

To facilitate intercomparison of model submissions, groups are requested to submit on the ISMIP6 Antarctica and Greenland regular grids at a resolution that is the closest to the modeled ice sheet grid. For Antarctica, the ISMIP6 grid is a polar stereographic projection, with a standard parallel at 71^∘ S and central meridian of 0^∘ W on datum WGS84. The lowest left corner is at $x=-\mathrm{3040}$ km and $y=-\mathrm{3040}$ km, while the upper right corner is at x=3040 km and y=3040 km. Acceptable resolutions are 32, 16, 8, 4, 2, or 1 km. Submissions will be stored on the submitted resolution for archiving and conservatively interpolated by ISMIP6 to the 8 km Antarctica grid for intercomparison. For Greenland, the ISMIP6 grid is a polar stereographic projection, with a standard parallel at 70^∘ N and central meridian of 45^∘ W on datum WGS84. The lowest left corner is at $x=-\mathrm{720}$ km and $y=-\mathrm{3450}$ km, while the upper right corner is at x=960 km and $y=-\mathrm{570}$ km. Acceptable resolutions are 20, 10, 5, 4, 2, or 1 km. Submissions will be stored on the submitted resolution for archiving and conservatively interpolated by ISMIP6 to the 5 km Greenland grid for intercomparison.

The ice sheet data request (Table B1) contains key characteristics needed to evaluate the ice sheet geometry and ice sheet flow. It also contains key ice-sheet-specific boundary conditions that may differ between models and a record of the forcing applied to the ice sheet model. To facilitate the analysis of the ice sheet contribution to sea level, a number of integrated measures (for example, ice sheet mass) are also requested. Two-dimensional state variables (ST) are requested as a yearly snapshot corresponding to the end of the year in a simulation for state variables (such as ice thickness) and as yearly average for flux variables (FL, such as surface mass balance). Fields such as surface mass balance flux should be what was applied as a boundary condition to the ice sheet model and may be different from the input forcing file.

Table B1 Data request for the dynamical ice sheet model submissions. These fields, if applicable to the model, are saved on the regular ISMIP6 ice sheet grid that is the closest to a model native grid and contain yearly output. Type indicates whether the variable is a state variable (ST) or a flux variable (SF).

Appendix C: Antarctic atmospheric forcing preparation

In general, all files were obtained from the CMIP distribution through the Earth System Grid Federation (ESGF). As indicated in Table C1, the first ensemble member for each model and CMIP experiment was selected. (The ensemble member is denoted in CMIP as “r1i1p1”, where “r” is a center-designated realization number, “i” is the initialization number, and “p” is the physics number. For CMIP6, a center-designated ensemble forcing number “f” is also used, but this mostly corresponds to the CMIP experiment.) The CESM2 datasets originated from the initial “MOAR” run and were provided directly to ISMIP6 prior to this dataset becoming available on the ESGF grid.

The primary atmospheric variables to be used are precipitation (pr), evaporation (evspsbl), runoff (mrro, mrros), and skin temperature (ts). As defined by the ISMIP6 protocol (Nowicki et al., 2016), surface mass balance (SMB) is the net of precipitation minus evaporation minus runoff. In general, these were taken from the ESGF atmosphere “Realm” of variables and are defined globally so as to accommodate an ISM grid extending beyond continental boundaries. In the CMIP5 ESM output, the runoff variable from the land surface Realm is often problematic over ice sheets. There are two runoff variables – “surface” and “total” – which are ambiguously defined. In many cases one or both of the variables incorporate a restoration to the ocean of the accumulated SMB for maintaining mass equilibrium in the absence of a dynamical ice sheet model. This is commonly referred to as a poor man's iceberg calving, or “frozen runoff”. Liquid runoff is currently a negligible term for continental-averaged Antarctic SMB but could conceivably become locally significant by 2100 on ice shelves and in coastal regions, particularly along the northern Antarctic Peninsula. Hence reasonable effort was made to incorporate a runoff variable into the forcing from the CMIP simulations where available.

The CMIP experiment output that was used is from the historical experiment, from the 21st Century Representative Concentration Pathways (RCPs) 2.6 and 8.5 from CMIP5, and Shared Socioeconomic Pathways (SSPs) 126 and 585 of CMIP6. A climatology of each variable is first constructed corresponding to the annual average over model years 1995–2014. In CMIP5, the historical experiment ends at model year 2004, so that the climatology is produced from both the historical and RCP simulation output. In CMIP6, the historical experiment ends at model year 2014, and the climatology is taken entirely from the historical simulation output. Using the climatology, annual anomalies are then computed for each variable over the period 1950–2100. The climatology and the anomalies are then regridded to an azimuthal equal-area grid designated for ISMIP6 with 8 km spacing. In general, the CMIP ESM native grid spacing at high latitudes is very fine zonally and relatively coarse meridionally, and this produces artifacts when applying conservative interpolation methods. Here we have used cubic spline interpolation in the transfer from native resolution to the 8 km grid and conservative interpolation from the 8 km grid to other ISMIP6 grids. The CMIP ESMs selected have a range in native latitudinal grid spacing from less that 1^∘ to more that 2.8^∘.

The SMB datasets are in units of kilograms per square meter per second water equivalent and need to be converted by users to meters per year ice equivalent via

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}\mathrm{aSMB}\left[\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}\right]& =\mathrm{aSMB}\left[\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}\right]\\ & \times \mathrm{31}\mathrm{556}\mathrm{926}\left[\mathrm{s}{\mathrm{yr}}^{-\mathrm{1}}\right]\\ & \times \left({\displaystyle \frac{\mathrm{1}}{\mathrm{1000}}}\right)\left[{\mathrm{m}}^{\mathrm{3}}{\mathrm{kg}}^{-\mathrm{1}}\right]\times \left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right),\end{array}\end{array}$$

where ρ_w and ρ_i are the densities of water and ice (typically 1000.0 and 917.0 kg m⁻³), respectively. The temperature datasets are provided in units of kelvin.

Table C1 CMIP models used to create surface mass balance and surface temperature forcing for ice sheet models. N/A indicates Not Applicable.

^a Runoff computed from daily files.

Appendix D: Implementation of Greenland tidewater glacier retreat parameterization in ice sheet models

We describe the method and implementation of tidewater glacier forcing used in large-scale ice sheet models with possibly relatively coarse resolution and with initial geometries that can differ from observations. The approach we are proposing is implementing a time-dependent set of retreat masks that define the maximum calving front position at any time during an experiment. Differences in initial ice sheet model geometry require that the retreat masks are calculated specifically for individual ice sheet models. The procedure is a mapping operation to translate the retreat, originally derived for the observed ice sheet, to the individual model geometry. Furthermore, coarse-resolution models have to consider a form of sub-grid implementation to reduce biases when the calving front retreats across grid cells of large horizontal extent. We assume in the following that the time-dependent ice sheet retreat around Greenland is known for groups of marine-terminating outlet glaciers in seven different regions, as described by Slater et al. (2019).

Figure D1 Distances from the ocean for all points identified as ice grounded on bed below sea level in potential contact with the ocean. The data are masked to areas less than 50 km from the nearest grounding line.

D1 Retreat masks

Here we first discuss the hypothetical case where an ice sheet model of very high spatial resolution has been initialized with ice front and grounding line positions in perfect agreement with observations. We will assume that the model grid (MG) is identical to a regular observational grid (OG), where the ice sheet geometry is defined (e.g., Morlighem et al., 2017). To determine the retreat masks we apply the following procedure.

• 1.

  Identify bed below sea level and in connection with the ocean on OG (Fig. D1).

  • a.

    Define mask of ice grounded below sea level.

  • b.

    Search all connected grid points starting at the marine margin.

• 2.

  Identify shortest distance from the ice front and grounding line for all points on OG identified in 1 along sections of bed below sea level (Fig. D1).

  • a.

    Use mask of connected points as defined in 1.

  • b.

    Define ocean mask.

  • c.

    Calculate distance to the nearest ocean point (2b) for all points in (2a).

  Define retreat masks by thresholding the distance map (2) for given retreat distances per region (not shown).

D2 Procedure for specific models

Because the modeled is generally different from the observed initial ice sheet mask, additional complications arise because the modeled and observed glacier fronts cannot be assumed to be closely corresponding. In the model the glacier fronts may lie further out or in and different glaciers may fall together. Furthermore, the model may not resolve individual outlet glaciers due to limited resolution of a coarse grid. We therefore perform the distance calculations on OG and determine connectivity according to the observed geometry, but based on the modeled ice mask. The procedure described in Section D1 is augmented with interpolation steps between MG and OG.

• Find the modeled ice front positions and interpolate to OG.

  • a.

    Define the mask of grounded ice on MG (threshold area fractions to get a binary mask).

  • b.

    Interpolate grounded ice mask to OG using binned regridding (see above).

• 1

  Identify bed below sea level and in connection with the ocean on OG for the modeled ice mask.

  • a.

    Use mask of ice grounded below sea level from 0.

  • b.

    Find all connected grid points starting at the modeled marine margin.

• 2

  Identify shortest distance from the ice front and grounding line for all points on OG identified in 1 along sections of bed below sea level.

  • a.

    Use mask of connected points as defined in 1a and b.

  • b.

    Define ocean mask.

  • c.

    Calculate distance to the nearest ocean point (2b) for all points in (2a).

• 3

  Remap distances found in 2c from OG to MG.

  • a.

    Use binned regridding with masked distances from 2.

• 4

  Identify grid points on MG intersected by the OG grid points in 1 and determine weights as area covered by OG points on MG.

  • a.

    Use binned regridding with mask from 1.

  • b.

    An additional weight calculation may be needed if mask in 0 contained partial cells.

Because of a general mismatch (in resolution) between OG and MG, with MG typically being coarser, we translate retreat on OG to partial thinning on grid MG according to the covered area. This sub-grid process is discussed in more detail next.

D3 Sub-grid implementation

The implementation of outlet glacier retreat in a coarse grid model requires a form of sub-grid process to take into account partial retreat. This is needed because the two end members (1 – retreat only full grid cells that are entirely ice free; 2 – retreat full grid cells already when becoming partially ice free) are under- and overestimating the retreat, respectively. This problem is illustrated in Figure D2a. The gray shading shows grid cells on the high-resolution OG that fall within the footprint of the coarser MG (orange shading).

Figure D2 Schematic of the sub-grid implementation to translate the retreat mask defined on a high-resolution observational grid (a) to the coarser resolution (b) of an ice sheet model.

The method we have tested is to translate partial retreat to partial thinning, relative to a reference thickness applied once a year. This gave comparable results for a test case of different grid resolutions from 5 to 20 km resolution (not shown). The limitation to apply the relative thinning once a year is required to avoid time step dependence for different models. The thinning relative to a reference thickness avoids a nonlinear thinning with time. Applying partial thinning without a reference thickness (once a year) has been shown to overestimate retreat, because the thickness is exponentially decreasing and approaching the upper end member of full retreat.

Figure D3 (a) Distance from the ocean for all ice sheet points identified as in potential contact with the ocean. (b) Ice thickness of all points connected to the ocean.

D4 Discussion

The method has been developed for tidewater glaciers, which are the predominant form of marine termination around Greenland. The few glaciers with floating ice tongues are not treated differently. Extending the framework for floating ice shelves would require prescribing grounding line positions, which is not possible in the present ice sheet models. For the few outlet glaciers with floating ice at the termini, the ice sheet response will likely be underestimated, because removing floating may be expected to be less effective to speed up glaciers upstream compared to directly removing ice at a calving front. Compared to the other uncertainties associated with this method, we considered this a minor effect. In our application, most of the outlet glaciers are expected to retreat in the future. However, the method allows for readvance of glaciers up to the initial ice mask. Nevertheless, readvance can only happen by the ice flow into formerly vacated grid cells, as the method does not “create” mass. A partial retreat mechanism has to be considered to avoid over- or underestimation of the retreat, in particular in coarse-resolution models. This sub-grid process is implemented as partial thinning.

Appendix E: Abbreviations

AOGCM:	Atmosphere–ocean general circulation model
CliC:	Climate and Cryosphere
CMIP5 or CMIP6:	Coupled Model Intercomparison Project Phase 5 or Phase 6
ESGF:	Earth System Grid Federation
ESM:	Earth system models
FL:	Flux variable
IMBIE2:	Ice sheet Mass Balance Inter-comparison Exercise
IPCC:	Intergovernmental Panel on Climate Change
ISMIP6:	Ice Sheet Model Intercomparison Project for CMIP6
ISM:	Ice sheet models
ISM–ESM:	Ice sheet models are fully coupled within Earth system models
MAR:	Modèle Atmosphérique Régionale
MeanAnt:	Calibrations based on observed mean sub-shelf basal melt over Antarctica
MEOP:	Marine Mammals Exploring the Oceans from Pole to Pole
MG:	Model grid
OG:	Observational grid
PIGL:	Pine Island grounding line calibration
RACMO2.3p2:	Regional Atmospheric Climate Model version 2.3p2
RCM:	Regional climate model
RCP:	Representative Concentration Pathways
SMB:	Surface mass balance
SSP:	Shared Socioeconomic Pathways
ST:	State variable
WCRP:	World Climate Research Programme
WGS84:	World Geodetic System 1984
WOA18:	World Ocean Atlas 2018