GREB-ISM v1.0: A coupled ice sheet model for the Globally Resolved Energy Balance model for global simulations on timescales of 100 kyr

We introduce a newly developed global ice sheet model coupled to the Globally Resolved Energy Balance (GREB) climate model for the simulation of global ice sheet evolution on timescales of 100 kyr or longer (GREB-ISM v1.0). Ice sheets and ice shelves are simulated on a global grid, fully interacting with the climate simulation of surface temperature, precipitation, albedo, land–sea mask, topography and sea level. Thus, it is a fully coupled atmosphere, ocean, land and ice sheet model. We test the model in ice sheet stand-alone and fully coupled simulations. The ice sheet model dynamics behave similarly to other hybrid SIA (shallow ice approximation) and SSA (shallow shelf approximation) models, but the West Antarctic Ice Sheet accumulates too much ice using present-day boundary conditions. The coupled model simulations produce global equilibrium ice sheet volumes and calving rates like those observed for present-day boundary conditions. We designed a series of idealized experiments driven by oscillating solar radiation forcing on periods of 20, 50 and 100 kyr in the Northern Hemisphere. These simulations show clear interactions between the climate system and ice sheets, resulting in slow buildup and fast decay of ice-covered areas and global ice volume. The results also show that Northern Hemisphere ice sheets respond more strongly to timescales longer than 100 kyr. The coupling to the atmosphere and sea level leads to climate interactions between the Northern and Southern Hemispheres. The model can run global simulations of 100 kyr d⁻¹ on a desktop computer, allowing the simulation of the whole Quaternary period (2.6 Myr) within 1 month.

1 Introduction

Understanding ice-age cycles in the Quaternary period requires an interdisciplinary research approach including the fields of astronomy, geology, physical geography, oceanography and atmospheric science. Geological proxy data show that sea level and surface temperature significantly oscillated with a preferred timescale of about 100 kyr during the last million years, indicating that large ice sheets and glaciers formed and retreated many times over this period (Imbrie et al., 1984; Shackelton, 2000; Short et al., 1991). These oscillations in the late Quaternary are known as the ice-age cycles.

By investigating ice-age cycles, researchers have identified many climate processes that generate long-term climate variability. Variations in Earth's orbit and resulting changes in solar forcing have been widely accepted as a major driver of ice-age cycles (Imbrie et al., 1984; Milankovitch, 1941; Short et al., 1991; Tabor et al., 2015; Wunsch, 2004). The Earth's axis and variations in orbital parameters, such as precession, eccentricity and obliquity, can effectively regulate the incoming solar radiation on the Earth's surface and season length for both hemispheres, leading to global temperature oscillations on timescales of 20 to 100 kyr (Huybers, 2011; Short et al., 1991). Additionally, greenhouse gases, especially CO₂, are considered as a critical forcing during the late Quaternary (Shackelton, 2000). Before the industrial revolution, atmospheric CO₂ varied as an internal climate feedback originating from the ocean, biosphere or lithosphere (Bauska et al., 2018; Hogg, 2008). This carbon cycle of the Earth system significantly changes the surface energy budget and affects climate variability.

The formation of large ice sheets is an important element of climate variability over the last million years (Bintanja and Van De Wal, 2008; Ganopolski and Brovkin, 2017). During ice ages, Northern Hemisphere ice sheets can cover a significant portion of the North American and European continents (Manabe and Broccoli, 1985; Mix et al., 2001), modifying climate through changes in the albedo of snow, low surface temperature, surface elevation and sea level change (Bintanja and Van De Wal, 2008; Felzer et al., 1996; Hock, 2005; Manabe and Broccoli, 1985; Mix et al., 2001; Overpeck et al., 2006). In addition, there are many other factors that potentially affected the ice-age climate system such as deep ocean temperatures, ocean and atmospheric circulation changes, vegetation cover, and atmospheric dust content. The interactions between these climate elements led to a complex picture during the Quaternary, and the details of these interactions still remain unclear.

Numerical modeling of the ice–climate coupled system is an important way to investigate the effect of ice sheets on the Quaternary climate system. In the early stage, climate models only simulated the atmosphere and ocean, and ice sheet variations were included as external forcing (Bush, 2004; Gates, 1976; Manabe and Broccoli, 1985; Webb et al., 1998). Most studies with numerical simulations focused on a specific period, like the last glacial maximum, and specific regions, like the Northern Hemisphere (Bush, 2004; Webb et al., 1998), due to limitations in computational resources. Ice sheet modeling at continental scale in response to orbital forcing requires the simulation of long (> 10 kyr) periods, due to the relatively slow ice sheet adjustment time to climate forcing. Numerical studies at large spatial and temporal scales therefore often use decoupled simulations with surface temperature and precipitation taken as boundary conditions for ice sheet models (Greve, 1997; Huybrechts, 2002; Payne et al., 2000). Fortunately, thanks to computer and model developments, progressively more studies apply coupled ice–climate simulations on timescale of 100 kyr to 1 Myr (Abe-Ouchi et al., 2013; Ganopolski et al., 2010; Tigchelaar et al., 2019; Willeit et al., 2019). However, as far as we know, there are currently no global, million-year, coupled ice-sheet–climate simulations available.

In this study, we introduce a fully coupled ice-sheet–climate model as a tool for paleo-climate research. The model is capable of simulating global, coupled ice–climate simulations of 100 kyr within 24 h on a desktop computer. It is designed for studies of global interactions between ice sheets and climate on timescales of 100 kyr to 1 Myr. The starting point for this development is the Globally Resolved Energy Balance (GREB v1.0) climate model, which simulates the fast climate feedbacks relevant for the climate response to external forcing, such as CO₂ concentration or variations in solar radiation, on timescales of up to 500 years (Dommenget and Flöter, 2011; Stassen et al., 2019). We introduce a new ice sheet model (ISM) into the GREB model, defining the new GREB-ISM model.

The study is organized as follows. First, we introduce the datasets used, followed by the core of the paper which describes the GREB-ISM. This section is organized in three parts: a short introduction of the original GREB model, followed by descriptions of the new ice sheet model and the changes made to the climate simulations in the GREB model to couple the climate system to the ice sheet model. In Sects. 4 and 5 we present a series of stand-alone ice sheet and fully coupled ice–climate simulations to evaluate the performance of the new model. The final section provides a short summary and discussion.

2 Data

Input values for most climatology for the GREB model, such as surface temperature, atmospheric humidity, horizontal winds and vertical air motion, are taken from the ERA-Interim dataset (Dee et al., 2011). Soil moisture is from NCEP reanalysis data from 1950–2008 (Kalnay et al., 1996), cloud cover climatology from the ISCCP project (Rossow and Schiffer, 1991) and ocean mixed layer depth climatology from Lorbacher et al. (2006). Precipitation data is from the Global Precipitation Climatology Project (GPCP, Adler et al. 2003), and for Antarctica we use the dataset from NCEP-DOE (Behrangi et al., 2020; Kanamitsu et al., 2002).

The modern observed bed topography and ice thickness data for Greenland and Antarctica are obtained from BedMachine (Morlighem et al., 2017, 2020), Greve (1997) and Martin et al. (2011). Ice surface velocity data come from the Making Earth System Data Records for Use in Research Environments (MEaSUREs) program (Joughin, 2017; Joughin et al., 2010; Mouginot et al., 2012; Rignot et al., 2011, 2017). In this study, the bed topography refers to all different types of ice basis. Figure 1 shows the global map in the GREB model resolution of the bed topography and observed ice thickness. Ice sheet calving rates are taken from Bigg (1999) for Greenland and Liu et al. (2015) for Antarctica. For paleoclimate proxies, the Greenland Ice Core Project (GRIP) (Greve, 1997; Johnsen et al., 1997) data are used to impose surface air temperature anomalies for the last 250 kyr. δ¹⁸O proxy from sea sediment (Imbrie et al., 1984; Imbrie and McIntyre, 2006) is used as a proxy for global sea level change for the last 250 kyr. The surface temperature and precipitation during the Last Glacial Maximum (LGM) for a forced transition experiment is obtained from CMIP6.PMIP.AWI.AWI-ESM-1-1-LR datasets (Shi et al., 2020) from the Paleoclimate Modelling Intercomparison Project (PMIP4, Kageyama et al., 2017).

Figure 1 Initial ice thickness (a, b) and bed rock (c–e) in GREB-ISM. Ice thickness less than 10 m is not shown.

Table 1 Symbol and parameters list for the GREB-ISM model.

3 Model description

Before we introduce the ice sheet model developed in this study, we give a short description of the GREB model. We discuss changes made to the GREB model to couple the climate variables of the GREB model to the ice sheet variables, introducing the new model: GREB-ISM. All variables of the GREB-ISM model, as discussed in this study, are listed in Table 1. A model schematic of the coupling between the ice sheet and the climate model is illustrated in Fig. 2.

Figure 2 Schematic illustrating the coupled GREB-ISM.

3.1 The Globally Resolved Energy Balance (GREB) model

