Coupling the U.K. Earth System Model to Dynamic Models of the Greenland and Antarctic Ice Sheets

Abstract

The physical interactions between ice sheets and the atmosphere and ocean around them are major factors in determining the state of the climate system, yet many current Earth System models omit them entirely or treat them very simply. In this work we describe how models of the Greenland and Antarctic ice sheets have been incorporated into the global U.K. Earth System model (UKESM1) via substantial technical developments with a two‐way coupling that passes fluxes of energy and water, and the topography of the ice sheet surface and ice shelf base, between the component models. File‐based coupling outside the running model executables is used throughout to pass information between the components, which we show is both physically appropriate and convenient within the UKESM1 structure. Ice sheet surface mass balance is computed in the land surface model using multi‐layer snowpacks in subgrid‐scale elevation ranges and compares well to the results of regional climate models. Ice shelf front discharge forms icebergs, which drift and melt in the ocean. Ice shelf basal mass balance is simulated using the full three‐dimensional ocean model representation of the circulation in ice‐shelf cavities. We show a range of example results, including from simulations with changes in ice sheet height and thickness of hundreds of meters, and changes in ice sheet grounding line and land‐terminating margin of many tens of kilometres, demonstrating that the coupled model is computationally stable when subject to significant changes in ice sheet geometry.

Introduction

Ice sheets play a major role in determining global- and regional-scale features of the climate and Earth System. They are interactive partners in climate processes encompassing surface radiation, atmospheric and ocean circulation, as well as influencing solid earth physics, gravitation and sea level. The growth of ice sheets is subject to several positive feedbacks with the climate - for example they commonly have a high albedo surface that reflects solar energy, promoting cooling and further growth as the ice expands in area. Large changes in ice sheet surface height also affect atmospheric circulation and precipitation patterns (Ridley et al., 2005), including those that lead to snowfall over the ice sheet itself. The supply of fresh water to the ocean from surface runoff, ice shelf basal melt, and calving of solid icebergs is likewise important for regional ocean circulation near the ice and the creation of bottom water in the global ocean circulation (e.g., Silvano et al., 2018), which distributes heat around the planet. The various interactions of ice sheets with the rest of the Earth System have helped amplify small changes in the solar radiation received by Earth due to variations in our orbit around the sun enough to cause the repeated glacial cycles of the last million years (e.g., Abe-Ouchi et al., 2013). Enough water to raise sea level by 65 meters is currently stored in the ice sheets of Antarctica and Greenland (Morlighem et al., 2017, 2019), but there has been a total variation of $\sim$140 m in sea level over these glacial cycles (Lambeck et al., 2014), implying a major role for ice sheets in radical changes in the climate and the location of habitable land on Earth. The recent review by Fyke et al. (2018) provides a detailed overview of the physical interactions between ice sheets and the rest of the Earth System. Global mean sea level (GMSL) rise is a significant consequence of the increasing concentrations of greenhouse gases in the atmosphere. Ice sheet mass loss accounts for around a third of the present rate of GMSL rise, and this contribution is expected to increase and eventually dominate GMSL change in the coming decades and centuries (Oppenheimer et al., 2019). Importantly, ice sheet mass change is also persistently the most uncertain term in the future GMSL budgets (Edwards et al., 2021; Meehl et al., 2007; Oppenheimer et al., 2019). It is thus important to understand the physical detail of how climate and ice sheets interact on sub-millennial timescales, and to be able to produce projections of the ice sheet contribution to GMSL rise that are consistent with the latest climate change projections used by the international community. The latter is the goal of the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6, Nowicki et al., 2016. CMIP6 is phase 6 of the Coupled Model Intercomparison Project, Eyring et al., 2016). The majority of the ISMIP6 projections are being carried out by ice sheet models (ISMs) run independently from the climate models used to create the climate projections, but the limitations of this uncoupled approach are widely recognised. For example, such studies omit the local feedbacks between ice sheet height and surface mass balance, the influence of changing ice sheet topography on large-scale atmospheric circulation, and struggle to provide physically justifiable basal melt forcing to floating ice shelves as they retreat (Jourdain et al., 2020; Slater et al., 2020). Such feedbacks can only be addressed and understood through simulations where the climate and ice components are allowed to evolve together.

Global models used for century-scale climate projections generally lack sufficient spatial resolution to explicitly reproduce the gradients in spatial patterns of accumulation and melt at the interfaces of the ice with the atmosphere and ocean, and have historically required that the locations of their solid boundaries remain fixed. For this reason it has been common to heavily parameterize and limit the ways in which ice sheets are included in Earth System models (ESMs). The focus of most coupled climate—ice sheet studies has been on the millennial-scale evolution of glacial paleoclimates, using ”intermediate complexity” ESMs (Claussen et al., 2002) that can be affordably integrated over relevant timescales (e.g., Ganopolski et al., 2010; Gregory et al., 2012; Goelzer et al., 2016), and similar models have been used to project future ice sheet mass loss and its implications for sea-level (e.g., DeConto & Pollard, 2016; Gregory et al., 2020; Huybrechts et al., 2002; Robinson et al., 2012). Empirical approximations are commonly used in such models to translate climate model fields into boundary conditions for ice sheets, such as the positive-degree-day method to relate surface air temperature to snow melt (Reeh, 1991) or transfer functions that use open-ocean conditions to derive melt rates for ice shelves without explicitly modeling any of the ocean domain local to the ice (Jourdain et al., 2020; Slater et al., 2020). It can be difficult to gain robust insight into the physics of climate—ice interactions when the components are coupled via proxy variables in this way, and these simpler ESMs also inevitably lack process detail in aspects of their general climate simulations compared to state-of-the-art climate models.

The additional computational expense of more complex climate models restricts the length of simulations - and thus the science questions - that could be addressed using interactive ice sheets in these models, and until recently even the most complex global models also did not have realistic enough simulations of polar surface physics to be directly coupled to ISMs. Ridley et al. (2005) was an early example of a simulation with a CMIP-class (at the time) climate model coupled to an ice sheet. However, improvements in both the spatial resolution and the process modeling in newer generations of climate models have led to an increase in work in this area, with, for example, Vizcaíno et al. (2013) using a subgrid-scale method to diagnose Greenland surface mass balance directly from the surface simulation of the Community Earth System Model (CESM, Hurrell et al., 2013). In the CMIP6-era, Muntjewerf et al. (2020) conducted ISMIP6 projection simulations using a dynamically evolving GrIS in CESM2 (Danabasoglu et al., 2020), and a number of groups are developing sophisticated coupling methods to enable the inclusion of ice sheets in their CMIP6 models (e.g., Gierz et al., 2020; Kapsch et al., 2021).

In this paper we describe the means by which models of the Greenland (GrIS) and Antarctic ice sheets (AIS) have been coupled to the UK Earth System Model (UKESM1, Sellar et al., 2019), producing a model suitable for investigating the physical feedbacks between modern ice sheets and the climate system consistent with its global climate projections. To calculate these feedbacks in a way that is closely driven by the physical processes in the climate model the ice sheets are simulated as interactive components of the ESM. They are directly coupled to the atmosphere—land—ocean general circulation model (AOGCM) via the exchange of fluxes of energy and mass, and information about the evolving shape of the ice sheet surfaces.

This work represents a very substantial model development and advancement in the scientific capabilities of UKESM1. The technical difficulties in evolving the shape of the solid boundaries in an AOGCM and keeping its climate simulation computationally stable are considerable, and no other CMIP6 ESM has yet been coupled to an interactive model of Antarctica. This paper will provide a technical description of our new capabilities, show that our coupled system evolves in a reasonable way and that our coupling techniques remain stable when subject to significant changes in ice sheet geometry. Our coupling methods are complex and deserve documentation and scrutiny in their own right to underpin future studies of ice sheet behavior and GMSL rise with this model. The results we show are thus intended primarily to demonstrate the function of our coupling methods, and detailed evaluation or analysis of the climates they represent is beyond the scope of this paper. Analysis of results of simulations using this model will be presented in forthcoming publications.

In Section 2, we will describe the basic model components and modifications to them relevant to ice sheet coupling in UKESM1. Section 3 presents an overview of the coupling framework, and Section 4 describes the land surface—ice coupling in detail. Section 5 describes how icebergs are handled in the coupled framework, and Section 6 details how ice shelves interact with the ocean component of UKESM1. We conclude with a discussion in Section 7.

Component Models: Descriptions and Modifications

The core atmosphere and ocean components of UKESM1 are derived from the GC3.1 configuration of the U.K. Met Office HadGEM3 AOGCM (Kuhlbrodt et al., 2018; Williams et al., 2018, GC3.1 hereafter). The elements of the ice sheet coupling system described in this paper can be used in either GC3.1 or UKESM1, depending on how they have been configured, and the example results shown in this paper have been drawn from a range of different setups. The details and naming conventions of the configurations that have been used are given at the end of this section after the components have been described.

As well as supplying the physical core of UKESM1 models, GC3.1 is an AOGCM used in its own right to contribute simulations to CMIP6. It comprises version GA7.1 of the Unified Model (UM) atmosphere, version GL7 of the Joint UK Land Environment Simulator (JULES) land surface model (Best et al., 2011), the GO6 configuration of version 3.6 Nucleus for European Ocean Modeling (NEMO) ocean model (Madec & Team, 2016) and version 5.1 of Community Ice CodE (CICE, Hunke et al., 2015). UKESM1 configures some of these core components slightly differently, and additionally includes components for land and ocean biology, and atmospheric and ocean chemistry. In this section, we note modifications we have made to these model codes and other relevant aspects that enable the coupling of interactive ice sheets to GC3.1 or UKESM1.

Core Physical Models UM - Atmosphere Model

The formulation of the UM atmosphere model itself is not modified to enable the ice coupling. The gravity-wave and land-surface drag schemes in the UM use measures of the subgrid-scale variation in surface topography, and these must be recalculated along with the topography when the ice sheet surface evolves. These fields are recalculated for the UM as part of the ice sheet coupling framework (see Section 4.1) after the ISM has run. In all the following model variants the UM is used with an “N96” horizontal resolution (1.$\mathrm{875}°$ longitude $\times$ 1.$\mathrm{25}°$ latitude) and with 85 vertical levels that cover the troposphere and stratosphere.

JULES - Land-Surface Model

The land surface model in GC3.1 provides the surface mass balance (SMB) to force the ice sheet model, and has been modified in several aspects from how it is used in GC3.1. It is important for ice sheet dynamics to model the geographical variation of SMB as accurately as possible, especially features such as the strong gradient at the margins of the Greenland ice sheet, from high-altitude net accumulation to net ablation at lower altitudes. The spatial resolution of the UM and JULES in UKESM1, like most global atmosphere models, is not sufficient to do this - at 70°N an N96 UM grid-box has sides ∼70 km long, whilst, for example, the recent ISMIP6 experimental protocol (Nowicki et al., 2020) chose to provide surface forcing data for ISMs on a 1 km grid - so some form of climate downscaling must be used to produce realistic ice sheet SMB.

JULES can represent subgrid-scale heterogeneity in surface characteristics through a system of tiles of different surface types. Each tile calculates its own exchange fluxes with the atmospheric boundary layer, with the average state of each grid-box being constructed for the atmosphere through an area-weighted sum of the tiles present. Building on this tile paradigm, Shannon et al. (2019) describes how JULES tiles can be extended to represent surfaces at different elevations within a grid-box and uses this capability to model subgrid-scale surface exchange on topographically steep glaciers in JULES, via a simple downscaling of the grid-box mean climate to each such “elevated tile”. GC3.1 does not use elevated tiles, but UKESM1 does, and we will here expand on the brief description given in Sellar et al. (2019) to motivate and demonstrate their use. Our use of elevated subgrid-scale tiles to model SMB on an ice sheet is similar to the approach used to model the surface of Greenland in CESM (Sellevold et al., 2019; Vizcaíno et al., 2013).

UKESM1 uses elevated tiles to divide each ice sheet grid-box into 10 elevation ranges, with finer resolution at lower elevations where temperatures are more likely to lead to surface melting. We use the same range division as CESM1 (Vizcaíno et al., 2013), with elevation class boundaries at 0, 200, 400, 700, 1000, 1300, 1600, 2000, 2500, 3000 and 10,000 m. The same class boundaries are used for Greenland and Antarctica. Coupled climate—ice configurations of UKESM1 have two tiles allocated to each elevation class in a grid-box: one for ice (e-ice) and one for rock (e-rock), with separately modeled snowpacks, and different values specified for the nominal depth and effective heat capacity of the subsurface. The area fraction of each of these tile types in an elevation class is set so as to represent the glaciated area at that height in the higher resolution ISM representation of the ice sheet. The e-ice and e-rock fractions in every JULES grid-box that has them must sum to 1. In coupled climate-ice configurations of UKESM1, e-ice tiles are, by construction, always covered with a snowpack of significant depth (here we use 100 m although the precise figure is somewhat arbitrary, see Section 4.1). Snow is not always present on e-rock tiles, and a nominal albedo of 0.1 for rock is assumed for the subsurface if it becomes exposed.

