Simulating interactive ice sheets in the multi-resolution AWI-ESM 1.2: A case study using SCOPE 1.0

Interactions between the climate and the cryosphere have the potential to induce strong non-linear transitions in the Earth’s climate. These interactions influence both the atmospheric circulation, by changing the ice sheet’s geometry, as well as the oceanic circulation, by modification of the water mass properties. Furthermore, the waxing and waning of large continental ice sheets influences the global albedo, altering the energy balance of the Earth System and inducing climate-ice sheet feedbacks on a global scale as evident in Pleistocene glacial-interglacial cycles. To date, few fully comprehensive models 5 exist, that do not only contain a coupled atmosphere/land/ocean component, but also consider interactive cryosphere physics. Yet, on glacial-interglacial and tectonic time scales, as well as in the Anthropocene, ice sheets are not in equilibrium with the climate, and prescribed fixed ice sheet representations in the model can principally be only an approximation to reality. Only climate models, that contain interactive ice sheets, can produce simulations of the Earth’s climate which include all feedbacks and processes related to atmosphere-land-ocean-ice interactions. Previous fully coupled models were limited either by low 10 spatial resolution or an incomplete representation of ice sheet processes, such as iceberg calving, surface ablation processes, and ocean/ice-shelf interactions. Here, we present the newly developed AWI-Earth System Model (AWI-ESM), which tackles some of these problems. Our modelling toolbox is based on the AWI-climate model, including atmosphere and vegetation components suitable for paleoclimate studies, a multi-resolution global ocean component which can be refined to simulate regions of interest at high resolution, and an ice sheet component suitable for simulating both ice sheet and ice shelf dynamics 15 and thermodynamics. We describe the currently implemented coupling between these components, present first results for the Mid-Holocene and Last Interglacial, and introduce further ideas for scientific applications for both future and past climate states with a focus on the Northern Hemisphere. Finally, we provide an outlook on the potential of such fully coupled Earth System models in improving representation of climate-ice sheet feedbacks in future paleoclimate studies with this model. 1 https://doi.org/10.5194/gmd-2020-159 Preprint. Discussion started: 24 July 2020 c © Author(s) 2020. CC BY 4.0 License.