The GREB model is developed and fully described in Dommenget and Flöter (2011), with the additional introduction of a new hydrological cycle model in Stassen et al. (2019). The model has three layers (atmosphere, surface and sub-surface ocean) with a global, horizontal grid spacing of 3.75^∘ × 3.75^∘ (96 point × 48 point). The GREB model simulates four prognostic variables: surface (T_surf), atmospheric (T_atmos) and subsurface ocean temperature (T_ocean), and surface humidity (q_air):

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\mathit{\gamma }}_{\text{surf}}{\displaystyle \frac{\mathrm{d}{T}_{\text{surf}}}{\mathrm{d}t}}=& F_{\text{solar}}+F_{\text{thermal}}+F_{\text{latent}}+F_{\text{sense}}\\ & +F_{\text{ocean}}+F_{\text{correct}},\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\mathit{\gamma }}_{\text{atmos}}{\displaystyle \frac{\mathrm{d}{T}_{\text{atmos}}}{\mathrm{d}t}}=& -F_{\text{sense}}+F_{{\mathrm{a}}_{\text{thermal}}}+Q_{\text{latent}}\\ & +{\mathit{\gamma }}_{\text{atmos}}\left({\mathit{\kappa }}_{\mathrm{a}}\cdot {\mathrm{\nabla }}^{\mathrm{2}}{T}_{\text{atmos}}-\mathit{u}\cdot \mathrm{\nabla }{T}_{\text{atmos}}\right),\end{array}\text{(3)}& {\displaystyle \frac{\mathrm{d}{T}_{\text{ocean}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\mathrm{\Delta }{T}_{{\mathrm{o}}_{\text{entrain}}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_{\text{ocean}}-{\mathit{\gamma }}_{\text{surf}}}}{F}_{{\mathrm{o}}_{\text{sense}}}+F_{{\mathrm{o}}_{\text{correct}}},\text{(4)}& {\displaystyle \frac{\mathrm{d}{q}_{\text{air}}}{\mathrm{d}t}}=\mathrm{\Delta }{q}_{\text{eva}}+\mathrm{\Delta }{q}_{\text{precip}}+{\mathit{\kappa }}_{\mathrm{a}}\cdot {\mathrm{\nabla }}^{\mathrm{2}}{q}_{\text{air}}-\mathit{u}\cdot \mathrm{\nabla }{q}_{\text{air}}+\mathrm{\Delta }{q}_{\text{correct}}.\end{array}$$

The main physical processes that control the surface temperature tendencies are solar (short-wave) and thermal (long-wave) radiation, the hydrological cycle (including evaporation, moisture transport and precipitation), horizontal transport of heat, and heat uptake in the subsurface ocean. GREB further simulates a number of diagnostic variables, such as precipitation snow/ice cover and sea ice, resulting from the simulation of the prognostic variables.

Atmospheric circulation (mean winds) and cloud cover are seasonally prescribed boundary conditions, and prescribed flux corrections (F_correct, $F_{{\mathrm{o}}_{\text{correct}}}$ and Δq_correct) are used to keep the GREB model close to the observed mean climate. State-independent flux corrections of surface temperatures or other variables allow a climate model to be close to those observed or any other state, while still being able to fully respond to external forcing or internal variability (Dommenget and Rezny, 2018; Irvine et al., 2013; Schneider, 1996). The flux correction terms are estimated by balancing the tendency in Eqs. (1)–(4) for observed boundary conditions to result in the observed T_surf, T_ocean​​​​​​​ and q_air for each calendar month (see Dommenget and Flöter, 2011, for details).

Since the GREB model does not simulate the atmospheric or ocean circulation, it is conceptually very different from coupled general circulation model (CGCM) simulations. The model does simulate important climate feedbacks such as the water vapor and ice–albedo feedback, but an important limitation of the GREB model is that the response to external forcing or model parameter perturbations does not involve circulation or cloud feedbacks. GREB does not have any internal (natural) variability since daily weather systems are not simulated. Subsequently, the control climate or response to external forcing can be estimated from one single year, assuming an equilibrium has been reached. The primary advantage of the GREB model in the context of this study is its simplicity, speed and low computational cost. The simulation of 1 year of global climate with the GREB model can be done in about 1 s (about 100 000 simulated years per day on a desktop computer), and model simplicity allows the user to straightforwardly investigate cause and effect in coupled simulations.

Table 2 Processes and their relevant numeric scheme for the ice sheet model.

3.2 Ice sheet model

The ice sheet model is a global thermomechanical ice flow model that comprises momentum balance, mass balance and energy balance modules with the prognostic variables: thickness and temperature, and diagnostic velocities. This subsection will describe the ice sheet model, including the model grid, dynamical methods used, parameterizations and approximations made. A short summary of the basic numerical schemes used are listed in Table 2.

3.2.1 Model grid

The ice sheet model uses the same horizontal grid as the GREB model. The Arakawa C scheme (Pollard and Deconto, 2012) is adopted for the simulation of velocities, with the ice thickness and temperature specified at the center of the grid, and zonal and meridional velocities are specified at the grid boundary midpoint. For the vertical coordinates, we apply a terrain-following coordinate, ξ, in the ice sheet model, where

$$\begin{array}{}\text{(5)}& \mathit{\xi }={\displaystyle \frac{z-\frac{H}{\mathrm{2}}}{\frac{H}{\mathrm{2}}}}.\end{array}$$

We chose the number of layers to be 4, to be close to the minimal number of layers which can still resolve the vertical velocity in the ice sheets: the surface layer (ξ=1), two Gaussian nodes ($\mathit{\xi }=\pm \frac{\mathrm{1}}{\sqrt{\mathrm{3}}}$, nodes for 2-point Gaussian quadrature; Hildebrand, 1987) and the base layer ($\mathit{\xi }=-\mathrm{1}$). The vertical integration in the model is based on Gauss–Jacobi quadrature (Hildebrand, 1987), where temperature vertical distribution is estimated by a polynomial curve fitting according to the four layers, which is expressed by

$$\begin{array}{}\text{(6)}& T\left(\mathit{\xi }\right)=c_{\mathrm{0}}+c_{\mathrm{1}}\mathit{\xi }+c_{\mathrm{2}}{\mathit{\xi }}^{\mathrm{2}}+c_{\mathrm{3}}{\mathit{\xi }}^{\mathrm{3}},\end{array}$$

where T is the temperature, c_i ($i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}$) are regression coefficients derived from the temperatures at the above four vertical nodes at each time step. The global, horizontal model grid has cyclic boundary conditions. For the grid points at the poles, we assume the poleward neighbor is the point at the same latitude, but shifted by 180^∘, following the approach in Allen et al. (1991). To avoid numerical instability in the polar regions, a zonal wave filter is applied from 76.875^∘ S to the South Pole (Lin and Rood, 1997; Suarez and Takacs, 1995).

3.2.2 Glacier mask

The GREB-ISM ice sheet evolution depends on whether the ice is grounded (land), floating (ice shelves) or if we have thin ice over the ocean (sea ice). The ice thickness, H, is used for both sea ice and ice sheet. For very thin ice cover, the gravity-driven ice flow is negligible and thus it does not follow ice sheet dynamics (e.g., snow or sea ice). To distinguish large ice mass from snow or sea ice, H must be above 10 m (Fyke et al., 2011). In detail, the points are as follows:

• Grounded ice (land) points. Ice sheet is grounded on bedrock, satisfying the condition (Larour et al., 2012):

  $$b+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{o}}}}H>\mathrm{0}.$$

• Floating ice (ice shelves) points. Ice thickness H ≥ 10 m and does not reach the bedrock, satisfying the floating condition:

  $$b+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{o}}}}H\le \mathrm{0}.$$

• Ocean points (all other points).. The ocean points here include sea ice grid (H≥0) as well.

The definition of this glacier mask does implicitly define groundling lines of glaciers by shifting points from grounded ice to floating ice according to the ice thickness, bed topography and global sea level (see also Sect. 3.3.8 for the sea level impact on the bed elevation).

3.2.3 Mass balance

The ice surface elevation, calculated from the mass balance equation, is the primary input from the ice sheet model to the GREB model, calculated for all global grid points. The mass balance equation is as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial H}{\partial t}}=s-a-\mathrm{\nabla }\cdot \left({\mathit{V}}_{\text{m}}H\right),\end{array}$$

where the accumulation of snow (s), ablation (melting) of ice (a) and ice transport ((∇⋅(V_mH)) control the mass balance. The surface mass balance terms (s,a) are calculated at the same time step as GREB (half-day), and thus we have seasonal ice thickness change. The ice transport term is calculated with an annual time step (Sect. 3.2.4).

The methods used to calculate the terms on the right hand side depend on whether ice is grounded (ice sheet), floating (ice shelves) or sea ice. The mass balance for sea ice is described in Sect. 3.3.4. For the ice sheet and shelves, the two local surface forcing terms for the ice mass balance from Eq. (7) are the source (accumulation) and sink (ablation) terms. The accumulation is due to snowfall:

$$\begin{array}{}\text{(8)}& s={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{o}}}{{\mathit{\rho }}_{\mathrm{i}}}}r\cdot p,\end{array}$$

with the snowfall ratio, r:

$$\begin{array}{}\text{(9)}& r=\left\{\begin{array}{ll}\mathrm{1},& T_{\text{atmos}}<T_{\mathrm{m}}\text{ and }{T}_{\text{surf}}<T_{\mathrm{m}}-\mathrm{2}{}^{\circ }\mathrm{C}\\ \frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}-\frac{T_{\text{surf}}-T_{\mathrm{m}}}{\mathrm{2}{}^{\circ }\mathrm{C}}\right),& T_{\text{atmos}}<T_{\mathrm{m}}\\ & \text{ and }{T}_{\mathrm{m}}-\mathrm{2}{}^{\circ }\mathrm{C}<T_{\text{surf}}<T_{\mathrm{m}}+\mathrm{2}{}^{\circ }\mathrm{C}\\ \mathrm{0},& \text{otherwise}.\end{array}\right.\end{array}$$

The ice ablation rate is due to surface melting by positive surface heat flux:

$$\begin{array}{}\text{(10)}& a=-{\displaystyle \frac{F_{\text{ice}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{m}}}},\end{array}$$

with the latent heat flux for melting ice, F_ice:

$$\begin{array}{}\text{(11)}& F_{\text{ice}}=\left\{\begin{array}{ll}{F}_{\text{surf}}& \text{partial\hspace{0.17em} melting}:F_{\text{surf}}\le F_{{max}_{\text{melt}}}\\ \frac{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{m}}H}{\mathrm{\Delta }t}& \text{complete\hspace{0.17em} melting}:F_{\text{surf}}>F_{{max}_{\text{melt}}}\\ \mathrm{0}& \text{no\hspace{0.17em} melting}:F_{\text{surf}}<\mathrm{0}.\end{array}\right.\end{array}$$

Here, the maximum heat flux for complete ice melting is $F_{{max}_{\text{melt}}}={\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{m}}\frac{H}{\mathrm{\Delta }t}$. The surface heat flux, F_surf, only considers the net surface heat flux beyond the freezing point:

$$\begin{array}{}\text{(12)}& F_{\text{surf}}={\mathit{\gamma }}_{\text{surf}}{\displaystyle \frac{T_{\text{se}}-T_{\mathrm{m}}}{\mathrm{\Delta }t}},\end{array}$$

and the estimated surface temperature without ice fusion is

$$\begin{array}{}\text{(13)}& T_{\text{se}}=T_{\mathrm{0}}+\mathrm{\Delta }t{\displaystyle \frac{F_{\text{net}}}{{\mathit{\gamma }}_{\text{surf}}}},\end{array}$$

where F_net is defined as the right hand side of Eq. (1).

We currently do not explicitly include an ocean basal melting scheme in our model. On the one hand, our ice shelf viscosity is tuned to fit the current day ice shelf thickness, which partially contains basal melting effects (see Sect. 3.2.4). On the other hand, including a basal melting scheme (Martin et al., 2011) does not contribute to a significant improvement in our model simulations (see Sect. 4.3).

Figure 3 GREB-ISM seasonal ice thickness (cm) during January–February–March (a, c) and July–August–September (b, d) from the coupled dynamic equilibrium experiment equilibrium state (200 kyr). The scale is chosen to highlight seasonal ice cover.

The snow accumulation and melting, as described above, control all land ice and snow cover and therefore also simulate the seasonal cycle of snow and ice cover over land. Figure 3 illustrates the seasonal cycle of ice cover in both hemispheres as simulated by the GREB-ISM model with present-day boundary conditions. The ice cover change for ocean points comes from sea ice changes, which is described in Sect. 3.3.4. The overall snow cover (land) distribution and seasonal cycle resemble observations (Robinson et al., 2012). Similarly, the mean sea ice extent and seasonal cycle are comparable with those observed (Rayner et al., 2003), with some overestimation of sea ice extent around Antarctica in summer.

Calving

A boundary condition for the mass transport equations is required at the ice front: here, ice from the ice sheet can be freely advected to the attached ocean grid and become sea ice (see Sect. 3.3.4 for the dynamics of sea ice). In this way, calving is diagnosed as transport from ground (land) or floating ice (shelves) onto ocean points.

3.2.4 Momentum balance

Ice flow on grounded ice points is solved based on the shallow ice approximation (SIA; Hutter, 1983; Morland, 1984) for momentum balance:

$$\begin{array}{}\text{(14)}& {\displaystyle }\mathit{V}={\mathit{V}}_{\mathrm{b}}-\mathrm{2}{\mathit{\rho }}_{\mathrm{i}}g\mathrm{\nabla }{z}_{\text{topo}}{\int }_{z_{\mathrm{b}}}^{z}A\exp \left({\displaystyle \frac{-Q}{RT}}\right){\mathit{\sigma }}_{\mathrm{e}}^{n-\mathrm{1}}\left(H-z^{\prime }\right)\mathrm{d}{z}^{\prime },\text{(15)}& {\displaystyle }{\mathit{V}}_{\text{m}}={\displaystyle \frac{\mathrm{1}}{z-z_{\mathrm{b}}}}{\int }_{z_{\mathrm{b}}}^z\mathit{V}\mathrm{d}{z}^{\prime },\end{array}$$

and the shallow shelf approximation (SSA; Macayeal, 1989) on floating ice points (solved in geo-coordinate latitude ϕ and longitude λ):

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\left({\mathit{\eta }}_{\text{SSA}}H\left(\mathrm{4}{\displaystyle \frac{\partial V_x}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}+\mathrm{2}{\displaystyle \frac{\partial V_y}{r_{\mathrm{e}}\partial \mathit{\varphi }}}\right)\right)\\ & +{\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\varphi }}}\left({\mathit{\eta }}_{\text{SSA}}H\left({\displaystyle \frac{\partial V_x}{r_{\mathrm{e}}\partial \mathit{\varphi }}}+{\displaystyle \frac{\partial V_y}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\right)\cos \mathit{\varphi }\right)\\ & ={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial z_{\text{topo}}}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}},\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\varphi }}}\left({\mathit{\eta }}_{\text{SSA}}H\left(\mathrm{4}{\displaystyle \frac{\partial V_y}{r_{\mathrm{e}}\partial \mathit{\varphi }}}+\mathrm{2}{\displaystyle \frac{\partial V_x}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\right)\cos \mathit{\varphi }\right)\\ & +{\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\left({\mathit{\eta }}_{\text{SSA}}H\left({\displaystyle \frac{\partial V_x}{r_{\mathrm{e}}\partial \mathit{\varphi }}}+{\displaystyle \frac{\partial V_y}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\right)\right)\\ & ={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial z_{\text{topo}}}{r_{\mathrm{e}}\partial \mathit{\varphi }}}.\end{array}\end{array}$$

The viscosity (η_SSA) in our model is larger than in other models (Bueler and Brown, 2009). This high viscosity is tuned by adjusting ice shelf thickness to observation, which may be impacted by uncertainties in the observations, other model fields or in physical processes such as ice shelf basal melting effects.

Vertical velocities are recovered through incompressibility:

$$\begin{array}{}\text{(18)}& w=-{\int }_{z_{\mathrm{b}}}^z\mathrm{\nabla }\cdot \mathit{V}\mathrm{d}z.\end{array}$$

The deformation of ice under stress is described by Glen's flow law (Glen, 1953, 1954, 1955):

$$\begin{array}{}\text{(19)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}EA{\mathit{\sigma }}_{\mathrm{e}}^{n-\mathrm{1}}}},A=A_{\mathrm{0}}\exp \left({\displaystyle \frac{-Q}{R{T}^{\prime }}}\right),\end{array}$$

where T^′ is the temperature corrected for the dependence of melting point on pressure:

$$\begin{array}{}\text{(20)}& T^{\prime }=T-\mathit{\beta }(H-z).\end{array}$$

In our model, the viscosity η_SSA (Table 1) has been set as a constant value to match with the observed ice surface velocity and calving in the stand-alone dynamic equilibrium experiment (Sect. 4.3). Each of Eqs. (14)–(18) above are expressed in z coordinates but are transformed into ξ coordinates for the model integration. Boundary conditions for the mechanical model are required at the ice sheet surface, base and at the ice shelf-ocean front. A stress-free ice surface is assumed:

$$\begin{array}{}\text{(21)}& \mathit{\sigma }\cdot \mathit{n}=\mathrm{0},\end{array}$$

where n is the normal unit vector at the ice surface.

At the base, the horizontal ice velocities follow the viscous-type sliding law defined in Greve (1997):

$$\begin{array}{}\text{(22)}& {\mathit{V}}_{\mathrm{b}}=-C_{\text{sl}}H\mathrm{|}\mathrm{|}\mathrm{\nabla }{z}_{\text{topo}}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}\mathrm{\nabla }{z}_{\text{topo}},z=z_{\mathrm{b}}.\end{array}$$

The value for C_sl is as in Greve (1997). In Sect. 4.3 we discuss to what extent variations in C_sl could improve the simulations.