Lower atmosphere air temperature and humidity are downscaled for elevated tiles following moist or dry adiabats as appropriate, as described in the JULES documentation for tiles with offset heights (Best et al., 2011). As in Shannon et al. (2019) we also adjust the downwelling longwave radiation seen by each tile consistent with the adjustment to the atmospheric temperature. In UKESM1 we further scale the resulting longwave radiation on tiles to preserve the original grid-box mean flux. In UKESM1 we do not downscale shortwave radiation, surface winds or precipitation components to the elevated tiles. These quantities have more complex relations to topography that cannot be adequately represented by dependence on altitude alone. Moreover, adjusting the phase of precipitation in JULES would make it inconsistent with atmospheric energy and moisture fluxes. Shannon et al. (2019) found that observations of accumulation and turbulent fluxes on certain individual glaciers could only be closely matched in JULES with considerable empirical tuning of precipitation and surface wind forcing, but such a fine-grained approach is not appropriate for a global-scale model that may be used to model climates and ice sheets that are quite different from present.

Each JULES tile uses a column of a multi-layer snowpack model. In GC3.1 the snowpack is represented by three thin surface layers sufficient to model seasonal snow, whereas in UKESM1 a number of modifications are made to make it adequate for modeling perennial snow and ice sheet SMB. The snowpack scheme conserves the mass of snow, with the depth being determined by the prognostic density of the snow. In UKESM1 the snow scheme is equipped with 10 levels, with layer interfaces at 0.04, 0.16, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0 and 6.0 m below the snow surface, providing resolution for the internal temperature of the snow and for percolation and refreezing of surface meltwater. The bottom level of the JULES snowpack expands in thickness to hold as much mass as the column accumulates, although in our implementation of ice sheet coupling the snowpack cannot grow indefinitely since mass in regions experiencing positive mass balance will be passed to the ISM (see Section 4.1). The SMB field we pass to the ISM as a surface boundary condition is simply the change in the total mass of the snowpack on the JULES e-ice tiles between coupling intervals. Our snowpack mass changes due to the accumulation of meteroic snow and losses from sublimation at the surface and the runoff of melted snow that has not refrozen inside the snowpack. However much the snowpack accumulates, the depth of snow in JULES does not directly affect the orographic height of the surface seen by the UM.

Refreezing and percolation of meltwater are modified with respect to their behavior in GC3.1, where snowpack pore space (the unoccupied space between snow grains that can accommodate liquid water) is purely a function of depth below the surface. In UKESM1, once compaction of the snowpack increases the density of a layer beyond a specified threshold (450 kg/m³), the pore space of the layer is gradually reduced as density increases further. A layer that achieves the density approaching that of solid ice in the snowpack at 850 kg/m³ is impermeable and thus forms an “ice lens”, diverting any percolating meltwater that reaches it from the snowpack above into runoff from the ice sheet. Meltwater that has percolated to the bottom of the snowpack without refreezing also runs off, as does all rain falling on the ice sheet. Rain in UKESM1 is not currently allowed to interact with the snowpack and plays no part in the SMB calculation. Not allowing rain to refreeze in the snowpack means that we may overestimate the potential for firn pore space to buffer meltwater losses in ice sheet ablation zones. Runoff is routed to outflow locations at the coast, predetermined by modern surface slopes (Oki & Sud, 1998), and passed to the adjacent surface ocean. Routing of runoff is not recalculated to take account of changes in ice sheet topography, nor does it contribute to subglacial hydrology in the ISM, which we do not model in UKESM1.

Radiation at the surface, and thus also the albedo in UKESM1, are modeled with a two-stream approximation where the shortwave is split into separate spectral bands for visible and near-infrared light. The treatment of snow and ice albedo in UKESM1 is also modified with respect to that in GC3.1, in which the snow albedo is always a function of grain size (Best et al., 2011). In UKESM1, on ice elevated tiles where surface snow density exceeds a threshold at 550 kg/m³, the albedo is instead a function of snow density, and declines toward low values typically observed for bare ice in these spectral bands as the snow ages or its interstitial pore space fills with meltwater. The snow scheme does not allow for deposited pollutants to affect the albedo, nor for meltwater to pond at the surface, although these phenomena are targets for future development.

The bottom boundary condition under the snowpack on each elevated tile is provided by a subsurface that has a heat capacity and provides a temperature boundary condition for the snowpack above. Unlike the soil subsurface of other JULES grid-boxes, the elevated tile subsurface is impervious to water and has no carbon content. For e-ice tiles, the subsurface represents the solid interior of the ice sheet, and surface melt may only percolate through and refreeze in the snowpack above. The subsurface is configured differently for e-ice and e-rock tiles, with different values specified for the nominal depth, effective heat capacity and the albedo should the subsurface be exposed to the atmosphere. When coupled to the ISM, our coupling implementation means that e-ice tiles always have mass in their snowpack and their subsurface can never be exposed. In areas of strong ablation in an ice-coupled configuration, ice-density snow will be relayered to the top of the e-ice snowpacks and this is effectively seen as bare ice in the climate model. The heat flux from the snowpack may also be passed as a surface boundary condition to the ISM, which then calculates a temperature for use as the bottom boundary condition in JULES. Our current model configurations, however, (including the examples shown in this paper) use a fixed-temperature assumption in the ISM and this coupling is effectively inactive.

Like the simple ice surface tile that has existed in all previous versions of JULES, the e-ice and e-rock tiles we have described cannot coexist in any given grid-box with the other JULES surface tile types, which represent various kinds of vegetation, soil, lakes or urban surfaces. A grid-box's e-ice and e-rock fractions can change, but it cannot switch to using normal JULES tiles and the JULES soil subsurface model. This means that the ice sheet cannot spread beyond pre-defined areas where grid-boxes with elevated tiles have been configured, and vegetation cannot recolonize areas from which ice has retreated. Although reestablishment of vegetation on bare ground occurs in only a few decades (Cannone et al., 2008; Tisdale et al., 1966), the millennial timescales of ice sheet dynamics involved in revealing significant areas of bare rock are too long for simulations with GC3.1 or UKESM1 for reasons of computational expense.

It is technically very complex to make the land-sea mask of the UM depend on time, so UKESM1 cannot currently represent inundation or exposure of land due to changing GMSL. One consequence of this, and the fact that JULES treats the surface of floating ice shelves as “land” (although we may in fact model ocean circulation underneath this surface, as in Section 6), is that the extent of Antarctic ice shelves cannot change in UKESM1. We will describe how this condition is practically enforced in Sections 2.1.4 and 6.1. All Greenland's outlet glaciers are treated as land-based glaciers in UKESM1, including those that are marine terminating in reality, as they are too small to be resolved on the ocean grid. As such, all Greenland glaciers may retreat from the coast in UKESM1 to expose unglaciated land, but they cannot advance past the coastline in JULES or past their initial calving front in the ISM. Our inability to represent true ice shelf retreat, advance or collapse is a significant limitation of our current scheme, but it does mean that JULES only needs to alter e-ice and e-rock tile fractions to represent the movement of the ice sheets that our system allows.

NEMO - Ocean Model

NEMO (Nucleus for European Modeling of the Ocean) provides the ocean component of both GC3.1 and UKESM1. Full technical details can be found in Madec and Team (2016), and Storkey et al. (2018) describe the GO6 version our configurations are based on. In this section, we describe specific aspects relevant to ice sheet coupling.

In GO6, NEMO uses a Lagrangian iceberg tracking scheme, fully described in Marsh et al. (2015). In response to a given calving flux in an ocean grid-box, icebergs are generated according to a climatological size distribution of length, width and volume. The icebergs drift according to the surface ocean state while their mass is steadily reduced (an iceberg cannot gain mass in this scheme) by parameterisations of basal melting, buoyant convection, and wave erosion. These icebergs cannot ground on shallow topography. In GC3.1 or UKESM1, when not coupled to an ice sheet, the spatial distribution of the calving flux is fixed although the magnitude can be allowed to vary according to the overall surface mass balance to maintain global mean ocean salinity if desired (Williams et al., 2018). Without coupling to an ice sheet, these models have no other way of returning snow that has fallen in the accumulation zone of an ice sheet to the ocean. With the ice sheet model coupled in, the calving flux is specified by the ISM directly in both spatial distribution and magnitude.

When the interactive AIS model is omitted, the ocean is inactive wherever there is an ice shelf and NEMO has a solid, vertical wall to the ocean floor at the ice shelf front, as is standard for GO6. In reality, there are ocean “cavities” with circulation beneath ice shelves and melting and freezing at the ice shelf base. When the dynamical model of AIS is used with UKESM1, we use the GO7 configuration of NEMO (Storkey et al., 2018), in which the active ocean can have an upper boundary of thick floating ice, represented at each location by a set of vertically contiguous inactive ocean layers including the top layer. GO7 normally uses a prescribed climatology of basal melting under ice shelves, but as in Mathiot et al. (2017) we allow UKESM1 to interactively calculate melting and freezing within the active ocean cavity using a 3 equation method (Holland & Jenkins, 1999) with a fixed depth boundary layer (Losch, 2008).

Two different mesh resolutions are commonly used in GC3.1 and UKESM1, based on the tripolar ORCA grid (Madec & Imbard, 1996) but with the definition extended under Antarctica to allow ocean to replace retreating ice: eORCA1 (nominally $1°$ longitude/latitude) and eORCA025 (nominally $\frac{1}{4}°$ longitude/latitude). The fjords that Greenland's marine-terminating glaciers flow into are all too small to be represented on either grid, and for this reason the ocean does not currently provide a marine boundary condition to GrIS in our coupling scheme. In the Antarctic seas and under floating ice shelves the eORCA025 grid has grid-boxes with sides between 4 and 14 km in length and eORCA1 grid has grid-boxes with sides between 16 and 58 km in length. Although neither of these is sufficient to explicitly resolve eddies or fine topographic features important in determining ocean dynamics close to Antarctica, these grids do nevertheless allow for the majority of the ice shelves around Antarctica to be represented in our ocean domain. Mathiot et al. (2017) shows that the simulation of basal mass balance under Antarctic shelves in a regional configuration of NEMO is acceptable at eORCA025 resolution.

Both GO6 and GO7 use a z* vertical coordinate (Adcroft & Campin, 2004), in which the vertical resolution and placement of the grid-box tops with respect to the depth of the ice—ocean interface (the ice shelf base) varies significantly with location. At the ocean floor and the ice—ocean interface the grid allows partial-height cells, each of which may be assigned any height within 20%–100% of its default thickness so that the vertical location of the boundaries can be more precisely represented. All UKESM1 configurations use 75 vertical layers in NEMO, with default thicknesses that increase from 1m near the surface to 203 m at depth. At shallower grounding lines (the boundary between ice that is floating and ice resting on bedrock) around 500 m deep the vertical resolution is 50 m, at 1000 m deep the resolution is 100 m, and by 1500 m vertical resolution is 130 m.

To enable ice shelf-ocean coupling, the ability to evolve the position of the ice—ocean interface and smoothly adapt the ocean state to match the new boundary locations has been added to NEMO. Changing the model's solid boundaries and the associated ocean state in a self-consistent way is complicated by NEMO's staggered Arakawa C-grid discretization for tracers and dynamics. Our techniques for doing this are described in detail in Section 6.

UniCiCles - Ice Sheet Model

The core of the ISM used in UKESM1 configurations with ice sheet dynamics is the Berkeley Ice Sheet Initiative for Climate at Extreme Scales model (BISICLES, Cornford et al., 2013). UniCiCles (Unified Model-CISM-BISICLES) is a package combining BISICLES with an interface that obtains boundary conditions from Unified Model or JULES data, using code derived from the Glint interface of the Glimmer-CISM ISM (Lipscomb et al., 2019; Rutt et al., 2009). We regard this interface code as part of the surface coupling, and describe its functions in Section 4.