ice sheet surface, is implemented by the Positive Degree Day (PDD) method within PISM, which is a temperature index scheme based on (Calov and Greve, 2005;Reeh, 1991). In this scheme, the total number of days multiply by the excess above the freezing temperature (273.15 K) are integrated. Degree-day-factors (DDFs), appropriate for snow (5 m K −1 day) and ice (8 m K −1 day), are applied in order to estimate the amount of surface ablation. As the climate model only saves monthly averaged climate characteristics, forcing for the ice sheet model adds white-noise in order to generate synthetic daily 130 temperatures as an ice sheet forcing.
In order to facilitate simulation of sea-proximal ice sheet regions and to realize the coupling between ocean and floating ice shelves, the setup additionally includes three calving parameterizations. First, thickness calving occurs at marginal ice shelf grid points with a thickness of less than 200 m. Second, continental-shelf calving is applied to ice that crosses the 2000 mdepth bedrock contour line (contemporary bedrock elevation of ETOPO1, Amante and Eakins, 2009). Finally, we employ 135 the Eigencalving parameterization, which uses stress field divergence with a proportionality constant of 1 · 10 17 , contributes to the calving. The grounding line position (that is the transition between the grounded ice sheet and floating ice shelves) is handled by a sub-grid scheme described by Feldmann and Levermann (2015). For the basal conditions of floating ice shelves, various parameterizations are available, which are described below in the section about the coupling between the ocean and ice sheet/ice shelf. Specifically, AWI-ESM-1-2-LR uses FESOM to drive the PICO box model developed by Reese et al. (2018), 140 which is built into PISM and paramaterizes ocean/ice shelf interaction as well as sub-shelf circulation.
The PISM setup presented in this study simulates the GrIS, and runs on a polar stereographic grid with a spatial resolution of 5 km (geographic coordinate system WGS 84, EPSG code: 3413 (https://epsg.io/3413)), following the newest ISMIP6 standard. As a topography boundary condition for the simulation of ice sheets we use the BedMachine v3 bedrock topography data set (Morlighem et al., 2017). While this study focuses on the interaction between climate and the GrIS, any other past and 145 present ice sheet domains are ignored. Further ice sheet domains for Antarctica and the entire Northern Hemisphere have also been implemented.

Scope: The Standalone Coupler
The coupling interface -SCOPE, the Standalone COuPlEr -facilitates coupling between the AWI-ESM-1-1 climate model component and the PISM component. It is a standalone Python program, operating sequentially between individual short term 150 integrations ("chunks") of each of the two model components. The coupler is modular. Interfaces to both ice sheet and climate components are implemented in such a way that arbitrary models can be connected to each other, irrespective of whether or not they provide coupling-relevant information in a format that fulfills requirements by the respective other model component.
Thereby, SCOPE introduces a generic layer which can be reused when introducing a model, that has been previously coupled via SCOPE, to a different climate/ice sheet model setup. A schematic of this setup is shown in 2. 155 In contrast to other coupling systems, such as OASIS (Valcke, 2013) or MeSSY (Garny et al., 2019), no actual modification to each model's source code is required to connect two model systems. Coupling is divided into two phases. In the first phase the first model (henceforth the sending model) provides relevant coupling information to the generic layer. This does not only include actual output fields, but also contains a description of the metadata that is relevant for appropriate interpretation and processing by the receiving model, such as variable names, units, and the simulation grid or mesh. The model information is 160 anonymized such that any information hinting at the field origins are removed, e.g. ECHAM6 information is simply referred to as "atmosphere" forcing, or PISM information is simply referred to as "ice" forcing.
After this step has been completed, in the second phase SCOPE is responsible for appropriately processing the data for suitable consumption by the receiving model. Per design SCOPE does not necessitate any modification to a particular model's source code as long as a suitable interface for external reading of forcing data exists. Consequently, the receiving model can 165 read the newly provided forcing data either during restart or during the run itself. As typical general circulation models include such an interface, connecting SCOPE to a new model component does not impose great technical challenges. A standard workflow for receiving data is therefore divided into 3 parts: 1. regridding the data to fit onto the receiving model's grid 2. performing any variable manipulations or corrections, including unit corrections or name changes 170 3. setting up the model-internal interfaces to actually read the transformed fields; thereby enabling communication between the receiving model and the newly generated forcing data SCOPE is configured via YAML files, allowing for a separation between user requirements and actual code logic. An example configuration file (Listing 1) for coupling between ECHAM6 (sending) and PISM (receiving) is shown below. Full examples of the coupling within the entire AWI-ESM-1-1 and PISM coupled model system are available with the source code, 175 please refer to the statement on code availability.

Computational Expense
Naturally, model throughput strongly depends on the machine on which the simulation runs and on the employed node configuration. Within the limits of parallel scalability, to some extent one can enhance throughput by employing more computing nodes, which leads to a trade-off in increased overall computational expense. For our model setup, which have balanced to 180 provide a compromise between throughput and computational expense, we derive a typical integration wall time for one model year (computed on compute nodes (2x 12-core Intel Xeon E5-2680 v3 (Haswell) @ 2.5 GHz) OR compute2 nodes (2x 18core Intel Xeon E5-2695 v4 (Broadwell) @ 2.1 GHz) "mistral" at the DKRZ in Hamburg, Germany) of ≈ 70 minutes (88% of integration time) for the climate system, ≈ 5 minutes (6% of integration time) for the coupling procedures, and ≈ 5 minutes (6% of integration time) for the ice sheet model, leading to a throughput of approximately 20 simulated years per day, ne-185 6 https://doi.org/10.5194/gmd-2020-159 Preprint. Discussion started: 24 July 2020 c Author(s) 2020. CC BY 4.0 License. such a procedure may or may not be warranted.

Climate/Cryosphere Interactions
In the following section, we describe the interactions between a climate model and an ice sheet model, using our new AWI-ESM-1-2 as an example. Additionally, we note some specific limitations and improvements that still need to be addressed.

195
In direction from ice sheet to atmosphere, the coupling procedure is as follows. After a PISM simulation chunk, coupling information is sent back to the atmosphere model. The data conveyed by PISM back to the atmosphere provides changes in surface elevation, surface ablation and runoff, and changes in the ice sheet extent. These fields are suitable to update ice sheet masks and orography in an atmospheric simulation, as well as representing ice sheet related freshwater fluxes in the climate simulation. The latter are routed within the atmosphere to the ocean via AWI-ESM-1-1's hydrology scheme.

200
When receiving this information, the ECHAM6 atmospheric model must perform several steps. First, relative changes to the orography and corresponding gravity wave drag parameters are applied corresponding to the height change sent by the ice sheet. Such changes cannot be directly implemented into the atmosphere model's restart files. Hence, new boundary conditions are computed, based on both the climate model state before the ice sheet simulation and the anomaly provided by PISM. The new boundary condition is prescribed as a target, which the atmosphere model linearly approaches over the next simulated 205 year of the climate simulation chunk. Next, changes in the ice sheet extent are used to update the glacial mask, which is shared between ECHAM6 and JSBACH. A thickness cutoff of 5 m of ice is used to determine areas which should be glaciated. As the JSBACH simulation also calculates soil layer moisture, the amount of soil moisture beneath a newly advanced ice sheet is stored and shielded of from the hydrological cycle; should the glacier retreat again, this water can be restored to maintain longterm conservation of mass. Vegetation beneath a newly advanced ice sheet is removed. In the case of a retreating ice sheet, 210 the vegetation module initializes deglaciated grid cells with bare soil, and the potential to grow tundra. This facilitates the emergence of plant life should climatic conditions favor vegetation development. Finally, ice mass loss from calving, grounded basal mass loss and surface mass loss is incorporated into the hydrology scheme and transported to the ocean. We note that the transformation of calved ice directly into liquid water is a model simplification, and in reality the freshwater flux from icebergs would instead be released after over time as the iceberg melts. In the case of ice sheet mass gain, a corresponding 215 amount of water is removed from the hydrology scheme in order to maintain mass conservation in the coupled climate/ice sheet system. Here, we can divide the total amount of precipitation over the ice sheet into two parts. Firstly, that part which can be used for ice sheet mass gain, and secondly, that part that remains liquid. In a model system without an included dynamic ice sheet component (in our case AWI-ESM-1-1 without PISM), the total precipitation over ice sheets is directly incorporated into the runoff scheme and transported back to the ocean, thus assuming any glaciated area in JSBACH is in steady state with 220 the climate forcing. However, this neglects the retention potential, i.e. a "buffer" on long time scales, that the ice sheets impose on the hydrological cycle by storing vast amounts of freshwater in the form of frozen ice. As the ice sheet grows, precipitation must be removed from the exchange between the coupled atmosphere-land/ocean system, and as the ice sheet shrinks this discharge must be increased accordingly. To correctly represent this effect, the amount of hydrological discharge is corrected in JSBACH's hydrology scheme during each coupling interval.

Ocean/Ice Sheet Coupling
In this section, we first describe a SCOPE internal solution to providing mass and energy fluxes to an ice sheet model; and next show how SCOPE can set up ocean forcing for use with PISM internal solutions for ocean forcing. Finally, we specify details that are relevant to the FESOM implementation.
Once received by SCOPE, the oceanic forcing has already been anonymized by the sending process. This generic ocean forc- The first equation (Eq. 1) computes the melting point temperature (T B ) at the ice shelf base/ocean interface via a linearized equation, where p B is the water pressure at the ice shelf base, and S B is the salinity of ocean water at the base of the ice shelf.
The parameters a and b represent the salinity and pressure coefficients of freezing equation while c is a small constant offset. the three-equation model provides the basal melting rates and the apparent temperature at the ocean-ice shelf interface. The latter is essentially the pressure and salinity-dependent freezing temperature (Eq. 1) as the interface is in equilibrium with the surroundings.
Instead of using this to directly provide heat and mass fluxes, SCOPE can also produce fields for the PISM internal parameterizations. If the ocean model is able to provide sub-shelf melting rates and temperatures directly, these can be read into PISM to force melting or refreezing and to drive the basal ice shelf temperature after SCOPE manipulates the variable metadata to

Model Tests
We test our new model by simulating warm climate states described in PMIP4. We choose model setups that correspond to PI-CTRL, MH, and LIG, as described in the introduction. We examine not only the mean climate states, as would be possible with a less comprehensive atmosphere/biosphere/ocean model, but also examine changes to the ice sheets and illustrate their 275 dynamics, with a particular focus on those quantities that may influence the remainder of the climate, and potentially tip the ice sheet (Notz, 2009) and, subsequently, the coupled cyrosphere/climate system into a new state: ice sheet orography, total ice sheet volume in sea-level equivalent, and ice sheet mass changes resulting from basal and surface melting, as well as iceberg discharge.

Spinup Procedure
In order to produce typical snapshots experiments as normally presented in PMIP, we spinup our coupled climate/cryosphere simulations using a multi-step approach: 1. We spin up the climate model using reconstructed ice sheet geometries and appropriately adjusted orbital parameters and greenhouse gas concentrations. In the simulations presented here, the climate-only spinups are each integrated for more than 3, 000 years.  3. At this point, we consider both climate and ice sheet systems to be in equilibrium with the prescribed orbital and greenhouse gas values that refer to a specific PMIP4 time slice. Ice volumes are stable, and upper ocean salinity and temperature trends are negligible. We now couple both ice sheet and climate systems together in an asynchronous fashion, first running small chunks of the climate model corresponding to 3 simulated model years, followed by 25 model years for the ice sheet model, exchanging coupling information after each chunk. This is to ensure that the climate model can 295 react to changes that the ice sheet has undergone relative to the reconstructed geometries used in step 1, and that the ice sheet in turn can react to any changes in the simulated climate forcing. We do this for at least 1000 ice years, and judge the system to be in equilibrium based upon fluctuations in the mass balance of the ice sheets; which should be close 0 in an equilibrated system. In this case, any changes that may take place in the deep ocean are still in disequilibrium. The coupling interval is asynchronous to allow the ice sheet geometry to quickly adapt to the new climate forcing. 4. Once the models have readjusted to the new geometries, temperatures, and melting rates, we begin to synchronously couple the models, running each system for 3 years before exchanging information.
We show changes in total solar incoming radiation at the top of the atmosphere relative to the Pre-Industrial control experiment for the two paleoclimate snapshots, and corresponding changes in sea-level equivalent ice sheet volume during the ice-sheet only spinups from step 2 in Figure 3. One can note that the changes in ice sheet volume are in line with the ex- ing the coupled AWI-ESM-1-2 generates a total GrIS ice volume that is the equivalent of 6.7 m of sea level. These simulation results confirm that the standalone PISM setup, forced with AWI-ESM-1-1-LR climate forcing, is able to reproduce general ice sheet characteristics that are in line with the assumptions made based upon sea level reconstructions (Dutton et al., 2015).
In particular, sea level reconstructions suggest a LIG GrIS being smaller than a MH GrIS, which in turn is smaller than a PI ice sheet.

315
Next, we examine the simulated ice sheet thicknesses following the ice-sheet-only simulations as these will then be used as new orographic boundary conditions for the subsequent climate simulation (4). We additionally demonstrate the dynamics simulated in our cryosphere simulations by showing the vertically integrated ice velocity. The Pre-Industrial control simulations show similar ice sheet thickness and velocities as are presently observed. The increase in the central ice dome's thickness seen in both the Mid Holocene and Last Interglacial timeslices can be linked to changes in the SMB, which in our case is primarily 320 driven by precipitation changes. However, the increased melting at the ice margins leads to a net decrease in overall ice volume.
Intensified melting occurs in our simulations due to the warmer temperatures simulated during the summer months for both the MH and LIG. Changes in the ice velocity may be related to increased orographic gradients, leading to faster ice flow. The ice divides also shift slightly due to changes in the location of the central ice dome.

325
In this section, we briefly present several examples of forcing sent through the SCOPE, and demonstrate the influence of several of the downscaling and surface mass balance choices provided within the atmosphere/PISM interface. We also discuss the influence of these choices on the ice sheet's surface mass balance.
Under a control simulation representative of the Pre-Industrial, we test our model's ability to simulate the surface mass balance of the GrIS. Shown in Figure 5 are the inputs for the PDD mass balance scheme 2 meter (a) surface air temperature Next, we examine the inputs and outputs to the ocean interaction component between FESOM and PISM, the PICO box 335 model. Shown in Figure 6 are temperature and salinity for the control state, as well as the corresponding melt rates as computed by PICO. Note that the values have already been interpolated to the PISM grid, and that FESOM offers a spatial resolution of between 20 − 30 km around the Greenland coast; in this case the mismatch between the grids of ocean and ice sheet is not as severe as between the atmosphere and ice sheet models. Nevertheless, since the coastline representation between the two models may differ, the ocean forcing data is extrapolated via a nearest-neighbor algorithm after interpolating it to the PISM

Applications: Potential for Paleoclimate Research
Given the significant interactions between the cryosphere and the remainder of the climate system, the development of AWI-ESM-1-2 represents a major step towards producing more comprehensive simulations of Earth's past and future climate. In this section, we present preliminary results for the Pre-Industrial and compare these against the PMIP timeslices for MH and LIG; both of which are, time periods in Earth's history where paleoclimate proxies indicate warmer-than-present temperatures (Otto-Bliesner et al., 2020;Brierley et al., 2020b). As opposed to the coupling tests shown above, these simulations are fully, 350 bi-directionally coupled. First, we show the accelerated simulations during which PISM and AWI-CM run asyncronously, followed by results for the synchronously coupled experiments.

Asynchronous Experiments
The asynchronous runs serve to allow the ice sheet to adapt to the climate forcing, which after re-coupling adjusts slightly, relative to the stand-alone case, as the standalone models do not interact. During the long ice sheet-only spinups from the 355 viewpoint of the ice sheet the climate state is fixed. Consequently, the ice sheet model does not experience any of the transient feedback-driven interannual and interdecadal variability that is an inherent characteristic of the climate system. As such, immediately after coupling, the mass balance of the ice sheet, and, as a result, the ice sheet geometry, adjusts to the changed forcing/boundary conditions. Furthermore, the climate model still needs to adjust to this new ice sheet geometry. In order to speed up model equilibration, we allow this adjustment to take place in an accelerated manner. Time series of ice sheet mass 360 balance and ice volume are shown in Figure 7 for each of the paleoclimate experiments during the accelerated coupling runs.
Once the ice sheet geometry has re-stabilized, we synchronously couple the models together for a further 100 years, which are used for evaluation of the quasi-equilibrated coupled climate/ice sheet system.

Synchronous Experiments
During the synchronous experiments, each 3 years of the climate simulation are used to force 3 years of the ice sheet model.

365
After that coupling information is exchanged. We find that the globally averaged near surface air temperature change only slightly relative to a control case without interactive ice sheets for each of the paleoclimate experiments (Figure 8). For the paleoclimate experiments, seasonal anomalies for near-surface air temperature, precipitation, sea surface temperature, and sea surface salinity are shown in Figure 9 for the MH experiment and in Figure 10 for LIG.
The globally averaged timeseries demonstrates that our new model AWI-ESM-1-2-LR, that includes interactive ice sheets, 370 performs similarly to the simulation with prescribed static ice sheets for the simulations representing the PI-CTRL case. Hence, the integration of an ice sheet does not degrade the model performance. However, the strongest temperature changes relative to the PI-CTRL simulation are seen in the DECK simulations for 1% CO 2 increase as well as in the abrupt 4CO 2 simulation.
Globally averaged 2m temperatures increase from 13 • C to approximately 21 • C and 19 • C by the middle of the century. The simulations for the past warm climates show slight cooling during the Mid Holocene and Last Interglacial (from a globally 375 averaged mean of 13 • C to 12.5 • C and 12.8 • C, respectively), which is in line with slightly lower concentrations of greenhouse 13 https://doi.org/10.5194/gmd-2020-159 Preprint. Discussion started: 24 July 2020 c Author(s) 2020. CC BY 4.0 License.
gases. The strongest climatological change occurs for both of these two periods seasonally, which is expected due to the changed orbital configuration. For the Mid Holocene, winters are on average colder; with local cooling of −2.10 • C over the Northern Hemisphere high latitudes (field average north of 60 • N). Summer values are comparable to the PI-CTRL simulation.
For the LIG simulation, summers are warmer in the Northern Hemisphere, with temperature anomalies up to +7 • C. Localized

395
In the present study, we have presented the new comprehensive Earth System Model AWI-ESM-1-2-LR, comprising of ECHAM6/JSBACH 6.3.01p4/3.11, FESOM 1.4, PISM 1.1.4, the OASIS3-mct coupler, and SCOPE, our newly developed, offline, script based coupling system. We showed the various coupling strategies implemented in our system. In particular, our model allows for a choice of several mass balance schemes; a Positive-Degree-Day approach can be used, or a more sophisticated model which implicitly resolves the diurnal cycle via a parameterization, the dEBM (Krebs-Kanzow et al., 2018b).
400 Additionally, our model is able to represent climate/ice sheet dynamics resulting from ice-shelf/ocean interaction. We tested our model with various DECK simulations described in the upcoming IPCC AR6 (CMIP6 DECK simulations), as well as two paleoclimate time slices used in PMIP4. In both sets of experiments, our results agree well with the respective ensemble mean, highlighting that ice sheet interaction does not deteriorate the mean climate state. Hence we expect, that common time slice climate state simulations deliver similar results. However, since the interaction between ice sheets and the climate system 405 in large can either dampen or amplify certain climate trends, results of transient simulations will probably differ in terms of timing. Our system is now able to resolve several crucial feedbacks which are believed to be of key importance when exploring either transient or equilibrium climate states among all climate components, as it has been shown that both freshwater influence (Barker et al., 2015) as well as orographic changes play a decisive role in the North Atlantic climate in other modelling studies.
as a result, our system can only be used to prognostically understand possible sea level changes. However, future versions of the model will resolve this limitation, thus allowing the exploration of a range of scientific questions on glacial/interglacial to tectonic timescales as well as for the near-term to long-term future. The inclusion of interactive ice sheets into the model Earth System represents a fundamental improvement of the process-representation in simulations, and thus is key in advancing our understanding of climate/ice sheet interactions as well as providing a first step towards the next generation of fully interactive 415 Earth System Models. Beyond the improvements presented here, several key features would be beneficial to be implemented in the future. Of particular interest for resolving climate-carbon cycle feedbacks during glacial-interglacial cycles is a fully interactive carbon cycle including biogeochemistry. Furthermore, incorporation of iceberg calving and feedbacks with the remainder of the ocean circulation would improve the representation of the evolution of the Earth system, both in the past as well as in the future.      Globally averaged 2m temperature ( C) PI (aso) 4xCO 2 (aso) 1% CO 2 (aso) PI (asoi) MH (asoi) LIG (asoi) 4xCO 2 (asoi) 1% CO 2 (asoi)