The stress conditions for the horizontal ice shelf velocities at the interface with the open-ocean points follow Greve and Blatter (2009), which in our model is expressed as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle }\begin{array}{rl}& \mathrm{4}{\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\left({\mathit{\eta }}_{\text{SSA}}H{\displaystyle \frac{\partial V_x}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\right)\\ & +\mathrm{2}{\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\left({\mathit{\eta }}_{\text{SSA}}H{\displaystyle \frac{\partial V_y}{r_{\mathrm{e}}\partial \mathit{\varphi }}}\right)={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial z_{\text{topo}}}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}},\end{array}\text{(24)}& {\displaystyle }\begin{array}{rl}& \mathrm{4}{\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\varphi }}}\left({\mathit{\eta }}_{\text{SSA}}H{\displaystyle \frac{\partial V_y}{a\partial \mathit{\varphi }}}\cos \mathit{\varphi }\right)\\ & +\mathrm{2}{\displaystyle \frac{\partial }{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\varphi }}}\left({\mathit{\eta }}_{\text{SSA}}H{\displaystyle \frac{\partial V_x}{r_{\mathrm{e}}\cos \mathit{\varphi }\partial \mathit{\lambda }}}\cos \mathit{\varphi }\right)\\ & ={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial z_{\text{topo}}}{r_{\mathrm{e}}\partial \mathit{\varphi }}}.\end{array}\end{array}$$

3.2.5 Energy balance

The ice temperature (energy) balance is calculated as follows:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\displaystyle \frac{\partial T}{\partial t}}& =-\mathit{V}\cdot \mathrm{\nabla }T-w{\displaystyle \frac{\partial }{\partial z}}T+{\displaystyle \frac{\partial }{\partial z}}{\displaystyle \frac{\mathit{\kappa }}{{\mathit{\rho }}_{\mathrm{i}}{C}_p}}{\displaystyle \frac{\partial }{\partial z}}T\\ & +{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{i}}{C}_p}}({\mathit{\sigma }}_{xz},{\mathit{\sigma }}_{yz})\cdot {\displaystyle \frac{\partial \mathit{V}}{\partial z}}.\end{array}\end{array}$$

The ice temperature balance at the surface is constrained by T_surf as computed in the GREB-ISM model (see Sect. 3.3.1):

$$\begin{array}{}\text{(26)}& T=T_{\text{surf}},z=z_{\text{topo}}.\end{array}$$

The geothermal heat flux is an important boundary condition for ice sheets. Previous studies show that the model with uniform geothermal heat flux is still able to reproduce the ice sheet evolution in the paleoclimate (Abe-Ouchi et al., 2007; Tigchelaar et al., 2019). For consistency, we therefore assume a globally constant bottom layer geothermal flux as in Huybrechts et al. (1996) and Payne et al. (2000):

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial T}{\partial z}}=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}{C}_{p}G}{\mathit{\kappa }}},z=z_{\mathrm{b}}.\end{array}$$

3.3 Coupling of the GREB model to the ice sheet

The introduction of an ice sheet model requires a number of changes to the original GREB model. In the following, we describe the changes made to the GREB model equations and illustrate how they affect the simulation of the GREB climate.

3.3.1 Energy exchange between GREB and ice sheet/sea ice

The introduction of a prognostic ice sheet model introduces the additional heat flux term, F_ice, for the T_surf tendency Eq. (1), resulting in the new equation:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{\mathit{\gamma }}_{\text{surf}}{\displaystyle \frac{\mathrm{d}{T}_{\text{surf}}}{\mathrm{d}t}}=& F_{\text{solar}}+F_{\text{thermal}}+F_{\text{latent}}+F_{\text{sense}}+F_{\text{ice}}\\ & +F_{\text{ocean}}+F_{\text{correct}}.\end{array}\end{array}$$

The calculations of F_ice are described in Sect. 3.2.3 (mass balance) and Sect. 3.3.4 (sea ice)​​​​​​​. The effect of F_ice can best be illustrated by a simple response experiment in which we add a 10 m ice cover and evaluate how surface temperature responds to it (Fig. 4). In this response experiment 10 m of ice cover is introduced over a large region of Europe (Fig. 4d, black box) at the start of the simulation, and then the fully coupled GREB-ISM model is run for 4 years to respond to this change.

Figure 4 GREB-ISM response to adding a 10 m ice sheet in surface temperature (units: ^∘C) and ice thickness (units: m). In (a)–(c) are the temperature and ice thickness evolution at three different locations. The black, red and blue curve represent control run surface temperature (without adding 10 m ice), scenario run surface temperature (with adding 10 m ice) and scenario run ice thickness. Panel (d) shows the temperature difference (units: ^∘C) between scenario and control at the end of the first simulation year. The black outlined region in (d) mark the area in which the initial 10 m ice sheet is added.

The introduction of the ice cover forces surface temperature below the freezing point at all locations, as long as the ice sheet is present (Fig. 4a–c). The atmospheric heat fluxes and sea ice dynamics force the sea ice to melt, which it does faster over the ocean points due to horizontal sea ice transport. Over land the ice cover melts after the first year and allows surface temperature to go back to the control run values. The atmospheric heat and moisture transport cause cooling in adjacent regions (Fig. 4d).

3.3.2 Surface heat capacity

The surface layer effective heat capacity (γ_surf) in the GREB model is equal to the heat capacity of a water column of the mixed layer depth over ice-free ocean points and equivalent to 2 m soil for all other points (e.g., land and ice covered). Thus, the formation of sea ice changes the heat capacity from that of the mixed layer depth to a 2 m soil column. This is unchanged from the original GREB model.

3.3.3 Precipitation correction

The hydrological cycle model in GREB developed in Stassen et al. (2019) simulates precipitation as a function of the simulated atmospheric humidity (q_air), the observed mean and standard deviation of the vertical air motion (ω_mean, ω_SD):

$$\begin{array}{}\text{(29)}& \begin{array}{rl}\mathrm{\Delta }{q}_{{\text{precip}}_{\text{S2019}}}& =r_{\text{precip}}\cdot q_{\text{air}}\cdot (c_{\text{rq}}\cdot \text{rq}+c_{\mathit{\omega }}\cdot {\mathit{\omega }}_{\text{mean}}\\ & +c_{\mathit{\omega }\text{SD}}\cdot {\mathit{\omega }}_{\text{SD}}).\end{array}\end{array}$$

This model aimed at a realistic simulation of precipitation with a focus on the regions of greatest precipitation, i.e., the tropical oceans. While the precipitation model is very good in these regions (Stassen et al., 2019), it only has limited skills over higher-latitude land regions, which are most important for the ice sheet mass balance of the GREB-ISM.

To allow the ice sheet mass balance to receive unbiased mean precipitation forcing under present-day conditions, we introduced a land precipitation correction in the GREB-ISM model. The new precipitation equation with flux correction is expressed as follows:

$$\begin{array}{}\text{(30)}& \mathrm{\Delta }{q}_{\text{precip}}=\mathrm{\Delta }{q}_{{\text{precip}}_{\text{S2019}}}+q_{\text{zonal}}\cdot p_{\text{correct}},\end{array}$$

where q_zonal⋅p_correct is the flux correction of the equation. The flux corrections are only active over land and are a function of calendar month. They are estimated in a way that the simulated Δq_precip matches the precipitation data in Sect. 2 for every calendar month of the year.

Here, we note that the $\mathrm{\Delta }{q}_{{\text{precip}}_{\text{S2019}}}$ model assumes that precipitation is proportional to the local humidity (q_air). Stassen et al. (2019) demonstrate that this assumption is less appropriate in higher-latitude land regions, as there is no clear relationship between the local q_air and Δq_precip. Due to the lack of a clear local relationship, we relaxed this constraint and assumed that the precipitation over land is a function of the zonal mean humidity, reflecting the mostly zonal structure of the atmospheric circulation. Therefore we set the correction term to be proportional to the zonal mean q_air defining q_zonal. Within 30^∘ of the poles q_zonal is estimated as the mean from the pole to 60^∘.

With this approach the precipitation over higher-latitude land responds to cooling or warming similarly to other regions (e.g., oceans for lower latitudes). We will discuss the precipitation response of the GREB-ISM further below in the context of the response experiments.

3.3.4 Sea ice

Sea ice is a diagnostic variable in the original GREB model but is now changed to be a prognostic variable in GREB-ISM. Over land and ice shelf points, ice thicknesses (H) follow the dynamics described in the ice sheet model Sect. 3.2. Over ocean points we use the same prognostic variable (H), but the sea ice thickness dynamics follow a different tendency equation, namely

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\partial H}{\partial t}}=\mathrm{\Delta }{H}_{\text{seaice}}-{\mathit{\kappa }}_{\text{si}}{\mathrm{\nabla }}^{\mathrm{2}}H,\end{array}$$

with the local sea ice growth,

$$\begin{array}{}\text{(32)}& \mathrm{\Delta }{H}_{\text{seaice}}={\displaystyle \frac{-F_{\text{ice}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{m}}}},\end{array}$$

and where the latent heat of ice fusion F_surf is defined by Eqs. (11)–(13):