BISICLES uses the adaptive-mesh Chombo libraries (Adams et al., 2019). Adaptive meshing allows areas of the ice sheet where the accurate representation of the dynamics is critically reliant on fine resolution in small areas (e.g., at grounding lines and shear margins of fast moving glaciers) to be modeled without incurring prohibitive computational expense across the slower-moving body of the ice sheet. All cells in the mesh are rectangles that may be recursively refined by subdivision into four smaller cells with the same aspect ratio, should accurate resolution of the dynamics warrant this.

Our configuration of BISICLES approximates the momentum equations using the ”shelfy-stream” approximation with simplified vertical shear strains included in the effective viscosity (Schoof & Hindmarsh, 2010), often referred to as SSA*. Basal traction is set to zero beneath floating ice and modeled using power laws beneath grounded ice. GrIS uses a linear drag law everywhere while AIS uses a cubic law far upstream from the grounding line, which tends to a Coulomb friction law near the grounding line (Tsai et al., 2015). BISICLES deals with partly floating cells primarily through its AMR, but also through three minor approximation choices: (a) one sided differences in the gravitational driving stress improve the numerical treatment of the stress balance at moderate resolutions; (b) oceanic melt rates are set to zero in partly floating cells to prevent unphysical grounding line retreat that otherwise occurs proportional to the mesh spacing; (c) basal traction can be reduced in partly floating cells. No adjustment is made to the bedrock or gravitational attraction anywhere in UKESM1 in response to changes in the ice. Ice in UniCiCles is assumed to have a fixed density of 910 kg/m³.

Our current configurations use fixed 3D internal temperature fields, and do not use interactive basal hydrology or subglacial outflow. We assume that geothermal heat flux beneath grounded ice is negligible. The uncertainty in the temperature field and temperature dependence of rate factor are dealt with by locally softening or stiffening the ice, using a coefficient of the effective viscosity known as the stiffening factor. Initial values for the stiffness factor, along with the basal traction coefficient, are inferred from observations of surface velocity and geometry (Cornford et al., 2015; Lee et al., 2015). This tuning method provides parameters that vary in the horizontal plane and compensates for some limitations of the model dynamics. The parameters have no dependence on time, although for timescales longer than a few hundred years atmospheric warmth penetrating through the ice by advection and diffusion would have a significant impact on flow rates (Huybrechts, 1996). For such runs the internal temperature evolution would need to be considered, as well as potential solid earth responses to changes in ice mass. Our coupling scheme is technically capable of passing a surface heat flux from the JULES snowpacks into UniCiCles, but not into the base of floating ice shelves from NEMO.

We use the version tagged as BISICLES-UKESM-ISMIP6, with no notable technical additions to the model used in Cornford et al. (2013). Specific settings for the GrIS and AIS domains, including the treatment of basal friction and grounding lines around Antarctica, follow the configurations used in ISMIP6 initMIP exercise (Goelzer et al., 2018; Seroussi et al., 2019). In our GrIS configuration, basal melting beneath floating ice is set to zero and the lateral margins are not constrained over land and are free to spread over unglaciated parts of Greenland. To allow for the lateral spread of ice, the basal drag coefficients determined for the grounded ice during initialization are spread to neighboring land areas.

In the example simulations used in this paper, GrIS is modeled with 9.6 km square base cells that may subdivide to 1.2 km and AIS with 8 km that may subdivide to 2 km (Figure 1). Although these meshes were computationally efficient during model development, further levels of refinement should ideally be used on AIS to more correctly model sensitive grounding line evolution (Vieli, 2005). The meshes are updated every 8 timesteps for GrIS and 4 for AIS allowing the resolution to evolve with the ice dynamics. BISICLES uses an explicit timestep scheme where the timestep decreases with increasing resolution. The timestep is set to 0.25 times the minimum grid spacing divided by the maximum ice speed. The timestep is reduced when necessary to produce output for diagnostics or coupling if requested at specific times.

Figure 1. Example Berkeley Ice Sheet Initiative for Climate at Extreme Scales mesh refinement and ice speeds on (a) Greenland and (b) Antarctica taken from the example UKESM1.0-ice simulation in Section 6.3.2.

A fixed-front calving model is used to maintain the position of marine terminating ice (see Section 5) for both GrIS and AIS. AIS is forced to maintain a minimum ice thickness of 10 m, which means that ice is artificially added to grid cells where necessary. This violates conservation of mass, but ensures that no new open ocean areas are created in our system. As noted in Section 2.1.2, we cannot allow ice dynamics to change the land-sea mask in UKESM1. For GrIS there is no need for a minimum floating ice thickness because fjords are not resolved on the NEMO or JULES grids, and they treat all GrIS glaciers as land-based. GrIS ice in UKESM1 is thus free to retreat inland, and in UniCiCles GrIS shelves can thin and disappear if the SMB or upstream flux of ice are low enough, but it cannot advance beyond its initial calving front position.

Components of GC3.1 and Different UKESM1 Configurations

UKESM1 and GC3.1 share a common codebase, meaning that all of the features and modifications described in Section 2.1 are available to all configurations. Here we describe which ice sheet features are actually used in each configuration (see Table 1. Henceforth we will use these specific UKESM variant names in the text, and use “UKESM1” to refer to topics that apply equally to all variants). Some of these configurations use our new features simply to generate additional information about the climate of glaciated regions and do not couple to a dynamic ISM. In both UKESM1.0 and GC3.1 configurations the initial snow mass in each grid-box sets the maximum that can be melted. Eventually this snow mass is exhausted in ablation areas, whereafter ice sheet runoff is unrealistically limited. This problem is not encountered in models with a dynamic ice sheet component.

Table 1 Summary of Ice Modeling Capabilities of Different AOGCM and ESM Configurations Referred to in This Paper

	Ice sheet snow physics	Elevated tile SMB	Sub-shelf ocean, shelf basal melt	Dynamic GrIS	Dynamic AIS	Notes
HadGEM3-GC3.1	-	-	-	-	-	Williams et al. (2018); Kuhlbrodt et al. (2018)
UKESM1.0	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	-	-	-	Sellar et al. (2019)
UKESM1.ice	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	-	[IMAGE OMITTED. SEE PDF.]	-	no UKCA, TRIFFID or MEDUSA
UKESM1.0-ice	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	no MEDUSA

GC3.1 models (Kuhlbrodt et al., 2018; Williams et al., 2018), across all the available resolutions, have none of the ice sheet modifications described in this section. They use the 3-layer seasonal snow configuration of the JULES snow scheme without elevated tiles, and have no ocean cavities under Antarctic ice shelves. They do not include a dynamic ISM.

UKESM1.0 (Sellar et al., 2019), as used for many of the UK CMIP6 simulations, is commonly run at N96 resolution for the atmosphere and eORCA1 in the ocean, the lowest of the supported horizontal resolutions for current UM-based AOGCMs (see Williams et al., 2018 for more details of the grids). This configuration has most of the land surface modifications described in Section 2.1.2, using elevated tiles with climate downscaling and the 10-layer snow configuration with modified albedo and percolation for both Greenland and Antarctica. This configuration does not use e-rock tiles, so cannot represent fractional glaciated area in a grid-box, and the ice sheet area and surface altitude are specified as in GC3.1, in which the GrIS is too extensive. This configuration uses a fixed iceberg calving flux, and has no ocean cavities under Antarctic ice shelves.

UKESM1.ice is used for the coupled climate—ice-sheet simulations in the UK submission to ISMIP6 (Nowicki et al., 2016). It has a coupled, dynamic representation of the Greenland ice sheet, but not Antarctica. UKESM1.ice is based on GC3.1, not UKESM1.0, because the latter had not been finalized when our ISMIP6 simulations were started. The ISMIP6 runs and examples in this paper used an N96, eORCA1 resolution. Unlike UKESM1.0, this configuration disables the interactive biology and chemistry modules in the atmosphere (UKCA StratTrop1.0, Archibald et al., 2020), land (TRIFFID, Cox, 2001) and ocean (MEDUSA-2.0, Yool et al., 2013), and limits vegetation in the land surface model to 5 plant functional types. It differs from GC3.1 in that ice sheet regions are equipped with both e-ice and e-rock tiles with climate downscaling and the additional snow physics described in Section 2.1.2 for both Greenland and Antarctica. This configuration does not have ocean cavities nor an interactive AIS. The UniCiCles ISM (see Section 2.1.4) simulates the evolution of the Greenland ice sheet and the calving flux for northern hemisphere icebergs. In the ISMIP6 simulations available through the Earth System Grid Federation, the calving flux magnitude is constant at a value that balances the simulated GrIS SMB for the 1970s, with a time-dependent pattern determined by UniCiCles.

UKESM1.0-ice is a configuration based on UKESM1.0, used in this paper but not for ISMIP6. This configuration includes UKCA, TRIFFID with 13 plant functional types (although vegetation cannot recolonize areas from which ice has retreated), e-ice and e-rock tiles and all the modifications in Section 2.1.2 for ice sheet surface fluxes. Because MEDUSA is currently technically incompatible with ocean cavities, ocean biogeochemistry is not included. UniCiCles simulates the dynamics of both GrIS and AIS, updating the orography in the ESM and controlling the iceberg calving flux distribution and magnitude. Different examples in this paper use either an eORCA1 or eORCA25 resolution ocean, both with an N96 atmosphere.

For each ice sheet domain, UniCiCles simulates 1 year in $\frac{1}{3}$—1 hour (depending on the degree of refinement required) on 36 cores. This is a computationally inexpensive part of UKESM1, which takes 6 hr on 720 cores to simulate 1 year of climate at N96eORCA1 resolution or 16 hr on 1,440 cores for 1 year of climate at N96eORCA025 resolution.

Coupling Framework

Of all the Earth System components in UKESM1, the atmosphere dynamics require the shortest timestep, at 20 min. The other sub-models, excluding UniCiCles, communicate once in every few of these timesteps, on timescales from minutes to hours. Structurally, this communication is achieved either by having the model subroutines compiled into the same executable (e.g., UKCA, JULES, TRIFFID and the UM are compiled together, and MEDUSA and CICE are part of the NEMO executable) or using the OASIS coupling software (Craig et al., 2017) to pass information between separate executables (Figure 2). On our computing platforms UKESM1 is run in individually submitted, sequential jobs simulating between one and three months each, and in each such job both OASIS and subroutine coupling pathways are used many times. These jobs, and coordination of many other related tasks (e.g., postprocessing and archiving of output) are automatically organized by Cylc (Oliver et al., 2019), a flexible scheduler inherent to the UKESM1 and GC3.1 infrastructure.

Figure 2. Coupling structure of the components of UKESM1.0-ice. Boxes joined by broken circles indicate that component routines are compiled together into one executable. Red arrows indicate hourly exchange between components running in parallel, gray arrows are annual exchange between components running serially.

Although ice sheet surface coupling with the atmosphere (at either the bare ice or firn surfaces) and ice shelf basal coupling with the ocean require as frequent communication as atmosphere surface exchanges do elsewhere, the dynamic flow of the ice sheet can be adequately modeled with much longer timesteps. Commonly used flow approximations are stable with timesteps of several months, and it is common to entirely neglect sub-annual variation in boundary conditions for ISMs. There is thus a clear need to separate the treatments of physical processes in and around the ice sheet that operate on these different timescales. In UKESM1 we incorporate the faster varying surface processes directly into other components of the ESM and reserve UniCiCles for ice sheet dynamics and thickness evolution, communicating with the ISM much less frequently and reducing the computational overhead of coupling. The resulting relative infrequency of coupling to UniCiCles means that many of the individually submitted UKESM1 jobs that simulate only a few months of climate do not need to communicate with UniCiCles at all.

Aside from the infrequency of ice communication, OASIS is also not suitable for coupling UniCiCles into UKESM1 because of our need for the coupling to carry out some specialized operations, described in the following sections. Among these is updating the geometry of the representation of the ice sheet in JULES and NEMO, as only UniCiCles is formulated to be able to move the ice sheet while running. Because of these considerations, we couple UniCiCles to the rest of UKESM1 by using Cylc to coordinate the exchange of boundary condition files in between the individual UKESM1 jobs in the sequence. The standard diagnostic systems in each component are used to accumulate and output files of time averaged fluxes or instantaneous locations of boundaries (see Table 2). Cylc runs scripts to turn these outputs into input files for the other components, which are read upon initialization of the next job when boundaries can be changed. The relative infrequency of the ISM coupling means that the computational overhead of this rather brute-force approach is negligible. All our ice-coupled UKESM1 configurations pass information from the atmosphere and ocean components to UniCiCles once per year of climate model simulation, in northern hemisphere winter. The time of year is an arbitrary choice, and no significant difference has been seen from more frequent coupling.

Table 2 Fields Passed Between Climate Models and Ice

Variable	Averaging	Domain
Land surface $\to$ ice
SMB	Annual average	e-ice tiles in JULES ice sheet grid-boxes, UniCiCles ice
Snow mass	Instantaneous	e-rock tiles in JULES ice sheet grid-boxes, UniCiCles unglaciated
Tile area fraction	Instantaneous	e-ice tiles in JULES ice sheet grid-boxes
Sub-snowpack heat flux	Annual average	e-ice tiles in JULES ice sheet grid-boxes, UniCiCles ice
Ice $\to$ land surface
Area of tile as fraction of grid-box	Instantaneous	Tiles in JULES ice sheet grid-boxes
Glaciated fraction of each tile	Instantaneous	Tiles in JULES ice sheet grid-boxes
Sub-snowpack temperature	Annual average	UniCiCles ice upper surface, e-ice tiles in JULES ice sheet grid-boxes
Changes to non-glaciated snowpack	Annual total	e-rock tiles in JULES ice sheet grid-boxes
Ice $\to$ atmosphere, via CAP
Surface height	Instantaneous	Mean over UM grid-box
Standard deviation of orography	Instantaneous	Mean over UM grid-box
Peak-to-trough orographic height	Instantaneous	Mean over UM grid-box
XX, XY, YY gradients		Mean over UM grid-box
Ocean $\to$ ice
Ice shelf BMB	Annual average	On ice shelf lower surface
Ice shelf basal heat flux	Annual average	On ice shelf lower surface
Ice $\to$ ocean
Column height of water under ice shelf	Instantaneous	Grid-boxes around Antarctica
Calved ice mass	Annual mean	Grid-boxes at surface marine boundary of ice

A useful advantage of having the ISM run as a quasi-standalone model for each period in our framework is that, if desired, the ice can be evolved for more than one year for every year of the climate boundary conditions calculated by UKESM1. Assuming that (a) the change in ice sheet geometry over $N$ years is too small to have a significant effect on the climate and SMB, and (b) that the global climate changes negligibly within $N$ years in a given simulation, we can allow our system to integrate the ISM by $N$ years for every 1 year of climate simulated by simply recycling the climate boundary condition files. It is clearly most accurate to use $N$ = 1, and evolve the climate and ice sheet synchronously, and this is how we usually simulate interactive ice with UKESM1. Using $N>1$ to accelerate the evolution of the ice with respect to the climate is a useful way to reduce computational expense, especially when investigating scientific questions about the far future or paleoclimatic changes with an expensive model such as UKESM1. An example of this asynchronous coupling for the GrIS is given at the end of Section 4.2.1 in which we use $N$ = 5, basing our choice of $N$ on Gregory et al. (2020) who use a similar scheme in a simpler coupled model. The use of $N>1$ distorts the evolution of the system when rapid changes in climate are occurring, especially if those changes occur as part of feedback loops with ice sheet dynamics. For example, when changes in ocean circulation are caused by addition of ice sheet meltwater, using asynchronous coupling means that ocean response will be too slow in ice sheet years, and water will not be conserved. The accumulation of meteoroic snow into a sufficiently large snowpack to initialize new glaciated areas in UniCiCles is also difficult to handle under such an asynchronous scheme. On multi-centennial timescales other assumptions used in our UKESM1 configurations, such as the fixed internal temperature of the ice sheet or lack of solid earth response, also become harder to justify and would need to be considered when interpreting the simulation results.

UniCiCles runs separate domains for each ice sheet it simulates, for example a relatively small regional Greenland domain at the same time as a much larger Antarctic one for a modern climate. This means it is simple to run simulations with just one ice sheet interactive or have both active at once. Domain-specific boundary files for each ice sheet are prepared by the coupling scripts as input to the ISM, and the outputs from each ISM domain are merged into global fields for input to the climate model. The BISICLES adaptive mesh itself does not need any special treatment, since the BISICLES interface accepts boundary information on regular grids whose cells coincide with any one of the sizes of grid-box allowed by mesh in use.

More details of the aims and implementation of each coupling task are given in the following sections, which are organized according to the different communication pathways used by UniCiCles: with the atmosphere above the ice sheet surface (Section 4); with icebergs at the surface of the ocean (Section 5); underneath floating ice shelves (Section 6).

Surface Climate and Mass Balance Implementation of Coupling

The surface boundary conditions that UniCiCles requires from the climate model for the ice sheet dynamics are the SMB and, if the ice sheet internal temperature were configured to evolve, a surface heat flux. Snowpack properties, SMB and its components - accumulation, sublimation, runoff of melted snow - and heat fluxes are diagnosed by JULES in ice sheet grid-boxes on every e-ice tile in the domain, including those with zero area (for such tiles this is a potential SMB) that may be needed during interpolation to the ISM mesh. Tiles with non-zero area are used to calculate grid-box-mean surface heat fluxes to the overlying atmosphere and surface runoff that is routed to the ocean. Fields required for coupling to UniCiCles (listed in Table 2) are output from the UM in diagnostic files as time-means over the ice coupling period. When Cylc determines that the ice sheet dynamics should be run, submission of the next segment of the climate run is delayed until the ISM segment has completed, and the first of the coupling scripts is triggered. This script selects the required variables from the UM and JULES output and reformats them so they may be read by UniCiCles.

The other steps that precede BISICLES dynamic time-stepping are all carried out inside UniCiCles itself, separately for each ice sheet domain. Snow and ice mass is conserved in the UM, whilst UniCiCles deals in ice volume, so in all cases the fixed ice density of 910  kg/m³ specified by BISICLES is used to convert between mass and volume representations. First, the SMB from the e-ice tiles and snow mass from the e-rock tiles are trilinearly interpolated from the three-dimensional latitude-longitude-elevated-tile fields provided by the JULES grid to the ice sheet surface on the regular Cartesian grid of UniCiCles, using the highest spatial resolution permitted by the adaptive mesh. The trilinear interpolation is not inherently conservative, so (optional) corrections have been implemented to adjust the variables on the ISM grid so that they either match the JULES originals when integrated over the whole ice sheet domain, or on the scale of individual JULES grid-boxes. In the examples in this paper these options are not used as they were found to sometimes introduce grid artifacts into the coupling fields that were less acceptable than the small amount of non-conservation they were correcting. Non-conservation may also arise when the extent of the ice sheet must be represented differently on the grid of each model, for example due to mismatches of coastline, but in those cases an adjustment to the SMB to conserve the area-integral might produce unrealistic results. Any unglaciated UniCiCles grid-boxes where the interpolated snow mass from e-rock tiles exceeds a threshold (we use 100 m ice-equivalent snow mass) are converted to glaciated grid-boxes with the corresponding depth of ice. Long-term accumulation of snow away from the ice sheet can thus result in the glaciation of new areas. On the other hand, the SMB from e-ice tiles is applied only to already glaciated UniCiCles grid-boxes, so this term cannot directly increase the ice sheet extent. With the boundary conditions prepared, BISICLES then evolves the ice sheet dynamics inside UniCiCles. It is assumed unnecessary to repeat three-dimensional interpolation of the UM boundary conditions to the ISM topography for every BISICLES timestep and the interpolation is only done once per coupling period.

After the ice dynamics have finished, the final subroutines of UniCiCles remap a number of fields from the ISM grid back to the elevated tiles of the UM (see Table 2). Glaciated UniCiCles grid-boxes where the ice thickness falls below a certain threshold (we use 50 m) are deglaciated, and their ice converted to a solid, ice-density snow mass that has the temperature of the ice sheet surface that is added to the corresponding e-rock tiles in JULES. This deglaciation threshold is imposed to ensure that a realistic negative SMB from an e-ice tile in JULES could not ablate all the ice in a glaciated UniCiCles grid-box within one coupling period, even under significantly warmer climates. The 50 m we use for this threshold is less than the 100 m used to initialize a new ice grid-box due to snow accumulation on e-rock tiles. Having different thresholds for glaciation and deglaciation introduces some hysteresis into the advance/retreat behavior in our model. In practice we found that if the two thresholds have the same value then some marginal areas can switch back and forth between glaciated and unglaciated representations every few coupling timesteps, which is undesirable as transferring mass/volume between the JULES tiles and higher resolution UniCiCles grid is not a perfectly reversible operation. The choice of these exact threshold values is somewhat arbitrary. To diagnose an accurate SMB we need sufficient snow on an e-ice tile to ensure that it cannot all melt in a year, and since we calculate ”potential” SMB on all e-ice tiles in the domain this no-melt criteria must be satisfied in some areas that would see very high melt rates if they existed in reality. A newly glaciated area with 100 m of ice does not cause an extreme perturbation to the flow field in BISICLES and still allows for a significant difference to the lower deglaciation threshold at 50 m, ensuring that JULES does not melt more ice than exists in UniCiCles at any location.

After UniCiCles has finished, the information required to reinitialize UM and JULES with the new ice sheet state is derived. From the highest-resolution grid of the final BISICLES diagnostic output a UM-format version of the regional topography is created, which is used to calculate the orographic height and statistics of subgrid-scale topographic variation on the JULES grid. The topography statistics are calculated by the program used to make standard UM boundary condition files, to ensure methodological consistency. Finally, the area fractions and snowpack properties of the elevated tiles in JULES are modified. The area fractions for the tiles reflect the horizontal and vertical distribution of glaciated and unglaciated grid-boxes in UniCiCles, remapped to the grid-boxes and elevation classes of the UM grid.

Snowpack properties on the JULES elevated tiles are modified in three ways. First, by adjusting the thickness of the bottom layer, the snowpack mass on each e-ice tile is reset to its initial value (equivalent to 100 m of ice-density snow). This step balances the mass changes specified by the SMB field that were applied to UniCiCles. Second, since JULES deals with snow mass per unit area, we account for the fact that changes in tile area resulting from the evolving shape of the ice sheet alter the area-integral mass of snow lying on the JULES tiles. Where e-rock tiles have grown in area, their snowpack mass in kg/m² is reduced to keep the area-integral in kg unchanged, and the snowpack layering is rearranged in line with this new depth. Where e-rock tile areas have shrunk, their mass in kg/m² is kept constant but the snow that has been lost is instead transferred to the ice sheet as an increment to the next year's SMB on the e-ice tile at the same elevation. Third, if any UniCiCles grid-boxes were glaciated following accumulation of snow, this mass is removed from the snowpack of the corresponding e-rock tiles in JULES. The changes described to the orography, tile areas and snowpack are applied to the restart file of UM and JULES. Once this completed the coupling cycle is finished and Cylc triggers the next segment of the climate simulation.

Example Results

This section provides technical illustrations of the surface climate—ice-sheet coupling methods in UKESM1. A more complete scientific evaluation of UKESM1 climate experiments with GrIS and AIS will be published in forthcoming work. These illustrative results for both ice sheets come from a UKESM1.ice (see Section 2.2) simulation forced according to the specification for a CMIP6 historical simulation (Eyring et al., 2016) in which the GrIS was coupled from 1970 onward, but the AIS was not coupled. The ice sheet topographies in the UM and the elevation tile fractions in JULES are taken from the initial UniCiCles states, which were derived as described in Section 2.1.4. Snowpacks on JULES e-ice tiles were initialized with 100 m of ice-density snow for compatibility with our SMB coupling method (see Section 4.1), e-rock tiles start without snow. Other aspects of the global ocean, land and atmosphere are taken from 1970 in instance r1i1p1f1 of the GC3.1 CMIP6 historical simulation (Andrews et al., 2020).

Greenland

UKESM1.ice surface climate results are shown averaged over 1970–1999, during which the ice sheet shape does not evolve significantly. They are compared with a simulation of the Modèle Atmosphérique Régional (MAR) regional climate model (RCM) configured for Greenland (Fettweis et al., 2017) over the same period that was forced at the boundaries with output from one of the members of the UKESM1.0 CMIP6 historical ensemble (Hofer et al., 2020). For presentation purposes, the figures show UKESM1.ice output from the matrix of JULES elevated tiles interpolated onto the same orography and grid as used in MAR.

The simulated SMB on e-ice tiles in UKESM1.ice compares well with the pattern and magnitude of SMB simulated by MAR (Figure 3). UKESM1.ice simulates an area-integrated total of 396 Gt/yr for the period 1970–1999 (ice sheet area 1.61 × 10¹²m²), compared to 515 Gt/yr in MAR (ice sheet area 1.78 × 10¹²m²). The spatial gradients of extreme positive and negative SMB at the margins of the ice sheet are not as sharp as in MAR, and the accumulation maxima on the south east coast are underestimated, likely due to UKESM1 being unable to fully resolve orographically forced precipitation. Integrated over the area of the ice sheet, both accumulation and runoff are slightly underestimated by UKESM1.ice at 525 Gt/yr and 127 Gt/yr respectively, compared to 655 Gt/yr and 163 Gt/yr in MAR.