$$\begin{array}{}\text{(33)}& \begin{array}{rlrlrl}& F_{\text{ice}}=F_{\text{surf}}& & \text{ice grows:}& & T_{\text{se}}<T_{\text{sm}},F_{\text{surf}}<\mathrm{0}\\ & & & & & \text{and }H<\mathrm{0.5}\mathrm{m},\\ & F_{\text{ice}}\text{ from Eq. (11)}& & \text{ice melts:}& & T_{\text{se}}>T_{\text{sm}},F_{\text{surf}}>\mathrm{0},\\ & F_{\text{ice}}=\mathrm{0}& & \text{no change:}& & \text{otherwise}.\end{array}\end{array}$$

The sea ice growth threshold of 0.5 m reflects the fact that sea ice is a very good insulator and subsequently does not transfer atmospheric heat fluxes very well once a certain ice thickness is reached. This in practice limits the growth of sea ice by atmospheric heat flux to less than 0.5 m typically. In this case F_ice=0 and it will no longer grow the sea ice, but only cool T_surf (Eq. 28).

Sea ice transport is estimated by isotropic diffusion (κ_si∇²H). This approximates the effect of turbulent winds and ocean currents transporting sea ice, leading to fast decay of sea ice near open ocean. The diffusion coefficient κ_si was chosen to roughly lead to a sea-ice-decaying timescale of about 1 month.

3.3.5 Albedo coupled to ice sheet

The surface albedo (α_surf) in the original GREB model was diagnosed as function of T_surf but is now diagnosed as a function of the ice thickness (H):

$$\begin{array}{}\text{(34)}& \begin{array}{rlrl}& {\mathit{\alpha }}_{\text{surf}}=\mathrm{0.1}& & H=\mathrm{0.0},\\ & {\mathit{\alpha }}_{\text{surf}}=\mathrm{0.1}+\mathrm{17.5}{\mathrm{m}}^{-\mathrm{1}}\cdot H& & H\in [\mathrm{0.0},\mathrm{0.02}\mathrm{m}],\\ & {\mathit{\alpha }}_{\text{surf}}=\mathrm{0.45}& & H>\mathrm{0.02}\mathrm{m}.\end{array}\end{array}$$

The linear relation between ice thickness and albedo in the GREB-ISM model was estimated from the assumption that for the observed northern hemispheric seasonal cycle of snow/ice cover over land the overall albedo matches the mean overall albedo of the original GREB model.

3.3.6 Topography coupled to ice sheet

The land topography (z_topo) in the original GREB model is a fixed boundary condition that influences a number of processes: thermal radiation, hydrological cycle, and the transport of heat and moisture by advection and diffusion. For GREB-ISM the land topography is now a function of the bed topography and ice sheet height:

$$\begin{array}{}\text{(35)}& \begin{array}{rl}& z_{\text{topo}}=b+H,\text{ for grounded ice},\\ & z_{\text{topo}}=\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{o}}}}\right)H,\text{ for floating ice}.\end{array}\end{array}$$

The GREB-ISM does not simulate any glacial isostatic adjustment.

3.3.7 Sensible heat flux between surface and atmosphere

The variable land topography (z_topo) should affect the sensible heat flux between T_surf and T_atmos, which was not simulated in the original GREB model. Here it needs to be considered that the GREB model does not resolve the vertical structure of the atmosphere, as it only has one atmospheric layer. However, in the real world T_atmos decreases with surface elevation, following a moist adiabatic lapse rate. We therefore change the sensible heat flux between T_surf and T_atmos, which was approximated in the original GREB model by Newtonian coupling between T_surf and T_atmos. In the GREB-ISM model this is now replaced with a Newtonian coupling between T_surf and an adjusted T_atmos:

$$\begin{array}{}\text{(36)}& F_{\text{sense}}={\text{ct}}_{\text{sense}}(T_{\text{atmos}}+\mathrm{\Gamma }\cdot z_{\text{topo}}-T_{\text{surf}}).\end{array}$$

Figure 5 GREB-ISM response to a lifting of the topography by 1000 m for surface temperature (a, units: ^∘C) and precipitation (b, units: mm d⁻¹). The response is defined as the scenario run (1000 m topography lifting) minus control run (no lifting) at the end of the first simulation year. The box represents the lifted region.

Here we choose a globally constant moist adiabatic lapse rate Γ = −6 K km⁻¹. The effect of this sensible heat flux is illustrated with a simple response experiment; see Fig. 5. For this experiment we increase z_topo and precipitation relative to a control simulation with no changes in z_topo (Fig. 5). T_surf decreases in response to the topographic perturbation, approximately linearly to the moist adiabatic lapse rate. The higher topography also affects the hydrological cycle, reducing the precipitation locally and also remotely through transport of relatively reduced atmospheric humidity.

3.3.8 Sea level and land–sea mask

A sea level subroutine is added in GREB-ISM. Only grounded ice thickness impacts the global sea level. Consequently, the sea level change slv is defined by

$$\begin{array}{}\text{(37)}& \text{slv}={\displaystyle \frac{{\int }_{\text{grounded}}(H-H_{\text{ref}})\mathrm{d}A}{A_{\text{ocean}}}},\end{array}$$

where H_ref is the reference ice thickness, A_ocean is total area of ocean grid and ∫_groundeddA is an integration over all grounded ice points. slv will be added to bed topography b, which eventually impacts the land–sea mask. The sea level and land–sea mask are updated every model year.

The soil moisture, which is a boundary condition for estimating surface evaporation, is initially set to observed values over land and then changes if land–sea distribution alters. If the sea level lowers and an ocean point turns into a land point (b > 0) then the land point has a soil moisture value of 0.3 (equivalent to the mean value for land points in Dommenget and Flöter, 2011). In turn, if the sea level rises and a land point turns into an ocean point (b < 0), then the soil moisture value is set to 1.0.

3.3.9 Meridional heat transport

The study by Dommenget et al. (2019) showed that the GREB model, without flux corrections for T_surf, has a high-latitude climate that is too cold and a tropical climate that is too warm, indicating that the meridional heat transport is too weak. The meridional heat transport in the GREB model results from the atmospheric heat transport by the mean advection due to the mean horizontal wind field and by isotropic diffusion. The latter depends on the diffusion coefficient ${\mathit{\kappa }}_{\mathrm{a}}=\mathrm{8}\times {\mathrm{10}}^{\mathrm{5}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in the GREB model. This value is not strongly constrained by observations and may effectively be different by an order of magnitude. Since the meridional heat transport may play an important role in the global ice-age cycle, we enhance this diffusion coefficient by a factor of 5. This reduces the mean T_surf bias in higher latitudes and the tropics in the GREB model without flux corrections, while at the same time does not increase biases in other locations, indicating it is a better approximation of the isotropic diffusion.

4 Model benchmark: ice sheet model stand-alone simulations

We start our evaluation of the new ice sheet model GREB-ISM with stand-alone ice sheet model simulations forced with idealized or observed boundary conditions. These simulations focus on the ice sheet simulation only. Section 4.1 and 4.2 use standard experiments from the European Ice Sheet Modelling Initiative (EISMINT) model intercomparison Phase I (Huybrechts et al., 1996) and II (Payne et al., 2000), which test the ice sheet model response to idealized mass and temperature forcing within a given horizontal resolution, with the ice mechanics decoupled from the thermodynamics in EISMINT I and coupled in EISMINT II. In Sect. 4.3, we discuss a simulation on the global GREB-ISM grid forced with observed boundary conditions to estimate the dynamically forced equilibrium of the ice sheet model. Finally, we discuss an idealized time-varying ice sheet response experiment, forced with temperature and precipitation similar to (Niu et al., 2019) over the past 250 kyr.

Table 3 Variables (upper) and parameters (below) list for EISMINT experiments.

Table 4 EISMINT I steady state experiment result comparison between GREB-ISM and the model ensemble from H96 for fixed-margin (F) and moving-margin (M) experiments.

4.1 EISMINT I

All simulations in EISMINT I (Huybrechts et al., 1996, H96 hereafter) are based on a regional grid in Cartesian coordinates that have higher resolutions than the GREB model grid (∼ 50 km). For a better comparison of the numerical schemes we changed the GREB-ISM grid (3.75^∘ × 3.75^∘) for these experiments to a model grid with 96 points in the zonal and 144 points in the meridional direction (3.75^∘ × 1.25^∘). Only the first 15 points in the meridional direction are used for the ice sheet simulation. The ice sheet divide in these simulations is the South Pole, and the length of the meridional grid is 50 km. The simulations are integrated for 200 kyr, but near equilibrium is reached after about 50 kyr.

The mass balance S and surface temperature T_surf forcings are given as follows:

$$\begin{array}{}\text{(38)}& {\displaystyle }\text{Fixed margin experiment:}\left\{\begin{array}{l}S=\mathrm{0.3}\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}\\ T_{\text{surf}}=T_{min}+S_{\mathrm{T}}{d}^{\mathrm{3}},\end{array}\right.\text{(39)}& {\displaystyle }\begin{array}{rl}& \text{Moving margin experiment:}\\ & \left\{\begin{array}{l}S=min\mathit{\{}{S}_{max},S_{\mathrm{b}}(R_{\text{el}}-d\left)\mathit{\right\}}\\ T_{\text{surf}}=(\mathrm{270}\mathrm{K}-S_{\mathrm{H}}H).\end{array}\right.\end{array}\end{array}$$

The parameters in Eqs. (38) and (39) are listed in Table 3. Table 4 shows the comparison between the new ice sheet model GREB-ISM and model results from H96. The GREB-ISM simulations of the ice thickness at divide and mass flux at midpoint are mostly similar to those found in H96 for both the fixed and moving margin experiments. The ice mass flux in the GREB-ISM is larger than in H96 for the moving margin experiment. An additional experiment (not shown) with the GREB-ISM in Cartesian coordinates as used in the EISMINT I simulation finds the ice mass flux close to H96, suggesting this result may be dependent on mesh shape.

Figure 6 Time evolution of ice thickness (a, c, unit: m) and homologous basal temperature (b, d, unit: K) in the EISMINT I fixed (a, b) and moving (c,  d) margin experiments with GREB-ISM with 20 and 40 kyr period forcing. R marks the range (maximum minus minimum in the last 50 kyr) of the simulated variables.

The transition experiments with oscillating forcing of temperature and mass balance with periods of 20 and 40 kyr are presented in Fig. 6. The GREB-ISM ice thickness simulation is similar to those of H96 for both fixed and moving margin experiments (Fig. 6). In both experiments, the basal temperature at the divide is about 1–2 ^∘C colder than in the H96 simulations, which is related to the coarse vertical resolution. This mismatch disappears if we increase the vertical resolution to 10 layers (not shown).

Table 5 Results for basic glaciological quantities in EISMINT II experiments after 200 kyr. Differences are defined as current experiment minus experiment A. Percentage changes are relative to experiment A. The results of P2000 are shown in the form of “mean ± range”. See text for details.

Figure 7 Time evolution of total ice volume (a, b, units: 10⁶ km³) and ice calving (c, d, units: km³ yr⁻¹) in Greenland (a, c) and Antarctica (b, d) from the forced stand-alone dynamic equilibrium experiment.

4.2 EISMINT II

EISMINT II experiments (Payne et al., 2000, P2000 here after) involve coupling between the mechanical and thermodynamical components of the ice sheet model. These experiments are designed to test how the ice sheet temperature variations interact with the ice sheet transport. The GREB-ISM model grid used is similar to in EISMINT I, but the number of points in the meridional direction is increased from 15 to 31 and the length of the meridional grid is set to 25 km. All experiments are integrated for 200 kyr. The boundary conditions for the first experiment (A) are as follows:

$$\begin{array}{}\text{(40)}& \left\{\begin{array}{l}S=min\mathit{\{}{S}_{max},S_{\mathrm{b}}(R_{\text{el}}-d\left)\mathit{\right\}}\\ T_{\text{surf}}=T_{min}+S_{\mathrm{T}}d,\end{array}\right.\end{array}$$

with the parameters given in Table 3. The results of experiment A are summarized in Table 5. The final GREB-ISM values for ice volume, area, divide thickness and basal temperature at the ice sheet divide are all within the range of the models in P2000, indicating a fairly good agreement. The basal melt fraction is underestimated by the GREB-ISM by about 30 %, which is related to a cold bias at the bed of the ice sheet.

Experiment B and C in EISMINT II are designed for testing the model sensitivity to various boundary conditions. T_min in experiment B is set as 5 K cooler than in experiment A, to evaluate the sensitivity of the model to the mean ice temperature. Table 5 depicts the difference between experiment B and A. The GREB-ISM shows, in general, similar changes in ice volume, ice divide thickness and ice divide basal temperature as in P2000. However, the basal melt fraction change shows a significant discrepancy, which is related to the cold bias of the basal temperature in experiment A.

For experiment C, S_max and R_el are set as 0.25 m yr⁻¹ and 425 km, respectively, to evaluate the impact of different mass balances. The results of experiment C are shown in Table 5. For the changes in ice volume, area, divide thickness and divide basal temperature, the response difference between experiment C and A in GREB-ISM is equivalent to results from P2000. The changes in melt fraction in the GREB-ISM deviate from those of P2000, which is again likely to be related to the cold bias in basal temperatures in the GREB-ISM in experiment A.

Overall, the model reproduces the total ice thickness and ice cover well in the idealized experiments of EISMINT I and II. Although there is a bias in the basal temperature estimation in GREB-ISM, this issue does not have a significant impact on the ice thickness and cover area, which suggests the model is appropriate for global climate and ice evolution simulations.

4.3 Globally forced dynamical equilibrium

We now focus on simulating the observed global ice sheets forced with present-day boundary conditions. Although we cannot assume that observed Greenland and Antarctic Ice Sheets are in equilibrium with present-day forcing, the dynamic equilibrium simulation should produce a global ice sheet distribution similar to the current observations.

Figure 8 Comparison of ice surface velocity (unit: m yr⁻¹) from observations (left) and the GREB-ISM forced stand-alone dynamic equilibrium experiment at equilibrium state (right).

Figure 9 Results from the GREB-ISM forced stand-alone dynamic equilibrium simulation at equilibrium state: annual mean ice thickness (a, c) and the ice thickness difference (b, d) between GREB-ISM simulation and the observation in Greenland (a, b) and Antarctica (c, d). The ice thickness observation is derived from the BedMachine dataset (Morlighem et al., 2017, 2020).

Table 6 Ice volume and boundary calving from the forced dynamic equilibrium experiment and observation.

Ice surface temperature and precipitation forcings in the experiment are set to the climatologies derived from ERA-interim, NCEP-DOE and GPCP data. GREB-ISM is run for 200 kyr, initialized with observed ice thickness. Figures 7–9 show results from this simulation, and Table 6 compares the simulation values of total ice volume boundary calving with observed values from the literature.

The model reaches an equilibrium after about 50 kyr for both the Northern and Southern Hemispheres. Greenland ice thicknesses and calving rates show only small differences compared with the initial values. They are also within the estimated calving values from observation (Bigg, 1999). The trends in Antarctica are larger, in particular over West Antarctica. Here we see a significant increase in ice volume and calving (Figs. 7d and 9d). The West Antarctic ice sheet thickness increase is inconsistent with the observed values, suggesting a model limitation.

We could not find the specific limitation that is causing West Antarctic Ice Sheet bias. The precipitation forcing does play a role in controlling the West Antarctic Ice Sheet, but we could not find any reasonable precipitation forcing that would result in significantly improved simulations of the West Antarctic Ice Sheet. The parameterization of the floating ice for ice shelves (SSA) also impacts the simulation of West Antarctic Ice Sheet. The ice shelf can grow and become grounded as an ice sheet with lower viscosity. However, again we could not find any reasonable value for the ice viscosity (η_SSA) that would significantly reduce this bias. We further tested different sliding law coefficient C_sl, ranging from 6 × 10³ to 6 × 10⁵ yr⁻¹. The result indicates that the varying coefficient values do not bring a fundamental simulation improvement. Similarly, a basal melting scheme (Martin et al., 2011) with different strength has also been tested, but improvement could not be found.

The simulated ice surface velocity for Antarctica and Greenland shows a reasonable pattern, capturing the main features of the transport (Fig. 8) and the mean values. For Antarctica the ice mean flow is 109 m yr⁻¹, faster than the observation (80 m yr⁻¹) from MEaSUREs data (Mouginot et al., 2012; Rignot et al., 2011, 2017), and slower in the interior and faster near the boundaries. The largest velocities (more than 1000 m yr⁻¹) appear in ice shelf regions (Ross and Filchner–Ronne Ice Shelf), which is due to the presence of the floating ice for ice shelves (SSA). Similarly, Greenland ice velocities are also in good agreement with observations (Joughin, 2017; Joughin et al., 2010) in terms of pattern and mean flow magnitude (57 simulated and 56 m yr⁻¹ observed).

4.4 Transition experiment

We next evaluate the capability of the global ice sheet model to respond to realistic changes in the boundary conditions. We therefore design an experiment, in which we force the GREB-ISM with surface temperature and precipitation over the past 250 kyr, similar to the one discussed in Niu et al. (2019) for the Northern Hemisphere, but extended to the whole globe to evaluate the response of the ice sheet on a global scale. The surface temperature and precipitation forcing for this experiment are as follows:

$$\begin{array}{}\text{(41)}& \left\{\begin{array}{ll}{T}_{\text{surf}}(\mathit{\lambda },\mathit{\varphi },t)=& T_{\text{today}}(\mathit{\lambda },\mathit{\varphi },t_{\text{day}})+\left(T_{\text{LGM}}\right(\mathit{\lambda },\mathit{\varphi },t_{\text{day}})\\ & -T_{\text{today}}(\mathit{\lambda },\mathit{\varphi },t_{\text{day}}))\frac{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}\left(t\right)-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{PD}}}{{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{LGM}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{PD}}}\\ S(\mathit{\lambda },\mathit{\varphi },t)=& min\left[S_{\text{today}}(\mathit{\lambda },\mathit{\varphi },t_{\text{day}})+\left(S_{\text{LGM}}\right(\mathit{\lambda },\mathit{\varphi },t_{\text{day}})\right.\\ & \left.-S_{\text{today}}(\mathit{\lambda },\mathit{\varphi },t_{\text{day}}))\frac{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}\left(t\right)-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{PD}}}{{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{LGM}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{PD}}},\mathrm{0}\right].\end{array}\right.\end{array}$$