Figure 3. 1970–1999 time mean GrIS SMB (m/yr LWE). (a) UKESM1.ice (interpolated from e-ice tiles onto the Modèle Atmosphérique Régional (MAR) orography, minimum = −1.36 m/yr; maximum = 1.28 m/yr), (b) MAR (minimum = −3.71 m/yr; maximum = 3.33 m/yr) forced by the UKESM1.0 historical climate simulation. The contour of SMB = 0 is marked in purple (UKESM1.ice) and green (MAR).

Retention of meltwater by refreezing in the snowpack is likely to be a significant influence on future mass loss of the GrIS (Noël et al., 2017). UKESM1.ice simulates a very similar pattern of refreezing to MAR (Figure 4). The proportion of meltwater that is refrozen in the highest melt area of GrIS in the south is higher in UKESM1.ice than in MAR, both on the east and west coasts. There is also refreezing in UKESM1.ice, but not in MAR, around the northern coasts. The magnitude of refreezing in the south in UKESM1.ice is supported by specialized snow models and observations (Steger et al., 2017), but the simulation of refreezing in the north is less likely to be realistic. There is more energy available for melt in the UKESM1.ice simulation than in MAR, especially higher up in the accumulation zone, which would account for some of the differences in melt and refreezing between the models. The positive feedback between melt, albedo and absorption of shortwave radiation makes determining a single cause for the differences challenging, and the horizontal resolution of JULES means that the representation of the available pore-space in the snowpack in UKESM1.ice is unlikely to capture realistic regional detail.

Figure 4. 1970–1999 time mean Greenland refreezing (m/yr LWE). (a) UKESM1.ice (interpolated from e-ice tiles onto the Modèle Atmosphérique Régional (MAR) orography, maximum = 0.95 m/yr), (b) MAR (maximum = 1.01 m/yr) forced by the UKESM1.0 historical climate simulation.

The surface albedo simulated for GrIS is an important determinant of summer melting. UKESM1.ice summer albedo is slightly lower than MAR over most of the central part of the ice sheet (Figure 5), and the lowest values are not so tightly confined to the margins. This lower albedo is reflected in the larger ablation area that can be seen in Figures 3 and 4. It is not simple to attribute these differences in the simulation of albedo to any one factor. The basic structures of the snow albedo schemes in JULES and MAR are similar but there are differences in the detail of formulation and tuning. As already noted, there is more energy available for surface melt on GrIS in UKESM1.ice than MAR, but this includes a significant component from turbulent heat fluxes as well as radiative fluxes. While the UKESM1.ice summer albedo is lower than MAR, it is not unrealistic compared to observations. The simulated 1982–2014 mean summer albedo averaged over the tile fractions compares well with the mean taken over the same mask from the CLARA-A2 (CM SAF cLoud, Albedo and Radiation, second edition) dataset based on satellite observations (Riihelä et al., 2019), with values of 0.77 and 0.74 respectively.

Figure 5. 1970–1999 time mean GrIS JJA albedo. (a) UKESM1.ice (interpolated from e-ice tiles onto the Modèle Atmosphérique Régional (MAR) orography, minimum = 0.58; maximum = 0.81), (b) MAR (minimum = 0.52; maximum = 0.83) forced by the UKESM1.0 historical climate simulation.

To demonstrate the ability of our coupled model to cope with large changes in the extent and shape of the GrIS we carried out a multi-centennial experiment with atmospheric carbon dioxide concentration instantaneously quadrupled in 1970. The ice sheet and climate are coupled annually until , and thereafter 5 years of ice dynamics are run for every 1 year of climate (N = 5, Section 3). After 1600 years of ice evolution the GrIS has lost around two-thirds of its area (Figure 6) with only a remnant dome on the east of Greenland. As well as providing a major change to the local surface conditions, this degree of ice mass loss also represents a significant perturbation to the orographic forcing of the atmospheric circulation, including impacts on the path and strength of the Atlantic jet stream. We are thus confident that our technical approach to coupled climate and ice sheet dynamics creates plausible interactions between model components and is robust for our intended uses.

Figure 6. Extent and height of the Greenland ice sheet after 1600 years under an idealized 4xCO2 climate. (a) Orographic height of glaciated grid-boxes on the UniCiCles grid (maximum = 3,265 m). The blue contour shows the original extent of the ice sheet as represented in UniCiCles, the gray contour shows the coastline; (b) The corresponding change in orographic height seen by the atmospheric dynamics on the UM grid (minimum = −2,548 m, maximum = 89 m). The blue contour shows the original extent of the ice sheet as represented in the UM, the gray contour shows the coastline.

Antarctica

Since we do not have an RCM simulation for Antarctica forced with UKESM1 climate boundary conditions, we compare the AIS SMB from this same historical UKESM1.ice simulation with MAR forced by ERA-interim (Dee et al., 2011) reanalysis output (Agosta et al., 2019). Biases in the background climate of UKESM1 thus influence this comparison. In particular all UKESM1 (and GC3.1) simulations have significant warm biases in the Southern Ocean and around the Antarctic coast.

The year-round cold surface temperatures on the AIS mean that SMB is overwhelmingly positive and generally quite small. Melting is rare and sublimation is more important for surface ablation (Agosta et al., 2019). Future climate warming over the AIS is likely to change this situation (Barthel et al., 2020), and the future SMB for this ice sheet is very uncertain. UKESM1.ice simulations have the same broad spatial pattern as MAR of greater accumulation near the coasts, especially in West Antarctica (Figure 7). UKESM1.ice simulates an area integrated total of 2,040 Gt/yr for the period 1970–1999 in this simulation, compared to 2,517 Gt/yr in MAR and 2,516 Gt/yr in RACMO2.3p2, another RCM forced by ERA-interim (van Wessem et al., 2018). UKESM1 simulates less snowfall over the whole of the ice sheet, and especially misses small-scale accumulation maxima associated with topographic features in MAR. Accumulation in Antarctica is controlled more by atmospheric circulation and less by local altitude than in Greenland, so our elevated tiles without any rearrangement of precipitation would be expected to be less successful for downscaling climate forcing to the AIS. Additionally, katabatic winds and blowing snow, which are important influences on the patterns of surface accumulation on Antarctica, are not well represented in UKESM1 due to the grid resolution.

Figure 7. Antarctic surface mass balance (m/yr LWE). (a) UKESM1.ice 1970–1999 time mean (interpolated from e-ice tiles onto Modèle Atmosphérique Régional (MAR) orography, minimum = −0.44 m/yr, maximum = 1.33 m/yr), (b) MAR forced by ERA-Interim reanalysis (Dee et al., 2011) 1979–1999 time mean (minimum = 0.0 m/yr, maximum = 2.41 m/yr).

MAR simulates almost no runoff because the small amount of meltwater all refreezes in the snowpack (Agosta et al., 2019). UKESM1.ice shows areas with a small net negative SMB because of greater melting, notably around the Amery ice shelf (Figure 7). Low-lying e-ice tiles on the Antarctic Peninsula, and those representing the Fimbul, Lazarev, Amery and Shackleton ice shelves in East Antarctica simulate both melting ($<$1 m/yr) and runoff ($<$20 cm/yr) (Figure 8). Almost everywhere else around the coast there is some melting, but it is smaller ($<$10 cm/yr) and refreezes, resulting in negligible runoff. Integrated over the whole AIS, UKESM1.ice produces considerably more melting and runoff than MAR (600 Gt/yr vs. 40 Gt/yr of melting, 100 Gt/yr vs. 1 Gt/yr runoff), likely because of the warm bias in the UKESM1 Southern Ocean. RACMO2.3p2 simulates 71 Gt/yr of surface melt and 3 Gt/yr of runoff.

Figure 8. Antarctic surface melt and run off (m/yr LWE) in UKESM1.ice 1970–1999 time mean (interpolated from e-ice tiles onto Modèle Atmosphérique Régional orography). (a) surface melt (maximum = 1.40 m/yr), (b) runoff (maximum = 0.63 m/yr).

Calving and Icebergs

Several different calving physics parameterisations are available in BISICLES to generate icebergs from the body of the ice sheet at marine-terminating boundaries. In UKESM1, where the extent of NEMO's floating ice shelves must be fixed because the atmosphere-ocean land-sea mask cannot evolve (see Section 2.1.4), we use a simple fixed-front calving scheme in which any ice that flows across the initial marine boundary of the ice sheet is treated as solid discharge to the ocean. If the ice thickness thins to zero before the flow reaches the boundary, or the ice speed at the boundary is zero, then no ice is calved. As noted in Section 2.1.4, the UKESM1 BISICLES configuration for AIS is required to maintain a minimum ice thickness of 10 m to maintain these ice shelves whereas there is no minimum thickness for GrIS ice. The case where calving ceases due to ice thickness becoming 0 thus only occurs for grounded ice flowing toward a calving front at the coast.

The flux of calved ice from UniCiCles is used to seed icebergs that are then handled by the Langrangian tracking scheme in NEMO, which allows them to drift and melt in the ocean (see Section 2.1.3). Our coupling simply converts the BISICLES diagnostic of calved mass into an iceberg seed boundary condition file for NEMO.

Implementation of Coupling

Each period of UniCiCles simulation results in a two-dimensional distribution of the mass of discharged floating ice beyond the ice margin, from which a coupling script creates two boundary condition files. One file is used in the next execution of NEMO as an annual-mean source of icebergs, which is applied spread over the next year of climate simulation. The resolution of the coastline is far more detailed in UniCiCles than in NEMO, and the regridding assigns each UniCiCles grid-box to the nearest open-ocean grid-box of NEMO, which in the case of Greenland can be a considerable distance away. The other file specifies a rate of ice mass removal on the UniCiCles grid, which is applied across the whole of the next execution of UniCiCles, removing the ice that has been used for icebergs in NEMO.

Example Results

Our UKESM1.ice simulation produces smaller iceberg fluxes per grid-box, from many more coastal locations of Greenland, than the prescribed flux used in UKESM1.0 historical simulations (Figure 9). The integrated discharge of 343 Gt/yr is smaller than 430 $\pm$ 40 Gt/yr estimated from observations for this period by Mankoff et al. (2019), but is greater than the 229 Gt/yr specified in UKESM1.0 simulations without interactive ice sheets, an amount that was calculated to minimize global salinity drift in the preindustrial simulation.

Figure 9. Iceberg seed magnitude and distribution (Gt/yr). Left: calving from BISICLES in UKESM1.ice 1970–1999 time mean (maximum in a 9.6 km UniCiCles grid-box = 27.34 Gt/yr); right top: BISICLES calving mapped to the NEMO grid to seed icebergs in UKESM1.ice 1970–1999 time mean (maximum in an eORCA1 grid-box = 52.8 Gt/yr); right bottom: time mean climatology of the Greenland iceberg seed specified in UKESM1.0 (maximum in an eORCA1 grid-box = 70.0 Gt/yr). The Greenland integral calving for each representation is given in each panel.

The difference between UniCiCles and NEMO coastlines is less important for Antarctica. In our simulations icebergs are calved almost everywhere along the marine boundary of the ice sheet, although the magnitude of calving depends strongly on the initialization of the AIS model and the state of the ocean under the shelves (see Section 6.3.2).

Ice Shelf Basal Mass Balance

As discussed in Section 2.1.3, even the higher resolution eORCA025 ocean mesh available to UKESM1 is too coarse to resolve the marine-terminating outlet glaciers of the GrIS, so basal mass balance (BMB, the net effect of ablation and accretion of basal ice) is simulated only for the AIS. At both eORCA1 and eORCA025 resolutions many Antarctic ice shelves are represented, along with a sub-ice-shelf ocean circulation enabling the calculation of basal melt. Smaller ice shelves such as Pine Island Glacier are poorly resolved on the eORCA1 grid, limiting the degree of ocean coupling available in UKESM1 with this configuration.

Implementation

The ocean circulation in the cavity beneath an ice shelf depends on the geometry of the cavity and on density-driven effects of ice shelf basal melting and freezing. The ocean circulation can fluctuate much more quickly than the cavity can change shape due to ice dynamics. Therefore the calculations of the melting and freezing at the ice—ocean interface are done within NEMO on the same timestep as the ocean dynamics (15 min for eORCA025, 45 min for eORCA1). This is similar to the reason for ice sheet snowpack processes being simulated in our land surface model rather than UniCiCles.