The surface temperature (T_surf) and ice mass balance (S) are present-day regional and seasonally varying climatologies (T_today, S_today) plus a seasonally changing (t_day) forcing pattern for T_surf and S that varies according to δ¹⁸O proxy data derived from the Greenland Ice Core Project (GRIP) dataset (Greve, 1997). ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{PD}}$ and ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{LGM}}$ represent δ¹⁸O at the present day and Last Glacial Maximum (LGM), respectively. The LGM reference climate forcing pattern $T_{\text{LGM}}(\mathit{\lambda },\mathit{\varphi },t_{\text{day}})$ is taken from the AWI Earth system model (AWI-ESM) (see Sect. 2 for detail), which results from a CGCM simulation forced by insolation, greenhouse gas and ice sheet. The main feature of this forcing pattern (not shown) is a much colder climate (more than 10 ^∘C) from North America to Central Asia and Antarctica, which contained or were surrounded by large ice sheets (Kageyama et al., 2017). The simulation is integrated between −250 kyr and the present and initialized with present-day observed ice thickness.

Figure 10 Time series of simulated sea level (left axis; units: m) from the stand-alone transition experiment and δ¹⁸O proxy data (right axis, the axis has been inverted).

The time series in Fig. 10 depicts the sea level change in this simulation from −200 kyr BP compared with a δ¹⁸O proxy time series from ocean sediments (Imbrie et al., 1984). The two curves show similar time series variations with a correlation of −0.67. This indicates that qualitatively the GREB-ISM ice sheet shows similar overall global ice sheet variations to those observed over the past 200 kyr. The GREB-ISM sea level varies by about 120 m, which is the same as observed sea level changes (Fairbanks, 1989; Lambeck et al., 2014), indicating that the simulated ice sheet volume variations are similar to those observed. The sea level is also 20 m lower than the present day due to the excess West Antarctic Ice Sheet volume that we also observed in the dynamical equilibrium simulation.

Figure 11 Global ice thickness (unit: m) distribution in the Last Interglacial (a, b), the Last Glacial Maximum (c, d) and present day (e, f) from the stand-alone transition experiment.

There are several significant extremes in the past 200 kyr simulation, which correspond to the Last Interglacial (LIG; −127 kyr), Last Glacial Maximum (LGM; −21 kyr) and the present day. The ice sheet thicknesses for these three time periods are shown in Fig. 11. During the LIG, only the Greenland Ice Sheet thickness exceeded 400 m in the Northern Hemisphere, and the Antarctic Ice Sheet thickness is similar to present day. During the LGM large European (e.g., Fennoscandia) and North American (Laurentide) ice sheets are reproduced with thousands of meters ice thickness, which is also what we expected according to previous studies (Clark et al., 2009; Velichko et al., 1997).

Table 7 Annual mean ice volume in the stand-alone transition experiment for different time periods from GREB-ISM simulation and from Fyke et al. (2011).

The estimate of ice sheet volume in Greenland and Antarctica for the Last Interglacial, Last Glacial Maximum and Late Holocene from GREB-ISM and from Fyke et al. (2011) are presented in Table 7. Overall, our simulation of the Greenland Ice Sheet is similar to Fyke et al. (2011) but with larger time variations. However, the simulation of Antarctica ice thickness shows very little to no variation between these three periods. The difference between the GREB-ISM model and Fyke et al. (2011) in Antarctica ice sheet may be due to different experimental setups. Fyke et al. (2011) varied and changed the ice shelf parameterization periods during their simulation, which was not done in our experiments. In summary, the results of this experiment indicate that the GREB-ISM ice sheet model does have realistic responses to time-varying boundary conditions.

5 Model benchmark: GREB-ISM coupled simulations

We now focus on the fully coupled GREB-ISM model, in which the ice sheet and other climate variables are interacting in both directions. In the following sections, two sets of experiments are presented. First a dynamic equilibrium experiment is conducted, which is similar to the experiment discussed in Sect. 4.3, but now fully coupled with fixed boundary conditions. Second, a set of experiments with shortwave radiation oscillating on periods of 20, 50 and 100 kyr for the Northern Hemisphere are conducted. Those two experiments are designed to evaluate how coupling influences the model's behavior and to what extent the ice sheet responds to periodic solar forcing. The discussion of these experiments will focus on the introduction of the GREB-ISM model. A more detailed analysis of the ice sheet dynamics coupled with climate dynamics is left for future studies.

Figure 12 Results from the fully coupled dynamic equilibrium experiment: evolution of global annual mean surface temperature (a, units: ^∘C), total ice volume (b, units: 10⁶ km³), annual mean precipitation (c, units: mm d⁻¹) and sea level change (d, units: m). The dashed lines are modern observation references.

Figure 13 Same as Fig. 9 but for the coupled dynamic equilibrium experiment.

5.1 Dynamic equilibrium for present-day conditions

In this experiment, the GREB-ISM model is fully coupled and forced with the fixed boundary conditions of present-day 340 ppm CO₂ concentration and solar radiation. T_surf and land precipitation are flux-corrected to the mean present-day values. However, those flux-corrected variables can respond to changes in the climate system, since the flux correction terms are state-independent (see Sect. 3.1). The simulation is 200 kyr long and results are shown in Figs. 12 and 13.

T_surf and precipitation show no long-term drift and are close to the observation (Fig. 12a and c). Both reach equilibrium after about 50 kyr. The global ice volume difference is mainly contributed by ice thickness difference in the Southern Hemisphere (Fig. 12b), which is similar to the one in the forced experiment discussed in Sect. 4.3 (Figs. 7 and 9). As the ice volume increases, the sea level shows a clear decrease tendency and reaches equilibrium after 50 kyr as well. The ice thickness spatial pattern in the coupled experiment is comparable to the stand-alone experiment (Figs. 13 and 9). Overall, this control run simulation shows that the coupled GREB-ISM system converges towards an equilibrium state close to the observed one. The simulated trends appear to be mostly due to the anomalous growth of the West Antarctic Ice Sheet.

Figure 14 Time evolution of change in total ice volume (black, unit: 7 × 10⁶ km³), surface temperature (red, unit: ^∘C), precipitation (blue, unit: 10⁻¹ mm d⁻¹), ice cover area (cyan, unit: 4 × 10⁶ km²) and solar radiation at 65^∘ N (orange, unit: 4 W m⁻²) from the shortwave oscillation experiment in Northern (upper) and Southern (lower) Hemisphere with forcing period of 20 kyr (a, d), 50 kyr (b, e) and 100 kyr (c, f). The control equilibrium state values from the coupled dynamic equilibrium experiment are removed to obtain changes.

Figure 15 Same as Fig. 14 but only for the last cycle of each run. The vertical dashed lines represent the solar forcing sine function phases of −90, 0 and 90^∘.

Figure 16 Ice thickness (unit: m) distribution in four phases for the forcing periods of 20 kyr (upper), 50 kyr (middle) and 100 kyr (lower) from the last cycle of the shortwave oscillation experiment. The corresponding −180, −90, 0 and 90^∘ phase of the solar forcing phases are marked in the headings.

Figure 17 Anomalies of surface temperature (upper; unit: ^∘C), precipitation (middle; unit:  mm d⁻¹) and glacier mask change (lower; brown, cyan and blue represent from ocean to land, from ocean to ice shelf and from land to ocean, respectively) in four phases during the last cycle of the 20 kyr shortwave oscillation experiment. The equilibrium state from coupled dynamic equilibrium experiment is removed to obtain anomalies. The corresponding −180, −90, 0 and 90^∘ phase of the solar forcing phases are marked in the headings.

5.2 Shortwave radiation oscillation experiment

In the following experiments we use the same setup as in the previous section but allow the Northern Hemisphere shortwave radiation, sw, to oscillate, taking the following form:

$$\begin{array}{}\text{(42)}& \text{sw}\left(t\right)=\left(\mathrm{1}+A_{\text{sw}}\cdot \sin \left(\mathrm{2}\mathit{\pi }{\displaystyle \frac{t}{\text{pd}}}\right)\right)\cdot {\text{sw}}_{\text{present}},\end{array}$$

where A_sw is the amplitude of the sw oscillations, which increases from 0 at 13^∘ N to 0.1 at 35^∘ N and maintains 0.1 northward of 35^∘ N. The oscillation period, pd, is set to 20, 50 and 100 kyr in three individual simulations. The sw oscillation is relative to the present-day solar radiation, sw_present. The shortwave maximum amplitude is about 20 W m⁻² at 65^∘ N in the annual mean (Fig. 14a–c) and varies with latitudes and seasons (not shown). The 20, 50 and 100 kyr oscillation periods are simulated for 210, 325 and 350 kyr. The time series for selected climate variables are shown in Fig. 14. The results are shown in reference to the final year of the control run, which is the coupled dynamical equilibrium simulation in Sect. 5.1. To illustrate ice form and retreat in one cycle, we show results from the last forcing cycle of each simulation in Figs. 15–17.

Starting with the 20 kyr oscillation run, there are a number of interesting aspects to point out (Figs. 15a, d, 16a, d and 17a–d). First, at the initial half-cycle, the ice volume is slightly lower than the reference state, indicating a warming period leads to deglaciation (Fig. 14a–c). Then, after the second cycle, the ice volume is always larger than in the control simulation and the cycles are very similar to each other. If we focus on the last cycle of the simulations (Figs. 15–17), we note that T_surf and precipitation are mostly in phase with each other and with the shortwave radiation forcing. The northern hemispheric T_surf oscillation amplitude is about ± 6 ^∘C, and the mean value is clearly below zero (the control run value). This is despite the fact that the mean shortwave radiation is the same as in the control run. This suggests that the oscillating shortwave radiation has a mean cooling effect. This overall cooling is related to the overall increase in the mean ice sheet volume and extent.

It is beyond this study to fully explore how this effect arises, but it is likely to be related to the ice–albedo effect. In the control run the northern hemispheric summer mean ice cover is nearly zero and, with increasing SW forcing, does not decrease much further. However, it can increase substantially for decreased SW forcing, leading to a mean ice cover in the oscillation run that is much larger than in the control. Subsequently, the northern hemispheric albedo is also much higher than in the control leading to a cooler northern hemispheric T_surf.

The ice sheet response to the 20 kyr shortwave oscillation has a number of interesting aspects. As mentioned above, the mean ice sheet volume is larger than in the control run. Indeed, it is never smaller than in the control run, not even at the minimum (compare Figs. 13a and 16a), with the exception of the first cycle. Ice-covered regions and ice volume are out of phase. The ice-covered regions (including land snow and sea ice) grow first and are nearly 180^∘ out of phase with the SW forcing. The ice sheet volume lags behind the ice-covered area and reaches its maximum nearly 90^∘ (a quarter cycle) after the minimum in shortwave radiation (Fig. 15a). This illustrates that the ice sheets have not had enough time to equilibrate with the sw forcing. Further, we can notice that the ice sheet growth and decay is asymmetric, with a slower buildup and faster decay in ice volume, with the reverse pattern in ice sheet area. In the buildup phase the ice sheet extends over large regions at lower latitudes but has relatively thin ice (Fig. 16b). In the decaying phase the ice sheets retreat to higher latitudes and the ice sheet is relatively thick (Fig. 16d).

The northern hemispheric sw forcing also leads to a response in the Southern Hemisphere climate (Fig. 15d). This is mainly due to the GREB-ISM atmospheric heat and moisture transport. It is also partly due to the change in global sea level induced by the northern hemispheric ice sheet changes. The Southern Hemisphere ice sheet changes are in-phase with the Northern Hemisphere climate. It is further noted that the amplitude of the Southern Hemisphere precipitation response relative to T_surf is bigger than in the Northern Hemisphere (compare Fig. 15a and d; given the same scaling factors). This suggests that the moisture transport is more affected by the northern hemispheric climate change than the heat transport.

The longer 50 and 100 kyr period runs show a number of changes relative to the 20 kyr run. First, the ice sheet volume amplitudes increase relative to the 20 kyr run, illustrating that the ice sheets are more sensitive to longer time period forcings (Fig. 15a–c). Second, we see a shift of the maximum ice volume closer to the phase of the minimum of the sw forcing, suggesting that the ice sheets become closer to equilibrium with longer period sw forcing. However, even the 100 kyr oscillation run still shows a significant delay in the ice sheet volume extrema relative to the forcing extrema, indicating that the ice sheets are not yet in equilibrium with the forcings. This illustrates that the intrinsic timescales of the northern hemispheric ice sheets are longer than 100 kyr. It is further interesting to note that the ice sheets can extend over shallow oceanic regions, like the Hudson Bay, Bering Strait or Arctic Sea in the Siberian sector (Fig. 16g and k), but at the same time do not extend into deep ocean regions (compare Fig. 1c with Fig. 16g and k).

The increase in ice thickness response for the longer 50 and 100 kyr period runs has, however, little impact on the amplitudes of the T_surf, precipitation and ice cover response in the Northern Hemisphere, which also occurs in the Southern Hemisphere (Fig. 15e and f). For ice sheets in the Southern Hemisphere, the ice thickness is almost keeping constant, which indicates the Antarctica Ice Sheet in the GREB-ISM is not very sensitive to the orbital forcing in Northern Hemisphere.

6 Summary and discussion

In this study we introduced a newly developed global ice sheet model coupled to the GREB model, defining the new model GREB-ISM. The ice sheet is simulated on the global grid fully interacting with the climate simulation on all grid points. The ice sheet mass balance is driven by accumulation of snow, melting by surface heat fluxes and changes due to ice transport. The ice transport follows the shallow ice approximation for grounded ice and shallow shelf approximation for ice shelves. Sea-ice–climate interactions are also included. The GREB-ISM climate simulation interacts with ice sheets through surface temperature, precipitation, albedo, land–sea mask, topography and sea level. To allow for these interactions, the original GREB model was changed by improving the precipitation simulation of land, including a prognostic sea ice thickness scheme; coupling the surface albedo to the ice thickness; allowing variable land topography as function of ice thickness; introducing global sea level variation and associated changes in land–sea masks; and improving the meridional turbulent, atmospheric heat transport. Thus, the new GREB-ISM is a fully coupled atmosphere, ocean, land and ice sheet model.

We evaluated the performance of the stand-alone ice sheet model in a series of idealized and realistic ice sheet model simulations. We conducted simulations following the EISMINT I and II idealized experiments and found that the GREB-ISM ice sheet model performs similarly to other models with some limitations in the simulation of internal ice temperature. In simulations with realistic climate forcing close to the present day, we found that the equilibrium in Greenland and most of the East Antarctic ice thickness distribution is very similar to that observed, but the West Antarctic Ice Sheet gains too much ice. The overall surface ice velocities and associated calving rates of this model are similar to those observed for both Greenland and East Antarctica.

We investigated the West Antarctic Ice Sheet thickness bias, by evaluating whether uncertainties in precipitation and the parameterization of the ice shelf dynamics (basal melting and viscosity) could cause this bias. However, we found that this bias is unlikely to be caused by these limitations alone, and it is likely to also result from other, so far unknown, limitations in the GREB-ISM model. A possible explanation could be the complexity of the topography and land–sea distribution of West Antarctica and Antarctic Peninsula, which is not well resolved in the current model resolution. So, the coarse grid resolution of this model is likely to play a role in this limitation (Cuzzone et al., 2019).

A time-dependent simulation with simplified surface temperature and precipitation forcing of the past 250 kyr illustrated that the GREB-ISM model can produce a realistic ice sheet response for Greenland, North American and Fennoscandian ice sheets, together with sea level variability. The results for the Antarctic Ice Sheet are less conclusive but may be due to the simplified setup of the experiment.

We further conducted a series of coupled GREB-ISM simulations to evaluate the full interaction of all climate elements in the model. The coupled model simulations produce global equilibrium ice sheets and calving rates very similar to those observed for present-day boundary conditions. Much of this success in creating a realistic global ice sheet is related to the fact that the GREB-ISM model works with flux correction of surface temperature and land precipitation. This leads to realistic mass balance estimates for the ice sheets even in a fully interactive coupled simulation.

When forced with idealized, oscillating solar radiation forcing on the Northern Hemisphere with different oscillation periods (20, 50 and 100 kyr) the model responds with growth of large continental ice sheets and clear interactions with the climate system in the Northern and Southern Hemispheres. The simulations illustrated asymmetries in the buildup and decay of large ice sheets in response to periodic forcing, showing that the ice sheets are more sensitive to longer-timescale forcings. These experiments illustrate the potential of this model for exploring such interactions in future studies.

The current version GREB-ISM is a useful tool to explore ice-sheet–climate global interaction in ice-age cycles. First, a globally fully coupled model enables us to explore the interaction between the two hemispheres. Most previous studies only simulate Northern or Southern Hemisphere ice sheets and take the others as prescribed boundary conditions (Ganopolski et al., 2010; Tigchelaar et al., 2019). Second, the model is very cheap and it has a high potential to do a fully coupled transition simulation for glacial cycles and sensitivity test. The previous studies pointed out that the ice sheet in paleoclimate has multiple stable equilibria (Abe-Ouchi et al., 2013). Therefore, the simulation of transitions are necessary for ice-age cycle processes. For instance, only in transition experiments can we see the ice sheet inertia effect that the longer forcing period leads to stronger ice thickness response (e.g., Fig. 15).

In summary, we presented a new model that is suited for the simulations of global-scale climate variability on timescales of 100 kyr and longer. Given the coarse resolution of the model, it may be less suitable for shorter timescale studies. The model is computationally efficient, calculating 100 000 model years of global simulations per day on a desktop computer, allowing the simulation of the whole Quaternary period (2.6 Myr) within 1 month. For simulations of climate and ice sheet variability over the Quaternary period the GREB-ISM model is, as presented here, a good starting point. Further development may include other relevant climate processes, such as the carbon cycle, deep ocean reservoirs, or the ability of the atmosphere and ocean circulation to respond to changes in topography and the climate state, as well as glacial isostatic adjustment. Such further developments are possible within the framework of the GREB-ISM model and will be addressed in future studies.