For coupling NEMO to UniCiCles, ice shelf BMB is calculated in NEMO as described in Mathiot et al. (2017). NEMO outputs the time-mean mass flux over the ice coupling period as a two-dimensional field along the ice—ocean interface and this is passed to UniCiCles. We do not yet pass temperatures or heat fluxes to the ice sheet model to use internally as part of our basal coupling (UniCiCles is run with fixed internal ice temperatures), although NEMO's parameterization of BMB does include an idealized diffusive heat flux at the base of the ice that is specified according to the thickness of the shelf and the temperature of the water.

After the climate segment has completed, a script bilinearly interpolates the BMB flux from the NEMO grid to the coarsest (8 km) mesh used in the UniCiCles AIS configuration and stores it in a BISICLES boundary file. As with the SMB, this bilinear interpolation is not conservative (conservation in the ice ocean coupling will be discussed in Section 6.2). In our configurations the ocean couples to UniCICles once every year of climate simulation. No significant difference in results has been seen from coupling more frequently. We have not tried BMB coupling in the framework of asynchronous coupling ($N>1$, Section 3), so cannot comment on that method in this context.

For coupling UniCiCles to NEMO, aside from the iceberg calving flux (see Section 5) information about the geometry of the floating ice shelves is the only other thing we pass. Discretising the new shelf geometry in NEMO is a complex procedure with a number of stages, partly done by scripts and partly by new routines inside NEMO, which we describe in the following. The scheme has proven to be computationally stable to a wide range of realistic and idealized perturbations to ice shelves that we have tested.

We diagnose the height of the ice shelf base above the sea floor bathymetry in UniCiCles, and regrid this field to represent the water column thickness (the depth of water between the bottom of the ice and the sea floor) in NEMO. A depth for the base of the ice shelf in NEMO can reliably be derived from this quantity without concern that any differences in the representation of bathymetry between UniCiCles and NEMO will affect the shape of the ocean cavity we model. A script then identifies the BISICLES grid-boxes on its highest resolution mesh (2 km) that overlap each NEMO grid-box. If at least 50% of them contain grounded ice (rather than floating ice or open water), the NEMO grid-box is deactivated, so that it becomes grounded ice. Otherwise, the NEMO water column thickness is made equal to the mean water column thickness from the BISICLES grid-boxes. NEMO has no means of calculating melt for a subgrid-scale distribution of floating and grounded areas. Some areas in UniCiCles near the grounding line may thus have no melt calculated for their location if the NEMO box they correspond to has been grounded, and in areas of steeply sloping or highly variable bathymetry the use of a simple grid-box mean water column thickness in NEMO may result in melt values that are not appropriate for every UniCiCles point they are applied to. If we use a lower threshold than 50% of grounded BISICLES grid-boxes as the criteria for ungrounding NEMO grid-boxes, the ocean cavity spreads under a greater area of ice and the area-integrated basal melting is correspondingly larger, but we see little difference in shelf stability that might result from NEMO explicitly calculating melt values for grid-boxes closer to the UniCiCles grounding line. Since the UKESM1.0-ice AIS configuration of UniCiCles is required to maintain a minimum ice thickness of 10 m and has a fixed-front calving model, the surface extent of the ice shelf on the NEMO mesh does not change through this process, and no new areas of open water can be created.

The depth of the ice shelf base on the NEMO grid computed by this script is an intermediate result. Further constraints apply to the number of active layers and partial layer thicknesses, and also to generating an acceptable ocean state in active grid-boxes for temperature, salinity and components of velocity, considering the Arakawa C-grid (Arakawa & Lamb, 1977) used in NEMO. These constraints are taken into consideration in an initialization step inside NEMO when the next segment of the climate model runs, according to the algorithms outlined in Appendix A.

Non-Conservation in NEMO From Ice Shelf Coupling

Our ice coupling scheme prioritizes model stability, but leads to potential conservation issues in the full ESM. The largest of these arises through changes in ocean volume in NEMO resulting from the ice-shelf—ocean coupling.

In reality, the conversion between freely floating solid ice and seawater at the ice shelf base makes no barystatic contribution to GMSL because the weight is supported by the ocean in either form. In contrast, in our model the shape of the shelf does not change while NEMO is running, but melt calculated at the ice ocean interface does add new water to the ocean. A volume $\mathrm{\Delta }{V}_{BMB}$ of water (positive for ice shelf loss) is added, which is rapidly redistributed by barotropic gravity waves, affecting GMSL.

UniCiCles updates the shape of the ice sheet as a result of SMB from the atmosphere model, ice shelf BMB from NEMO, and its own ice dynamics. The ice draft in NEMO is then updated accordingly, changing the seawater-equivalent volume of the ice shelf by $\mathrm{\Delta }{V}_{shelf}$ (negative for ice loss), the change in ice volume below the NEMO sea-surface calculated from ice thickness using Archimedes' principle. When we modify the NEMO state to match this new cavity under the ice shelf, we do not affect GMSL in NEMO, since the ice removed is replaced with seawater of exactly that volume (see Appendix A2). In a hypothetical case where a reduction of ice shelf draft were entirely due to BMB, this technique involves ”double counting” the volume of melt water added to the ocean, although GMSL is, correctly, changed only once. In this case water would be added first to the ocean as $\mathrm{\Delta }{V}_{BMB}$ as NEMO was running (changing GMSL because the ice shelf draft does not reduce as it does so), and then added again as $\mathrm{\Delta }{V}_{shelf}$ during the update of the ice draft (not changing GMSL, because ice shelf draft changes by exactly the opposite amount, $-\mathrm{\Delta }{V}_{shelf}$).

UniCiCles also supplies a seawater-equivalent volume $\mathrm{\Delta }{V}_{calving}$ (positive) of ice to NEMO as icebergs coming from the flow of ice across the calving front, whose position is fixed. When icebergs melt in NEMO, this water is added as a barystatic contribution to GMSL in the same way as ice shelf melting. In the real ocean, however, iceberg calving and melt has no effect on GMSL, since the ice concerned is already floating. Including contributions from the seawater-equivalent volume of ice transferred by dynamics across the grounding line ($\mathrm{\Delta }{V}_{ground}$) and the seawater-equivalent volume of mass change due to ice shelf surface mass balance ($\mathrm{\Delta }{V}_{SMB}$, positive for accumulation), in both reality and the model, the net change of the seawater-equivalent volume of the ice shelf is [Image Omitted. See PDF]

In the real ocean, the barystatic contribution to GMSL from changes in floating ice shelf volume actually comes from $\mathrm{\Delta }{V}_{ground}$ and $\mathrm{\Delta }{V}_{SMB}$, when the weight becomes supported by the ocean. Thus, the barystatic contribution to GMSL from processes that change a floating ice shelf is $\mathrm{\Delta }{V}_{ground}+\mathrm{\Delta }{V}_{SMB}$ in reality, but $\mathrm{\Delta }{V}_{BMB}+\mathrm{\Delta }{V}_{calving}$ in our model, and GMSL rise due to our model ocean—ice coupling only agrees with reality if these four terms balance. These conditions are not generally met, especially not in any individual year, meaning that on updating the ice shelf draft in NEMO UKESM1.0-ice produces a conservation error ($\mathrm{\Delta }{V}_{err}$) in the seawater-equivalent volume of ocean and ice shelf combined of [Image Omitted. See PDF]

Combining Equations 1 and 2, it can be seen that our $\mathrm{\Delta }{V}_{err}=-\mathrm{\Delta }{V}_{shelf}$, regardless of the physical processes acting on the shelf. Consider the case where both surface and basal mass balance are 0: $\mathrm{\Delta }{V}_{BMB}=\mathrm{\Delta }{V}_{SMB}=0=>\mathrm{\Delta }{V}_{shelf}=\mathrm{\Delta }{V}_{ground}-\mathrm{\Delta }{V}_{calving}$. In that case, barystatic GMSL rise would be $\mathrm{\Delta }{V}_{ground}$ in reality and $\mathrm{\Delta }{V}_{calving}$ in the model, and the change in volume would only be conservative if $\mathrm{\Delta }{V}_{shelf}=0$ i.e., the ice shelf were in balance. In a balanced state without melting, the same ice flux crosses the grounding line and the calving front and it makes no difference to GMSL where it is counted. This principle still holds when melting is present, and our coupling produces volume fluxes with the correct effect on GMSL only when the ice shelf is in steady state. Since $\mathrm{\Delta }{V}_{shelf}$ lost by the ice shelf is exactly equivalent to the volume gained by the ocean during the ice draft update, our non-conservation can be seen to be exactly the same as the total water added or removed from NEMO when we modify the ocean state to fit the new ice draft (see Appendix A2).

Another way of viewing this is to note that, by formulation, water is globally conserved in NEMO (and JULES) under an assumption that land ice is in a perfect steady state. For example, NEMO's BMB paradigm assumes that loss of volume from the ice shelf to the ocean from basal melt is immediately replaced by ice flowing from upstream. Any attempt to impose changes in ice sheet geometry without reformulating this basic paradigm must thus trigger some non-conservation, and that non-conservation is, by definition, the degree to which the ice has been moved away from its initial, assumed-steady state. Although the absolute magnitude of the error in our simulated Earth System is simply equivalent to the change in volume of the shelf, the conceptual cause of this error lies not just in how we realise that change in volume, but comes also from the unrealistic assumption in NEMO and the other climate models that the ice shelf geometry does not evolve.

Heat is also not well conserved in our implementation of the ocean-ice coupling. A conductive heat flux between the ocean and ice shelf is prescribed in NEMO, although this is not used yet as a basal boundary condition in UniCiCles, nor is it passed through to the atmosphere model above. Icebergs calved from the shelf also do not carry a temperature from UniCiCles through into the Lagrangian scheme in NEMO. Icebergs in NEMO are assigned an internal temperature of −4$°\mathrm{C}$, but the empirical equations determining iceberg melt are only weakly influenced by ocean-iceberg temperature difference. Melt added to the ocean from icebergs is assumed to be at 0$°\mathrm{C}$, and only the latent heat required for the phase transition is removed from the ocean. The non-conservation in volume and heat caused by ice shelf coupling could in principle be removed from NEMO via local or global corrections, and we are exploring methods that do this, but they are not mature enough for use at this point.

The non-conservation of ocean volume in our implementation of ice shelf coupling has little practical impact in UKESM1.0-ice, since our inability to move the land-sea mask means that any changes in GMSL are not able to be expressed in our model. The correct contribution to GMSL from the change in ice volume above flotation in our simulations can be diagnosed accurately from UniCiCles. However, non-conservation is an unwanted side-effect of the ice coupling in this version of UKESM and should be monitored during simulations, especially those on centennial timescales or with very large changes in ice shelf volume.

Example Results

As examples of the ice-shelf—ocean coupling, this section will include a configuration used for the idealized MISOMIP1 intercomparison (Asay-Davis et al., 2016), and a more realistic UKESM1.0-ice AIS in a global model. As noted in Section 5, UKESM1.0-ice does not include ocean coupling for Greenland. Even more so than for the surface coupling examples in Section 4.2, the basal ice shelf coupling examples here are intended as technical exercises, illustrating that our coupling techniques work and can be used in the context of a global ESM.

MISOMIP1

Initial development of the ice-shelf—ocean coupling for UKESM1.0-ice was done in the context of MISOMIP1, one of a coordinated suite of multi-model intercomparison exercises (Asay-Davis et al., 2016). MISOMIP1 uses a small, idealized basin geometry, with an ice shelf attached to an ice sheet grounded on a sloping bathymetry that allows unstable retreat of the grounding line (Schoof, 2007) to occur in part of the domain.

Participants were required to conduct a set of simulations with common resolutions and parameterisations (denoted COM), and allowed to contribute a less-constrained second set using their own preferred grid and science configurations (TYP). The 2km-resolution COM simulation provides a useful benchmark when compared with TYP results to demonstrate the impact of ocean resolution in coupled ice—ocean simulations. In an initial calibration step, the transfer coefficients for heat and salt in the boundary layer under the ice in each configuration are tuned to achieve a common melt rate in a simple reference configuration to ensure that the different model setups start in similar regimes. For our TYP simulations we used the same scientific configuration of NEMO as in UKESM1.0-ice at eORCA025 resolution ($\sim$8 km at the nominal latitude of the domain). For MISOMIP1, outside of the context of the full UKESM1 infrastructure, we used a different file system coupling to pass information to and from BISICLES which nevertheless used the same approach as implemented in UKESM1.0-ice.

In the MISOMIP1 “IceOcean1” experiment, strong warming is applied at the far end of the ocean domain for the first 100 years of the simulation and the ice shelf retreats. In the second half of the experiment the ocean is cooled again, allowing the ice sheet to re-advance. The COM and TYP configurations of our coupled NEMO-BISICLES model show a similar overall sensitivity to the warming, with the grounding line retreating $\sim$18 km in COM and $\sim$20 km in TYP in the first part of the experiment (Figure 10). Similar reductions in ice thickness in all parts of the shelf are seen in both models with the maximum melt rates concentrated at the grounding lines. However, these similar overall results come from different realizations of ocean physics in places: the circulation along the shelf is much stronger in the 2 km ocean model, with a sharp boundary current that forms where water exits the cavity. In contrast, currents in the eORCA025 configuration are slower and broader. The tuning of the heat transfer coefficient in the initial calibration step makes shelf melt in TYP more sensitive to the thermal gradient across the boundary layer, but the same amount of melt and retreat occurs in both simulations.

Figure 10. Simulated properties under the ice shelf in the MISOMIP1 IceOcean1 experiment after 100 years of retreat under warm ocean forcing. Top (a–c): 2 km reference simulation, bottom (d–f) eORCA025 resolution. (a and d) melt rate (m/yr, negative values indicate refreezing), (b and e) friction velocity (mm/s), (c and f) change in ice draft (m, negative values indicate thinning). Axis labels show distance in km relative to the bottom left of the combined ice sheet—ocean domain. The grounding line is located at around x = 120 km at the start of the experiment.

UKESM1.0-Ice With Interactive Antarctica

The idealized, gently sloping boundaries of the small MISOMIP1 domain do not provide a sufficient test of the robustness of an ice—ocean coupling scheme. Our first attempts at moving the ice shelf position in NEMO and updating the ocean state worked well in the MISOMIP1 experiments but proved prone to instabilities in many cases when applied to the more realistic representation of the Antarctic coast and ice shelves in UKESM1. A realistic simulation of the behavior of Antarctic ice shelves in an ESM requires much more than a well functioning coupling system. The ocean, ice and atmosphere models must all individually represent their parts of the system well, and should not feature biases that disrupt the other components unacceptably when coupled. On grounds of resolution alone, we would want to use the eORCA025 ocean configuration when simulating interactive Antarctic ice shelves in UKESM1.0-ice, since the eORCA1 mesh provides very few grid-boxes under some crucial outlet glaciers - for example, eORCA1 has fewer than 10 open columns under Pine Island Glacier. However, it is common for global climate models to feature significant biases in the Southern Ocean and Antarctic Seas (Beadling et al., 2020), and the eORCA025 configurations of UKESM1 and GC3.1 have particularly large biases. Storkey et al. (2018) note a very weak Antarctic Circumpolar Current that is linked to an unrealistic eastward reach of an extended Ross gyre, and anomalous transport of cold and fresh water masses onto the continental shelves of the Amundsen and Bellingshausen Seas. These poorly simulated shelf waters are preexisting features of the climate model simulation, not set up by the ice shelf melt or coupling, but they mean that realistic simulations of Antarctic shelf—ocean interactions cannot be carried out with the eORCA025 version of UKESM1.0-ice. The more widely used configuration of UKESM1 with an eORCA1 ocean has much more acceptable Southern Ocean water properties, so this is the version we are currently developing for ice coupling despite the poor resolution of smaller ice shelves. Even with these drawbacks, output from development simulations of UKESM1.0-ice using these ocean resolutions illustrates that our ice shelf—ocean coupling can cope with complex topographies, a wide range of ocean conditions and a significant amount of grounding line retreat.

Initializing consistent states across the components of ESMs is challenging. The millennial timescales of ice sheet adjustment, along with the fact that neither the GrIS or AIS are in equilibrium for the preindustrial period that is commonly used as a stable baseline for climate simulations, makes initializing ice sheet components particularly difficult. For our proof of concept UKESM1.0-ice simulations, drift that results from inconsistencies in the initial states is immaterial, as we are more interested in demonstrating technical capability than producing accurate climate projections. We initialize the ocean from rest with temperature and salinity derived from the EN4 dataset (Good et al., 2013). The initial atmosphere state comes from the UKESM1.0 preindustrial run (Kuhlbrodt et al., 2018). In the simulation with an eORCA1 ocean, UniCiCles is initialized directly from the modern state of a standalone BISICLES simulation (Martin et al., 2019). For the run with an eORCA025 ocean, the atmosphere and ocean were run together under constant 1970 greenhouse gas and other forcings for 5 years to generate average SMB and BMB fields that could be used with UniCiCles during a short additional standalone ISM spinup of 10 years, to overcome an initial coupling shock (measured by a spike in the continental average rate of ice thickness change). Both UKESM1.0-ice versions were then run for several decades with all components coupled under 1) constant 1970 greenhouse gas and other forcings and 2) instantaneously quadrupled 1970 CO₂ concentrations (“abrupt4xCO2”), to assess the sensitivity and computational stability of the coupled system.

The simulations with constant 1970s forcings can be used to illustrate problems when coupling Antarctic ice shelves to eORCA025 configurations of UKESM1. After 60 years the cold, fresh bias that has developed in the Amundsen region in the eORCA025 model has put the cavities around West Antarctica into an unrealistic melting regime, and there is a significant warm bias in Weddell Sea and around East Antarctica (Figure 11). By contrast the eORCA1 model has much smaller biases in temperature and salinity and higher area-average melt rates under shelves along the Amundsen coast than the eORCA025 configuration, even though the grid-boxes are much larger.

Figure 11. Drift in ocean conditions around Antarctica and simulated BMB under ice shelves in UKESM1.0-ice development simulations after 70 years. (a, c, and e) use an eORCA025 ocean, (b, d, and e) use an eORCA1 ocean. (a–d) are decadal average, 300–1000m depth averaged anomalies with respect to initial conditions (EN4, Good et al., 2013) for potential temperature (a) minimum = −3.93 °C, maximum = 2.28 °C; (b) minimum = −2.19 °C, maximum = 2.12 °C) and salinity (c) minimum = −31.44 psu, maximum = 0.805 psu; (d) minimum = −31.3 psu, maximum = 0.578 psu). (e and f) are decadal average ice shelf BMB (e) minimum = −51.2 m/yr, total = −829 Gt/yr; (f) minimum = −32.7 m/yr, total = −1,003 Gt/yr. Below 5 m/yr the scale interval increases to 20 m/yr). Black contour shows surface extent of shelf, gray contour shows southerly extent of ocean cavity in this depth range.

The ability to stably simulate ocean cavity circulation in regions where automated changes are being made to complex topography is fundamental to a model like UKESM1.0-ice. The sudden jump in atmospheric greenhouse gas concentrations in our configuration quickly leads to warming of the ocean under the Ross ice shelf and significant changes to ice thickness under the eastern part of the shelf. This includes many kilometres of grounding line retreat and the ungrounding of the ice sheet from several pinning points (Figure 12). The response of the AIS to the warming in this simulation is very much greater than the drift in the simulation with constant 1970 forcing, indicating that these changes are not simply a consequence of the initialization of the ice sheet. Although the eORCA grids do not include all the features of the cavity topography in UniCiCles, the irregularity of the grid and convergence of meridians close to the south pole of the ocean grid means that many features of the ISM geometry change that are smaller than the nominal $1°$ (eORCA1) or $1/4°$ (eORCA025) longitude/latitude grid scale are reproduced in the ocean. eORCA025 grid-boxes can have sides as small as 4 km, with 16 km for eORCA1. Accordingly, the patterns and range of shelf thinning simulated in UniCiCles are reasonably well reproduced in the corresponding NEMO discretization, although the eORCA1 version does miss more detail. These simulations suggest that the coupled ice sheet in the eORCA025 configuration may be more sensitive, as the Ross ice shelf thins more and the grounding line retreats further over the course of the simulation than when forced with the eORCA1 ocean. It should be noted, however, that the two configurations start with different ice sheet states and each coupled simulation has both a forced response and different rates of numerical drift away from the initial conditions so it is difficult to draw robust interpretations just from these simulations.

Figure 12. Change in the floating ice area and water column thickness under the Ross ice shelf after 70 years of abrupt4xCO2 climate as represented on the native grids of the ocean and ice models. Colors show the water column thickness (m) for regions that have ungrounded during the course of the experiment, black contour shows the extent of floating ice. Light blue area is open ocean, cream is grounded ice, white is ice that was already floating at the start of the experiment. (a and b) are from a simulation with the eORCA025 ocean, (c and d) are from a simulation with an eORCA1 ocean. Note the different color scales for the different model configurations. (a and c) show data on the high resolution BISICLES grid (maximum in a) is 771 m, in (c) is 528 m); (b and d) show the same changes after discretization onto the respective NEMO meshes in use (b) eORCA025 maximum = 653 m, (d) eORCA1, maximum = 239 m). A distance scale (km) is given on the regular BISICLES grid, the irregular NEMO grid is projected onto longitude/latitude axes.

In Section 6.2, we showed that the non-conservation of volume added to GMSL in NEMO from the way we handle the ice—ocean coupling in NEMO is equivalent to the change in the ice shelf draft over the simulation, adjusted for the relative densities of ice and seawater. These abrupt4xCO2 experiments simulate large and rapid changes in Antarctic ice shelves, so this continental-scale nonconservation would be expected to be large. In the simulation with the eORCA1 ocean the error is $\sim$1.5 mm/yr out of a total of $\sim$9 mm/yr from all mass and steric contributions to GMSL rise in NEMO. In the simulation with the eORCA025 ocean, shelf retreat is larger and the error amounts to $\sim$3.6 mm/yr of the 16 mm/yr total. This represents a significant degree of non-conservation in the NEMO global water budget, although it has a limited practical impact on our ESM component models, as discussed in Section 6.2.

The ocean-forced thinning of the Ross ice shelf and retreat of the grounding line has implications for the dynamical flow of the ice. The loss of buttressing from the ice shelf would be expected to speed up flow of the grounded ice whilst flow on the shelf slows down as the ice thins (e.g., Gagliardini et al., 2010). In our simulation, the glaciers that feed the shelf more than double in speed whilst flow reduces toward the shelf edge, reducing the solid ice mass calved from the shelf (Figure 13). Although the broad response is the same across the configurations with different ocean resolutions, with the eORCA1 ocean there is less acceleration of the grounded ice, consistent with the grounding line response. The acceleration on the western side of the Ross ice shelf in the simulation coupled to the eORCA1 ocean comes from small areas of high SMB forcing in that area.

Figure 13. Net change in ice velocity and thickness on the Ross ice shelf after 70 years of abrupt4xCO2 climate. Colors show change (year 70 - year 0) in velocity, black contours show magnitude of thinning, blue contour shows the retreated grounding line at the end of the simulation, gray contour shows the edge of the floating ice. (a) For the simulation with an eORCA025 ocean (extremes are a slowing of 637 m/yr, a speed up of 1,238 m/yr, 702 m of thinning, 348 m of thickening); (b) the simulation with an eORCA1 ocean (extremes are a slowing of 309 m/yr, a speed up of 4,299 m/yr, 596 m of thinning, 585 m of thickening). The Cartesian BISICLES grid is projected on longitude (°E) latitude (°N) axes.

As shown, when UKESM1.0-ice is run with high CO₂ concentrations it simulates significant ocean warming and ice shelf retreat in the Ross Sea. Analyzing the physics and realism of this change are not within the scope of this paper, but they are under further study. On the basis of these examples, we are confident that our ocean—ice coupling scheme technically functions well in the context of a realistic Antarctic geometry, subject to changes that UKESM1 will likely need to simulate for future climate projections.

Discussion and Future Work

We have created a practical system by which dynamical models of the Greenland and Antarctic ice sheets may be interactively coupled into the UKESM1 Earth System model. Mass fluxes due to ablation and accumulation at the atmosphere-ice interface (surface mass balance, SMB) and ocean-ice interface (basal mass balance, BMB) are calculated by the climate model. Through the coupling these fluxes can drive growth and retreat of the ice sheets. Land-based ice sheet margins, surface topography, grounding line and ice shelf thickness can all evolve consistently in the ice and climate models. Solid ice discharged at the marine margins of the ice shelves produce icebergs, which drift and melt in the ocean.

Our system allows us to address a new range of science questions with UKESM1. These concern not just the processes that contribute to GMSL but also the impact of changes in ice sheet topography on atmospheric circulation and SMB, the effect of ice sheet meltwater on ocean stratification and circulation, and the impact of coupled changes on ice sheet instabilities (e.g., Schoof, 2007). All of these effects could feed back onto ice sheet evolution, and hence affect ice sheet mass loss and future GMSL rise as a consequence of anthropogenic climate change. This is an emerging area for study in global models, and other ESMs with interactive ice sheets or ocean cavities under ice shelves are also starting to investigate these questions (e.g., Jeong et al., 2020; Muntjewerf et al., 2020). Detailed results and analysis from simulations that use our new capabilities in UKESM1 will be presented in forthcoming papers.

Two-way coupling of ice sheets in ESMs is a new field that must deal with a number of scientific and technical difficulties. Our system suffers from a number of limitations that we wish to address in future work, especially regarding conservation, the inclusion of marine influence on the margins of the GrIS and allowing Antarctic ice shelves to advance or retreat their ice fronts. Our current scheme does represent changes in shelf grounding lines and buttressing effects on the grounded ice sheet but does not allow us to fully simulate the impact of shelf collapse on the local Earth System. It is especially difficult to alter the predefined surface coastlines of the land-sea mask in the UKESM system and the technical barriers to removing this limitation in our model are significant, but it is an area of active development.

The spatial resolution of our global ocean model is not fine enough for it to represent the marine-terminating glaciers of GrIS, the turbulent eddies and plumes of buoyant meltwater that affect Antarctic ice shelf BMB, or the narrow canyons and steep slopes that determine the flow onto and off the Antarctic continental shelf. Parameterisations for translating distant open ocean properties into marine boundary forcing conditions for the ice sheet grid have been proposed (e.g., Jourdain et al., 2020; Slater et al., 2020) and we are considering how we might adapt these for UKESM1 in ongoing work. Without basal melting, the modern-day simulated ice discharge from GrIS in UKESM1.ice is lower than observed, and evolves differently in the coming century when compared to more detailed standalone BISICLES simulations that do include marine forcing in their parameterization of calving. Large-scale retreat of the ice sheet from the coast is likely to limit the role that marine forcing and calving plays in the future evolution of GrIS (Gregory et al., 2020), although some studies suggest that ocean influences on outlet glacier dynamics may still be influential on centennial timescales (Aschwanden et al., 2019; Choi et al., 2021).

The low spatial resolution of global atmosphere and land models causes problems for surface mass balance simulation, which we mitigate by implementing subgrid-scale elevation tiles. However, this scheme does not help with boundary variables that do not vary consistently and monotonically with altitude. Improvements are possible, such as schemes used in other models to take account of the reduction of precipitation at high altitude in Greenland, and repartitioning between liquid and solid phases (e.g., van Kampenhout et al., 2020). Resolution-dependent atmospheric phenomena, such as synoptic features in the Atlantic and katabatic winds in Antarctica, are more difficult to address. Other approaches, such as nested higher resolution domains over the ice sheets have been explored (e.g., van Kampenhout et al., 2019), which suggest that simply increasing atmospheric resolution over the ice does not automatically result in an improved simulation of ice sheet boundary conditions.

Many scientific questions concerning the co-evolution of ice sheets and climate require simulations of many centuries or millennia. The computational cost of the atmosphere and ocean models in UKESM1 make such simulations prohibitively expensive and time-consuming. The same applies where ensembles of simulations are required to deal with uncertain initial conditions or parameter settings. The asynchronous coupling scheme described in Section 3 can only be used when the atmosphere and ocean climates are changing much more slowly than the ice sheet, which is not true when radiative forcing changes rapidly or the ice sheets drive rapid changes in climate. In addition, the non-conservation inherent to our ice coupling scheme may start to materially affect the simulation of the model on long enough timescales. For such questions, something like UKESM1 is unlikely to be an ideal choice of model, and simpler or lower-resolution climate models coupled to an ISM are still necessary (e.g., Smith et al., 2021).

Despite all the above, the main difficulties we have found with global ESM simulations including interactive ice sheets is the unrealistic simulation of key physical processes in the climate components - especially in the ocean surrounding Antarctica - rather than in the ISM or the coupling itself. These problems are not confined to UKESM1 (e.g., Beadling et al., 2020). Insufficient resolution in global climate models is certainly part of this problem, but it is clear that simply increasing the resolution does not provide a simple solution. What is needed, rather, is more attention to be paid to the simulation, parameterization and evaluation of climate physics in these remote and poorly observed regions of the Earth.

Conclusions

UKESM1 has been coupled to dynamic models of Greenland and Antarctic ice sheets, making it the first CMIP6 model with an interactive, dynamic Antarctic ice sheet component. This work represents a major model development, requiring modification to science and infrastructure code of the land surface, ocean and ice sheet components of UKESM1 as well as the coupling framework itself. Our techniques, and their advantages and drawbacks are described and discussed. We demonstrate that our scientific approaches to modeling processes at the ice—land and ice—ocean interfaces produce plausible results in the context of a global ESM. We demonstrate that our techniques for moving the solid boundaries of the grounded ice sheet and floating ice shelves are robust and remain computationally stable when simulating extreme ice sheet retreat.

A AppendixGenerating a Valid NEMO State When Ice Shelf Thickness Changes

In the following, we consider vertical stacks (”columns”) of NEMO grid-boxes, bounded by the solid sea floor below and the bottom surface of the ice shelf above. Tracer (here, temperature and salinity) values are defined at the centers of these boxes, and velocities on the faces of the boxes.

A1 Discretising the Shelf Geometry

In UKESM1.0-ice, the algorithm that discretises the ice shelf cavity geometry when NEMO is initialized described in Mathiot et al. (2017) has been changed to allow the grounding line position to be more consistently determined and to improve model stability. Unlike Mathiot et al. (2017), our procedure is not indeterminately iterative, and instead has two stages.

Stage one adjusts the ice shelf draft using depth criteria to ensure that it is compatible with the bathymetry being used. Every water column under the ice shelf is checked and adjusted if necessary. The sea floor bathymetry is never altered during a run.

• a minimum ice shelf draft $D_{min}$ is enforced by thickening the draft to $D_{min}$ where necessary.

• where the water column under the ice shelf is less than a height $H_{ground}$, the ice draft is thickened so that the ice becomes grounded, closing this column.

• where the water column under the ice shelf is greater than $H_{ground}$ but less than $H_{open}$, ice draft is reduced until the bottom of the shelf is $H_{open}$ meters above the sea floor.

• where the ice shelf draft has been reduced, the first test for $D_{min}$ is repeated and anywhere that is now thinner than the minimum draft is grounded.

In UKESM1.0-ice we currently use 10 m for $D_{min}$, $H_{ground}$ and $H_{open}$. Having $H_{open}=H_{ground}$ makes the third step above inactive in our examples, but our scheme allows for the possibility of setting different thresholds for situations where one might want to ground the ice shelf because it is too close to the bathymetry and where the cavity should be opened up to allow the circulation of water to be more accurately modeled.

The second stage addresses the discretization of all ungrounded water columns under the ice shelf:

• on every side where an ungrounded water column has an ungrounded neighbor, there must be at least two open velocity faces in the column to connect them. The ice shelf draft can be reduced to ensure this condition.

• a flood-fill algorithm is used to find columns that are not connected to the main body of the ocean, and these columns are grounded. Extreme conditions may develop in these isolated areas that lead to instabilities if they later become part of the main ocean as the grounding line retreats.

• grid-boxes that are surrounded by ice and share no open velocity faces with a neighbor on any side are filled with ice. Narrow chimneys in the ice are examples of this.

The resulting discretisations of the realistic ice shelf geometries we have tested allows for stable flow between all ocean grid-boxes under the shelf, without narrow areas where the ocean circulation is not resolved that could lead to the development of instabilities. The resulting representation of the ice shelf in NEMO does not imply a conservation of ice mass with respect to the ice shelf in UniCiCles, although we are concerned with the non-conservation of water in our system (see Section 6.2) and are developing means to correct it.

A2 Modifying the NEMO State to Fill the New Domain

Once the discretization of the ice shelf geometry has been determined, the ocean state needs to be adjusted to fit the new domain. This task includes initializing tracer values and velocities in newly opened wet grid-boxes, removing values from grid-boxes that have become solid ice (”dry” grid-boxes), or modifying the state of partial-height grid-boxes whose volume has changed. This state adjustment needs to be done in a way that provides stable, continuous behavior with respect to the previous period of ocean integration.

Tracer values are adjusted following Favier et al. (2019):

• If a grid-box is wet both before and after the ice shelf has been moved, its properties remain unchanged, even if it is a partial cell whose thickness has changed.

• If a new grid-box is created in a column, then tracer values are obtained by extrapolating from neighboring cells.

• If a new water column is created, tracers values and a sea surface height are obtained by extrapolating from neighboring columns.

• If a wet cell is closed then the cell is masked and its properties are lost.

The algorithm and extrapolation described above is applied to tracers and sea surface height as in Favier et al. (2019). This procedure minimizes the creation of new density gradients that could disrupt the circulation.

Velocity values must also be adjusted. It is not sufficient to simply initialize new wet grid-boxes with zero velocity and remove ones that have closed. This can result in large instantaneous changes in the local horizontal divergence, producing high vertical velocities and unphysical sea surface height tendencies in the first model timestep. Favier et al. (2019) address this problem by adjusting velocities in each column to conserve the barotropic transport. This method cannot be applied when an entire water column is grounded and we find it often leads to unstable numerical artifacts when used with realistic AIS shelf geometries in UKESM1.0-ice. Instead, for just the first timestep, we artificially force the three-dimensional divergence field to be unchanged across the change in discretisations by adding artificial volume fluxes where necessary. This prevents the formation of instabilities from the instantaneous change in geometry. At the end of the first timestep NEMO recomputes a full new velocity field itself, consistent with both the new boundary locations and tracer values and the artificial additions are no longer required. The heat and salt flux associated with this added volume is computed using the local temperature and salinity. We have found this method to be extremely computationally stable, even when faced with unrealistically large changes in ice shelf geometry. As with the shelf discretization above, our manipulation of ocean tracers and dynamics through these processes does not in general conserve heat, salt or mass with respect to the NEMO state that existed before the coupling with the ice sheet.

Acknowledgments

The authors would like acknowledge Erica Neiniger for helping construct the UKESM1 ice task cycle in Cylc, Christopher Bull for collaborating on the MISOMIP1 simulations, and the UKESM core development team on whose work this model builds. The authors also thank Xavier Fettweis for providing the MAR simulation data used here for comparison. Model development and simulations were carried out on the Monsoon and NEXCS collaborative High Performance Computing facilities funded by the Met Office and the Natural Environment Research Council. The authors also acknowledge the UK JWCRP Joint Marine Modeling Programme for providing support and access to ocean model configurations. R. S. Smith, A. Siahaan, V. Lee, A. J. Payne, P. R. Holland and C. G. Jones were funded by the National Environmental Research Council (NERC) national capability grants for the UK Earth System Modeling project, NE/N017978/1 and NE/N01801X/1. C. G. Jones also acknowledges funding from the European Commission under the H2020 Research grant no. 641816 (CRESCENDO). J. Ridley and P. Mathiot were supported by the Met Office Hadley Center Climate Programme funded by BEIS and Defra. The authors would also like to acknowledge the efforts and constructive comments of our reviewers and editor, who helped to improve this paper.

Data Availability Statement

UKESM1 source code and model configurations are available via the Met Office Science Repository Service (https://code.metoffice.gov.uk/trac/UKESM), UniCiCles coupling scripts and wrappers via the NCAS-CMS PUMA service (https://puma.nerc.ac.uk/trac/UniCiCles) and BISICLES source from the Lawrence Berkeley National Laboratory ANAG repository (https://anag-repo.lbl.gov/svn/BISICLES). Ocean model source code is additionally available from the Institut Pierre Simon Laplace Forge respository (http://forge.ipsl.jussieu.fr/nemo/browser/NEMO). A document listing the configurations and individual code branches used in this work is available at https://puma.nerc.ac.uk/svn/UniCiCles_svn/UniCiCles/branches/rss_ukesm/ukesm1.0-ice_JAMES21_details.txt. The simulation data (apart from MISOMIP1 results) used in this study are archived at the Met Office and are available for research purposes through the JASMIN platform (www.jasmin.ac.uk) maintained by the Centre for Environmental Data Analysis (CEDA); for details please contact [email protected] referencing this paper. MISOMIP1 simulation results are archived at the Open Science Foundation (https://osf.io/3cmuw).