Articles | Volume 11, issue 11
Development and technical paper
22 Nov 2018
Development and technical paper |  | 22 Nov 2018

Application of HadCM3@Bristolv1.0 simulations of paleoclimate as forcing for an ice-sheet model, ANICE2.1: set-up and benchmark experiments

Constantijn J. Berends, Bas de Boer, and Roderik S. W. van de Wal

Fully coupled ice-sheet–climate modelling over 10 000–100 000-year timescales at high spatial and temporal resolution remains beyond the capability of current computational systems. Forcing an ice-sheet model with precalculated output from a general circulation model (GCM) offers a middle ground, balancing the need to accurately capture both long-term processes, in particular circulation-driven changes in precipitation, and processes requiring a high spatial resolution like ablation. Here, we present and evaluate a model set-up that forces the ANICE 3-D thermodynamic ice-sheet–shelf model calculating the four large continental ice sheets (Antarctica, Greenland, North America, and Eurasia) with precalculated output from two steady-state simulations with the HadCM3 (GCM) using a so-called matrix method of coupling both components, whereby simulations with various levels of pCO2 and ice-sheet configuration are combined to form a time-continuous transient climate forcing consistent with the modelled ice sheets. We address the difficulties in downscaling low-resolution GCM output to the higher-resolution grid of an ice-sheet model and account for differences between GCM and ice-sheet model surface topography ranging from interglacial to glacial conditions. Although the approach presented here can be applied to a matrix with any number of GCM snapshots, we limited our experiments to a matrix of only two snapshots. As a benchmark experiment to assess the validity of this model set-up, we perform a simulation of the entire last glacial cycle from 120 kyr ago to present day. The simulated eustatic sea-level drop at the Last Glacial Maximum (LGM) for the combined Antarctic, Greenland, Eurasian, and North American ice sheets amounts to 100 m, in line with many other studies. The simulated ice sheets at the LGM agree well with the ICE-5G reconstruction and the more recent DATED-1 reconstruction in terms of total volume and geographical location of the ice sheets. Moreover, modelled benthic oxygen isotope abundance and the relative contributions from global ice volume and deep-water temperature agree well with available data, as do surface temperature histories for the Greenland and Antarctic ice sheets. This model strategy can be used to create time-continuous ice-sheet distribution and sea-level reconstructions for geological periods up to several million years in duration, capturing climate-model-driven variations in the mass balance of the ice sheet.

1 Introduction

Sea-level rise due to the large-scale retreat of the Greenland and Antarctic ice sheets poses one of the main long-term risks of climate change (Church et al., 2013). However, accurate projections of the magnitude and rate of retreat are limited by our understanding of the feedback processes between global climate and the cryosphere on centennial to multi-millennial timescales. One way to test the performance of ice-sheet models that are used for these future projections is to apply these models to ice-sheet evolution in the geological past, both during glacial periods with more ice than present day and warmer periods with less ice (e.g. Bamber et al., 2009; Pollard and DeConto, 2009; de Boer et al., 2013; Dutton et al., 2015).

Ideally, such a model set-up would consist of a general circulation model (GCM) fully coupled to an ice-sheet model, exchanging information every model time step. However, whereas the computational load of typical ice-sheet models allows simulations of 10 000–100 000 years to be carried out within a reasonable amount of time, GCMs are much more computationally demanding, limiting simulation time to decadal or centennial timescales. Fully coupled ice-sheet–climate modelling of complete glacial cycles is therefore not feasible with the current state of model infrastructure.

In order to gain insight into the long-term interactions between the climate and the cryosphere despite these computational limitations, different solutions have been proposed in the past. Several studies of past glacial cycles using ice-sheet models (Bintanja et al., 2002; de Boer et al., 2014) apply a present-day climate with a uniform temperature offset based on a “glacial index”, usually from ice-core isotope records, adapting precipitation based on a Clausius–Clapeyron-type relationship. Others have used a similar glacial index to create a linear combination of output of different GCM time-slice simulations (Marshall et al., 2000, 2002; Charbit et al., 2002, 2007; Tarasov and Peltier, 2004; Zweck and Huybrechts, 2005; Niu et al., 2017). Both types of studies share the shortcoming of having no clear physical cause for the prescribed climatological variations and no explicit feedback from the cryosphere back onto the prescribed climate. Stap et al. (2014, 2016) used a zonally averaged energy balance model coupled to a one-dimensional ice-sheet model to simulate the behaviour of global climate and the cryosphere over millions of years, trading regional details for the ability to simulate long-term feedback processes. Others used dynamically coupled ice-sheet models to Earth system models of intermediate complexity (Charbit et al., 2005; Ganopolski et al., 2010). This approach comes closer to the ideal case of an ice-sheet model fully coupled to a GCM, but since EMICs typically have a coarse spatial resolution, processes influencing the surface mass balance variably over the different parts of the ice sheet (e.g. precipitation, ablation) still need to be parameterized. Other studies have asynchronously coupled ice-sheet models to GCMs (Herrington and Poulsen, 2012) or used fully coupled ice-sheet–GCM set-ups with low-resolution GCMs for shorter periods of model time (Gregory et al., 2012), all showing that non-linear and non-local processes, particularly atmospheric stationary waves, surface albedo, and altitude feedback, can significantly affect the behaviour of ice sheets under a changing climate. Although such studies explicitly describe many more physical processes and feedbacks, computational resources quickly become a limiting factor in the length and number of simulations than can be performed. Abe-Ouchi et al. (2013) performed a very detailed decoupling of the effects on climate of changes in pCO2, albedo, surface elevation, and atmospheric circulation based on several GCM snapshots and used these to force an ice-sheet model in a manner similar to both the glacial index method and the method described in this paper, highlighting the importance of the isostatic adjustment of the lithosphere in producing the 100 kyr glacial cycles. By using precalculated GCM output, this approach makes it possible to run many different simulations and investigate the effects of different physical processes.

The “matrix method” of hybrid ice-sheet–climate modelling (Pollard, 2010; Pollard et al., 2013) is based on a collection of steady-state GCM simulations in which different values for one or more parameters such as pCO2, insolation, or global ice coverage are used to construct a so-called “climate matrix”. By varying these parameters continuously over time and interpolating between these precalculated climate states, a time-continuous climate history can be constructed, which can be used to force an ice-sheet model. Pollard et al. (2013) used this method to simulate the evolution of the Antarctic ice sheet during the early Oligocene for 6 million years using a 40 km resolution ice-sheet model forced with output from the GENESIS version 3 GCM. They concluded that the method had some drawbacks, including a crude albedo feedback and inability to smoothly track orographic precipitation, but that it was adequate for studying the large-scale ice-sheet evolution in which they were interested.

An important difference between the glacial index approach and the matrix method is the latter's more explicit description of the feedback of an expanding or retreating ice sheet on local, regional, and global climate. In a glacial index model, the temporal evolution of the prescribed climatology is determined by an external forcing record (typically pCO2, benthic δ18O, or ice-core isotopes). The matrix method combines this external forcing with one or more internally modelled parameters (typically ice volume or extent) to determine the applied climatology, thus allowing changes in ice-sheet configuration to feed back on climate. Although this approach does still not explicitly describe all the feedback processes that can be included in a fully coupled ice-sheet model – AOGCM set-up, such as the influence of a growing ice-sheet dome on atmospheric circulation and stationary waves and the influence of freshwater fluxes on ocean circulation – it at least partially captures the feedbacks which are not accounted for in a glacial index model and it does not require much more computational resources.

In this study, we constructed a model set-up with a climate matrix consisting of two simulations with the HadCM3 GCM. The climate that is obtained from this matrix, based on the prescribed atmospheric CO2 concentration and internally modelled ice sheets, is applied to the mass balance module of the ANICE ice-sheet model, which simulates the evolution of all four major continental ice sheets (North America, Eurasia, Greenland, and Antarctica) simultaneously. Difficulties in bridging the differences in model resolution and differences in ice-sheet configuration between GCM and ice-sheet model state, especially regarding the orographic forcing of precipitation resulting from ice-sheet advance, are addressed and overcome. As a benchmark experiment, a simulation of the entire last glacial cycle, from 120 kyr to present day, was performed with this model set-up. We show that, because of several improvements to the way changes in albedo and precipitation are handled by the model, we simulate ice sheets at the Last Glacial Maximum (LGM) that agree very well with geomorphology-based reconstructions for Eurasia and better than previous ANICE versions for North America.

Previous work with the ANICE ice-sheet model (de Boer et al., 2013, 2014) used an inverse coupling method, whereby a global temperature offset is calculated in every model time step such that the resulting deep-water temperature, combined with simulated global ice volume, matches a prescribed δ18O record. This approach essentially determines how global climate should have behaved in order to produce the observed δ18O record – regardless of what, if anything, could have caused the resulting strong, rapid climatic variations. Instead of working back from the a posteriori result of benthic δ18O, the new approach presented here starts with the a priori forcings of insolation and pCO2 and determines what global climate should have looked like based on the forcings and the modelled ice sheets. Although this still does not solve the discrepancy between the rapid cooling and sea-level drop suggested by the δ18O record and sea-level data, on the one hand, and the much more gradual decline in pCO2 and surface temperature shown in the ice cores on the other that was observed by other studies (Bintanja and van de Wal, 2008; van de Wal et al., 2011; de Boer et al., 2014; Niu et al., 2017), it might provide new insights on the cause of this discrepancy.

Figure 1LGM ice thickness distributions from the ICE-5G reconstruction (Peltier, 2004) for (a) the Northern Hemisphere and (b) Antarctica. Contour lines for the Northern Hemisphere show ice thickness, and contour lines for Antarctica show surface elevation. Bedrock elevation not covered by ice is shown by colours, with present-day shorelines shown in blue.


2 Methodology

2.1 Climate model

HadCM3 is a coupled atmosphere–ocean general circulation model (Gordon et al., 2000; Valdes et al., 2017). It has been shown to be capable of accurately reproducing the heat budget of the present-day climate (Gordon et al., 2000) and has been used for future climate projections in the IPCC AR4 (Solomon et al., 2007) as well as paleoclimate reconstructions such as PMIP2 (Braconnot et al., 2007) and PlioMIP (Haywood and Valdes, 2003; Dolan et al., 2011, 2015; Haywood et al., 2013). The atmosphere module of HadCM3 covers the entire globe with grid cells of 2.5 latitude by 3.75 longitude, giving a north–south resolution of about 278 km, whereas east–west resolution varies from about 70 km over northern Greenland (80 latitude) to about 290 km over southern Canada (45 latitude, the southern-most area covered by the ANICE grid). The ocean is modelled at a horizontal resolution of 1.25 by 1.25, with 20 vertical layers.

In their 2010 study, Singarayer and Valdes used HadCM3 to simulate global climate during the LGM, the pre-industrial period (PI), and several time slices in between. Orbital parameters representative of the era are used according to Laskar et al. (2004), atmospheric CO2 concentration is prescribed according to the Vostok ice-core record (190 ppmv at the LGM; Petit et al., 1999; Loulergue et al., 2008), and orographic forcing follows the ICE-5G ice distribution reconstruction by Peltier (2004), shown in Fig. 1. Temperature and precipitation fields resulting from these two experiments are shown in Figs. 2 and 3.

Figure 2Annual mean 2 m temperature for the Northern Hemisphere (a) and Antarctica (b) and the total annual precipitation (c, d) calculated with HadCM3 in the PI_Control experiment (Singarayer and Valdes, 2010).


Figure 3Annual mean 2 m temperature for the Northern Hemisphere (a) and Antarctica (b) and the total annual precipitation (c, d) calculated with HadCM3 in the LGM experiment (Singarayer and Valdes, 2010).


The modelled glacial–interglacial global mean temperature difference is 4.3 K, which is in good agreement with results from other model studies (Hewitt et al., 2001; Braconnot et al., 2007), as well as reconstructions from multiple proxies (Jansen et al., 2007; Annan and Hargreaves, 2013). Comparisons of the model results with ice-core isotope temperature reconstructions from Greenland (GRIP; Masson-Delmotte et al., 2005) and Antarctica (EPICA Dome C; Jouzel et al., 2007), as well as borehole-derived surface temperature reconstructions (Dahl-Jensen et al., 1998), indicate that glacial–interglacial temperature changes at these high latitudes are slightly underestimated by the model, by up to 1.5 K over Antarctica and up to 4 K over Greenland.

2.2 Ice-sheet model

To simulate the ice evolution on Earth we use ANICE, a coupled 3-D ice-sheet–shelf model (Bintanja and Van de Wal, 2008; de Boer et al., 2013, 2014, 2015). It combines the shallow ice approximation (SIA) for grounded ice with the shallow shelf approximation (SSA) for floating ice shelves to solve the mechanical equations and incorporates a thermodynamical module to calculate internal ice temperatures. In ANICE, the applied mass balance is calculated using the parameterization by (Bintanja et al., 2005; Bintanja and van de Wal, 2008), which uses present-day monthly precipitation values for which changes in precipitation follow from a Clausius–Clapeyron relation as a function of free atmospheric temperature. Time- and latitude-dependent insolation values according to the reconstruction by Laskar et al. (2004) are used to prescribe incoming radiation at the top of the atmosphere. Ablation is calculated using the surface temperature–albedo–insolation parameterization by Bintanja et al. (2002). In the transition zone near the grounding line, SIA and SSA ice velocities are averaged using the approach by Winkelmann et al. (2011), as explained by de Boer et al. (2013). Sub-shelf melt is calculated based on a combination of the temperature-based formulation by Martin et al. (2011) and the glacial–interglacial parameterization by Pollard and DeConto (2009) tuned by de Boer et al. (2013) to produce realistic present-day Antarctic shelves and grounding lines. A more detailed explanation is provided by de Boer et al. (2013) and references therein. Ice calving is treated by a simple threshold thickness of 200 m, whereby any shelf ice below this thickness is removed. ANICE calculates ice-sheet evolution on four separate grids simultaneously, covering the areas of the large Pleistocene ice sheets: North America, Eurasia, Greenland, and Antarctica. The areas covered by the four model domains are shown in Fig. 4. Horizontal resolution is 20 km for Greenland and 40 km for the other three regions. Splitting North America and Greenland into separate model domains means the Laurentide and Greenland ice sheets can no longer merge in the north, which they might have done during the LGM. However, we assume this to be not important for the large-scale evolution discussed in this study.

In their 2013 study, de Boer et al. simulated global ice distribution and sea-level variation over the last 1 million years, forcing ANICE with the LR04 benthic δ18O record using an inverse routine. Their simulated LGM ice sheets are shown in Fig. 5. They showed that their results are in good agreement with existing independent literature in terms of sea-level contributions (Rohling et al., 2009; Thompson and Goldstein, 2006), seawater heavy isotope enrichment (Duplessy et al., 2002; Lhomme and Clarke, 2005), and other modelling studies (Huybrechts, 2002; Bintanja et al., 2005; Bintanja and van de Wal, 2008; Pollard and DeConto, 2009), although ice-sheet location and extent do not agree well with evidence from geomorphology (Ehlers and Gibbard, 2007; de Boer et al., 2013, and references therein). The latter is likely a result from the absence of feedback from the growth of large ice sheets onto large-scale atmospheric circulation patterns in the model, e.g. failing to reproduce the decrease in precipitation over the Barents Sea–Kara Sea area caused by the appearance of the large Fennoscandian ice dome, resulting in the appearance of an unrealistically large ice dome over the Barents Sea. The highly parameterized climate forcing and resulting computational efficiency of ANICE allow these transient simulations of multiple glacial cycles to be carried out within 10–100 h on single-core systems, making ensemble simulations feasible.

Figure 4The areas of the world covered by the four model domains of ANICE2.1. In the North America and Eurasia domains, Greenland is omitted.


2.3 Climate matrix forcing

A climate matrix, as defined by Pollard (2010), is a collection of output data from different steady-state GCM simulations that differ from each other in one or more key parameters or boundary conditions, such as prescribed atmospheric pCO2, orbital configuration, or ice-sheet configuration. At every point in time during the simulation, the location of the model state within this matrix is extracted from the matrix by interpolating between its constituent precalculated climate states. The pair of climate states generated by Singarayer and Valdes (2010) using HadCM3 is based on otherwise identical input parameters that differ in two respects: pCO2 and ice-sheet coverage. These climate states can be viewed as points on a two-dimensional climate matrix, with pCO2 constituting one dimension and ice-sheet coverage constituting another. In order to calculate a climate state for intermediate pCO2 and ice-sheet coverage values, simple weight functions yielding linear interpolation in this climate phase space will yield the corresponding monthly temperature and precipitation fields. The weighting factor wCO2 is calculated as

(1) w CO 2 = p CO 2 - p CO 2 , LGM p CO 2 , PI - p CO 2 , LGM ,

with pCO2,PI=280 ppmv and pCO2,LGM=190 ppmv. Although the dependence of radiative forcing on pCO2 is logarithmic rather than linear, preliminary experiments showed that changing this in the calculation of the weighting factor did not result in significant changes in modelled sea level at the LGM considering the uncertainty from other model parameters.

To determine the position of the model state along the pCO2 dimension of the climate matrix, we use the EPICA ice-core record by Lüthi et al. (2008). However, the ice-sheet coverage dimension of the matrix, described by wice, is more complicated and cannot be adequately described by a single scalar weight function. Since a continental-sized ice sheet affects both local and global temperature mainly because of the increase in albedo, we chose to represent this process in the model by making the ice-sheet coverage dimension of the climate matrix a spatially variable field wice(x,y) calculated by scaling between the local absorbed insolation at present day and at the LGM. In this way the albedo feedback is captured more realistically. The absorbed insolation Iabs is calculated by multiplying incoming insolation at the top of the atmosphere QTOA (from Laskar et al., 2004) with the surface albedo α, the latter being calculated internally by ANICE:

(2) I abs x , y = 1 - α x , y Q TOA x , y .

The weighting field is calculated by scaling between the PI and LGM reference fields,

(3) w ins x , y = I abs , mod x , y - I abs , LGM x , y I abs , PI x , y - I abs , LGM x , y ,

running from 0 at the LGM to 1 for the PI. To account for both local and regional effects, a Gaussian smoothing filter F with a radius of 200 km and a total average value are added to the weighting field:

(4) w ice x , y = 1 7 w ins x , y + 3 7 F w ins x , y + 3 7 w ins ,

with the weights of the respective unsmoothed, smoothed, and average values determined experimentally such that the precipitation on the ice-sheet flanks, resulting from applying the Roe precipitation model, has values similar to those on the flanks of the ice sheets in the reference GCM snapshots. The value of 200 km for the smoothing radius is based on de Boer et al. (2014), who used a similar smoothing procedure in their precipitation model. Preliminary experiments showed that changing this value did not result in significant changes in modelled LGM sea level, within the uncertainty arising from other model parameters. For all four ice sheets, these spatially variable ice-weighting fields are combined with the scalar pCO2 weight wCO2 to yield the final weighting parameter wtot:

(5) w tot = w CO 2 + w ice 2 ,

which is used to linearly interpolate between the states in the climate matrix and calculate the reference temperature, precipitation, and orography. Preliminary experiments showed that changing the distribution of contributions from wCO2 and wice did not result in significant changes in modelled LGM sea level within the uncertainty arising from other model parameters. Since the two variables generally show coherent temporal behaviour, the two weighting factors are usually close together, meaning wtot does not change much when altering the distribution. When too much weight is given to wice (between 2 and 4 times more than wCO2), eventually a threshold is reached at which the drop in pCO2 during the early phase of the glacial cycle does not result in a strong enough cooling to trigger the growth of ice, thus breaking down this similarity.

Figure 5Ice sheets (white) and shelves (light blue) at the LGM over (a) the Northern Hemisphere and (b) Antarctica, as simulated with the default ANICE version from de Boer et al. (2014). Contour lines for the Northern Hemisphere show ice thickness, and contour lines for Antarctica show surface elevation. Bedrock elevation not covered by ice is shown by colours, with present-day shorelines shown in blue and the ICE-5G ice margin shown in red.


Precipitation is customarily interpolated logarithmically to accurately reflect relative changes and to prevent the occurrence of negative values.


Being linear combinations of output data from a relatively low-resolution GCM, these three data fields necessarily have a lower resolution than the ice-sheet model to which they will be applied. To correct for this, the temperature and precipitation are adapted based on the difference between the interpolated reference orography href,GCM and the actual model orography using the approach by de Boer et al. (2013) described in Appendix A.

Since the relative changes in ice-sheet size for Greenland and Antarctica are much smaller than those for North America and Eurasia, the relative changes in absorbed insolation in those regions are proportionally smaller and should therefore have had less impact on local climate. For example, for North America the total absorbed insolation over the model grid at the LGM is 32 % lower than at present day, whereas for Antarctica this change is only 5 %. This is reflected in the model by giving more weight to the pCO2 parameter:

(9) GRL , ANT : w tot = 3 w CO 2 + w ice x , y 4 .

Preliminary experiments showed that here, too, the sensitivity of the modelled ice volume to this distribution is relatively low.

2.4 Lapse rate

One of the major simplifications in the ANICE mass balance model is the assumption that temperature decreases linearly with altitude – the spatially and temporally constant lapse rate of −8 K km−1. As has been shown by de Boer et al. (2014), the methodology of combining this constant lapse rate with a global temperature offset derived from external forcing produced realistic results in terms of global and regional ice volume when simulating Pleistocene glacial cycles. However, even though the reference orography field obtained from the climate matrix is already close to the model orography and the correction applied to the GCM reference temperature field is therefore much smaller, preliminary experiments showed that even making these relatively small corrections using a constant lapse rate resulted in distorted results.

The limitations of this constant lapse rate procedure can be seen over the western part of Canada, an area that is hypothesized to have remained ice-free for the larger part of the last glacial cycle until a few thousand years before the LGM. Here, results from the LGM experiment with HadCM3 (Singarayer and Valdes, 2010) indicate mean annual surface temperatures of around 235 K, or −38C. When calculating this surface temperature following the approach by de Boer et al. (2014), starting with the present-day surface temperature at bedrock and scaling with the constant lapse rate of −8 K km−1 to the ice-sheet surface (with an ice thickness of up to 5000 m, as indicated by ICE-5G), the resulting value is about 220 K, or −53C, which is about 15C colder than calculated with the GCM, as shown in Fig. 6. A problem occurs during the inception and the subsequent build-up towards the LGM, when this area is still ice-free in the model. Using the GCM-generated temperature field as a reference and scaling this down to bedrock level will then result in surface temperatures that are actually warmer than present day. This is unlikely and results in overestimated melt rates near the ice margins.

Figure 6Mean annual surface temperature at the LGM over North America as generated with HadCM3 by Singarayer and Valdes (2010) (a) versus the temperature field generated for these conditions using a constant lapse rate approach (b). GCM temperatures are substantially higher over the main dome of the ice sheet (area indicated by black circle).


A solution to this is to slightly adapt the constant lapse rate approximation. Assuming the GCM-generated temperature field at the LGM is still based upon the present-day temperature field plus a global offset and a (local) lapse rate correction, similar to the old ANICE method, this local lapse rate correction field is then calculated as

(10) λ LGM x , y = - T LGM x , y - T PI x , y + Δ T LGM h LGM x , y - h PI x , y ,

and the downscaling from the GCM grid to the ice model grid, previously described by Eq. (A1), now being calculated as

(11) T x , y = T ref x , y + λ LGM x , y h x , y - h ref x , y ,

where the local lapse rate at the LGM, λLGM, is calculated by dividing the difference between the local GCM-calculated surface temperature, TLGM, and the extrapolated temperature at local bedrock altitude, Tbedx,y,t=TPIx,y,t+ΔTLGM, by the change in local orography, hLGM, with respect to present day (hPI). The temperature offset ΔTLGM is the mean difference in GCM-calculated temperature between the LGM and PI fields over the ice-free area in the respective model region (either North America or Eurasia) at the LGM. For North America, this results in a value of ΔTLGM=-14.9 K. This methodology ensures that when the modelled ice sheet is identical to the ICE-5G ice sheet at the LGM and the CO2 concentration is at the LGM value (190 ppmv pCO2), the temperature field that is used to calculate the mass balance is still identical to the GCM-calculated temperature field. It also guarantees that, when pCO2 is at 190 ppmv but no ice is present in the model, mean annual surface temperatures are uniformly lower than present day by ΔTLGM.

Of course, the latter scenario only occurs during non-physical steady-state experiments such as forcing ANICE with the LGM GCM climate but initializing with present-day ice cover. During transient experiments, the modelled ice sheets generally resemble those “expected” by the mass balance model through the climate state on which it is based, so the applied lapse rate correction is generally small. This variable lapse rate solution is used in the surface mass balance models for North America and Eurasia, since those regions see the dramatic changes in orography that require this correction. For Greenland and Antarctica, where the changes in ice cover are relatively small even during glacial cycles, the constant lapse rate is still applied with a value of 8 K km−1 based on earlier work with ANICE by Helsen et al. (2013) and de Boer et al. (2014).

2.5 Precipitation

Present-day observations from Greenland indicate that the effect a continental-sized ice sheet has on local precipitation is mostly due to geometry; more precipitation falls on the flanks due to orographic forcing, and as a result the dome becomes a plateau desert (Roe and Lindzen, 2001; Roe, 2002). The different character of this process calls for a different representation in the model than the absorbed insolation-based temperature calculation. In order to calculate monthly precipitation values, for North America and Eurasia we use the “local ice-weighting” method described by Pollard (2010). For every element of the spatial grid, ice thickness relative to the ice thicknesses at that element for the different reference GCM states, limited by the total volume of the ice sheet, is used to obtain the interpolation parameter for the ice dimension of the climate matrix. Although physically, precipitation is influenced by surface altitude and not ice thickness, the fact that the weight is calculated based on scaling the model state between two extremes means the end result is the same as long as the rate of change of ice thickness and surface altitude is the same. The discrepancy between the two is caused by isostatic adjustment. During the inception phase of the glacial cycle, the ice grows slowly enough that there is hardly any discernible time lag between ice thickness and surface altitude. During the deglaciation this is not true anymore, but since ice-sheet evolution during that phase is dominated by ablation rather than precipitation, a parameterization based on elevation instead of ice thickness yields similar results. The interpolation parameter for the “ice” dimension of the climate matrix wice is expressed as

(12) w ice x , y = Hi mod x , y - Hi PI x , y Hi LGM x , y - Hi PI x , y V mod - V PI V LGM - V PI ,

where Himod is the modelled local ice thickness and HiPI and HiLGM are the local ice thickness values in the reference fields from the GCM states. Vmod, VPI, and VLGM are the modelled and reference ice-sheet volumes. For Greenland and Antarctica, only the total ice volume limitation is applied and the interpolation weight is calculated as

(13) w ice x , y = V mod - V PI V LGM - V PI .

The first term in Eq. (12) describes the local ice-weighting method by Pollard (2010), whereas the second term describes the total ice volume scaling. Combining these two terms ensures that precipitation prescribed to the model only decreases over areas where the model actually simulates ice and that the drop in precipitation caused by the ice–plateau–desert effect scales appropriately with ice-sheet size. Since the thickness of a growing ice sheet levels off much earlier than its horizontal extent, an ice sheet only a quarter of its LGM extent can already have nearly the same maximum thickness. Scaling precipitation based on local thickness alone will therefore result in the ice plateau becoming too dry too early in the growth phase, limiting further growth. Preliminary experiments showed that including the total ice-sheet volume in the calculation of the weighting factor solved this problem, resulting in a growth rate more in line with expectations from sea-level records.

The reason that the local ice thickness term is absent in the calculation for Greenland and Antarctica shown in Eq. (13) is that the ICE-5G LGM ice sheets that were used to calculate the corresponding GCM states are, in many places, thinner at the LGM than at present day, even though the total volume of the ice sheet is larger. This would mean that an increase in modelled ice thickness would lead to an increase in applied local precipitation, causing unrealistic ice growth. Therefore, in order to prevent such unrealistic scenarios, precipitation is scaled only by the total ice-sheet volume.

For Greenland and Antarctica, the reference GCM precipitation field PGCM,ref is downscaled from the GCM to the ice-sheet model resolution based on the difference in temperature between the model state Tmod and the reference GCM state TGCM, as shown in Eq. (14), according to a Clausius–Clapeyron relationship similar to the approach by de Boer et al. (2014) described in Appendix A. This ensures that smaller-scale topographical features present in the model but not in the lower-resolution GCM have an influence on local precipitation through their effect on local surface temperature.

(14) P mod x , y = P GCM , ref x , y 1.0266 T mod x , y - T GCM x , y

Similarly, for North America and Eurasia, precipitation is adjusted using the Roe (2002) parameterization for the wind-orography-based correction of precipitation as described in Eqs. (A3)–(A6), but now by using the GCM-generated precipitation and orography as reference fields instead of their ERA-40 equivalents. This allows for a better representation of the orographic forcing of precipitation on the migrating ice flanks as these ice sheets advance and retreat, an effect that cannot be captured by interpolating by different GCM snapshots alone.

3 Results

3.1 Last glacial cycle benchmark

As a benchmark experiment, the new model set-up was used to perform a simulation of the last glacial cycle. The climate matrix for this experiment consists solely of the PI_Control and LGM experiments by Singarayer and Valdes (2010). Following the approach by Bintanja et al. (2002), the model was tuned by adjusting the ablation parameter c3 in Eq. (A9) individually for all four ice-sheet regions such that their modelled sea-level contribution at the LGM matched the values postulated by ICE-5G (Peltier, 2004). The resulting c3 values, which are hereafter kept fixed, are shown in Table 1. This 120 kyr simulation took about 12 h to complete on a single-processor system, meaning it is feasible to use this model set-up to perform ensemble simulations without demanding excessive amounts of computation time.

Table 1Tuned values of the ablation parameter c3 as used in Eq. (A9).

Download Print Version | Download XLSX

Figure 7Global mean sea-level contributions over time for the four individual ice sheets, as well as the global total, for the LGC benchmark experiment (green) and the default ANICE control run (red) compared to the ICE-5G sea level at the LGM for the four individual ice sheets and throughout the last glacial cycle for the global sum (dashed line).


Figure 8Ice sheets (white) and shelves (light blue) at the LGM over (a) the Northern Hemisphere and (b) Antarctica, as simulated with the new model set-up. Contour lines for the Northern Hemisphere show ice thickness, and contour lines for Antarctica show surface elevation. Bedrock elevation not covered by ice is shown by colours, with present-day shorelines shown in blue and the ICE-5G ice margin shown in red.


Shown in Fig. 7 are the results of this experiment in terms of the global mean sea-level contributions of the four separate ice sheets over time, as well as the total global mean sea level, together with the same values from a simulation of the same period of time with the default ANICE model forced with the LR04 benthic δ18O record using an inverse routine. As can be seen, the new model set-up obtains a close match to the postulated ICE-5G LGM ice volume for all ice sheets except Greenland. The resulting ice sheets at the LGM are shown in Fig. 8. As can be seen, the north-west Canadian corridor is now blocked by ice, which was still open in the default ANICE simulation shown earlier in Fig. 5. Although the main dome of the ice sheets is not as thick as in the ICE-5G reconstruction, it now lies more westward than in the simulation with the default ANICE model, forming a ridge running from midwest Canada to the eastern shores of Hudson Bay, which is in better agreement with the reconstruction. The southern margin lies too far to the north, varying from 400 km near the Atlantic coast to up to 950 km in the midwest. The Antarctic ice sheet now shows a much stronger increase in ice volume around the LGM, matching the 16 m of eustatic sea-level contribution postulated by ICE-5G (Peltier, 2004). Most of the ice mass increase takes place in West Antarctica; as can be seen, both the Ross and Ronne shelves become fully grounded. The Greenland ice sheet does show some minor growth over the glacial cycle, though not as much as postulated. It must be noted that several modelling studies of Greenland using the ANICE model (de Boer et al., 2013, 2014) have had trouble in this regard, mostly because of the difficulty in simulating the ice shelves that might have formed around the continent at the time but are not there now (Bradley et al., 2018).

The simulated Eurasian ice sheet is now in better agreement with the consensus regarding the Fennoscandian dome, as well as with the total ice volume or sea-level contribution. When simulated with the default ANICE version, the main dome of the Eurasian ice sheet forms over the Barents Sea, extending eastward to about 70 E. The new model set-up results in a dome over Fennoscandia and a smaller dome over the Barents Sea. The present-day southern North Sea area, formerly Doggerland, remains ice-free, in agreement with paleo-data (Hughes et al., 2016). Compared to the recent DATED-1 reconstruction of the Eurasian ice sheet (Hughes et al., 2016) at the LGM shown in Fig. 9, the modelled ice sheet does not extend as far south over northern Germany, Poland, and Lithuania. The simulated Atlantic side of the ice margin agrees well with the reconstruction, reaching the edge of the continental shelf everywhere.

Figure 9Comparison of the simulated Eurasian ice sheet at the LGM with the DATED-1 reconstruction (Hughes et al., 2016). Contour lines show ice thickness. The modelled ice sheet has a volume of 17 m sea-level equivalent, in agreement with the 17 m of the ICE-5G reconstruction, whereas the DATED-1 ice sheet is equivalent to 24 m sea level.


Peltier (2004) provides an ice volume of the Eurasian ice sheet of about 17 m sea-level equivalent based on GPS observations of isostatic rebound, whereas Hughes et al. (2016) state a volume of 24 m based on geomorphological evidence of the extent and a logarithmic linear regression between ice-sheet area and volume. By slightly increasing the ablation tuning parameter, thus decreasing ablation and increasing ice volume, we were able to produce a Eurasian ice sheet with a volume of 24 m sea-level equivalent that matches the DATED-1 horizontal extent very well, as shown in Fig. 10. However, we believe that a “chain” of simulations such as this (an ice-sheet reconstruction, forcing a GCM, forcing an ice-sheet model) should aim for consistency first, meaning that the ice sheet produced at the end of the chain should match the one that was used as forcing at the start of the chain. Although there are more recent, more extensive data available for the volume and extent of the Eurasian ice sheet, prescribing to the ice-sheet model a climate that was calculated based on the presence of a different ice sheet would make it much more difficult to determine the cause of any model–data mismatch in the final results. We therefore did not use this new, probably more physically realistic Eurasian ice sheet as our benchmark.

3.2 Sensitivity to forcing and model parameters

In order to estimate the uncertainty in modelled global mean sea level following from the uncertainty in the EPICA pCO2 record, we performed simulations with the forcing record adjusted to its respective upper and lower bounds based on an LGM uncertainty of 10 ppmv (Lüthi et al., 2008). Additionally, we investigated the model sensitivity to the four ablation tuning parameters c3 for the different ice sheets mentioned earlier by performing simulations in which these parameters had been either increased or decreased by 10 % relative to their benchmark value. We also assessed model sensitivity to the SSA and SIA flow enhancement factors, with the upper and lower limits determined by Ma et al. (2010), in order to test the sensitivity to the ice-sheet dynamics. Results from these different sensitivity tests are shown in Fig. 11. The resulting uncertainty in simulated LGM ice volume amounts to about 6 m sea-level equivalent in either direction, about 6 % of the total signal, for both the CO2 and ablation parameter experiments. Sensitivity to the flow enhancement factor ratio is lower at about 4 % of the total signal.

Figure 10Comparison of the larger simulated Eurasian ice sheet at the LGM with the DATED-1 reconstruction (Hughes et al., 2016). Contour lines show ice thickness. The modelled ice sheet has a volume of 24 m sea-level equivalent, in agreement with the DATED-1 ice sheet.


Figure 11Modelled sea level contribution over time for all four individual ice sheets and the total sum. The ±2σ confidence interval is shown for the ensemble of simulations from the sensitivity analysis.


3.3 Benthic oxygen isotope abundance

Included in ANICE is a module that tracks the oxygen isotope abundances of the ocean (δ18Osw), precipitation, and the ice sheets. In the default ANICE version, an inverse routine is used to calculate a global temperature offset using the difference between modelled and observed benthic oxygen isotope abundance, implying that modelled and observed are per definition in agreement. In our new model set-up, the isotopic content of the ice sheets is still tracked, but now the global mean temperature anomaly from the climate matrix is used to determine a deep-water temperature anomaly (ΔTdw) and hence a modelled value for benthic δ18O. This deep-water temperature anomaly is calculated from the modelled mean annual surface temperature anomaly over the ocean following the approach by de Boer et al. (2014) using a 4000-year running average and a scaling factor of 0.25. As opposed to the approach by de Boer et al. (2014), in which an inverse method was used to match modelled benthic δ18O to an externally prescribed record, modelled δ18O can now be independently compared to such a record in order to test the performance of the matrix method.

We compared our modelled benthic oxygen isotope abundance and the relative contributions to this signal by seawater heavy oxygen enrichment and deep-water temperature change to the LR04 benthic oxygen isotope stack (Lisiecki and Raymo, 2005) to data by Shakun et al. (2015), who analysed 49 ODP drilling locations at which both surface-dwelling planktonic and benthic foraminiferal oxygen isotope abundance data were available, thereby allowing them to make a data-based decoupling of the contributions from ice volume and deep-water temperature to the benthic oxygen isotope signal. This model–data comparison is shown in Fig. 12. As can be seen, the results from the LGM benchmark experiment are in good agreement with the data, similar to the default ANICE model. The drop in benthic δ18O at the LGM of about 1.7 ‰ is reproduced comparably well by both the inverse-method-forced model by de Boer et al. (2014) and the new matrix-method-forced model set-up. The contribution from the change in deep-water temperature is slightly smaller in the new model set-up, though still in good agreement with the calculated global mean offset of 2 to 3 K at the LGM. The new model set-up fails to reproduce the strong drop in benthic δ18O during the inception of the glacial cycle, “catching up” at only 75 kyr.

Figure 12(a) Modelled benthic oxygen isotope abundance from the default ANICE model (de Boer et al., 2014) and the LGM benchmark experiment compared to different datasets (LR04, Shakun et al., 2015). (b) δ18O of seawater due to depletion of heavy isotopes. (c) Contribution to benthic oxygen isotope abundance due to changes in deep-water temperature. (d) Derived deep-water temperature anomaly.


3.4 Ice-core temperature reconstructions

Shown in Fig. 13 are the modelled mean annual surface temperature anomalies over the Antarctic and Greenland ice sheets for the simulation with the default ANICE version and for the LGC benchmark experiment compared to the EPICA Dome C reconstruction by Jouzel et al. (2007), a stack of the GISP2 reconstruction by Alley (2000), and the NGRIP reconstruction by Kindler et al. (2014). As can be seen, both model versions agree well with each other and reasonably well with the Greenland isotope-based reconstructions (Alley, 2000; Kindler et al., 2014) regarding Greenland surface temperature anomalies. The Greenland records have been smoothed with a 4 kyr running mean to filter out Dansgaard–Oeschger events, which are not present in our model forcing or climate reference runs and are also not included as feedback mechanisms in our model physics. Regarding Antarctic surface temperature anomalies, the new model set-up agrees particularly well with the EPICA isotope-based reconstruction (Jouzel et al., 2007), showing almost no significant deviations except for the first 20 kyr of the inception, during which the model fails to reproduce the observed rapid cooling.

Figure 13Modelled versus reconstructed temperature anomaly for Antarctica (EPICA Dome C; Jouzel et al., 2007) and Greenland (GISP2; Alley, 2000; NGRIP; Kindler et al., 2014).


4 Conclusions

We have presented and evaluated a hybrid ice-sheet–climate model set-up that combines results from precalculated GCM simulations to force an ice-sheet model. Using the matrix method of GCM-ISM coupling, the impacts upon global climate of changes in atmospheric CO2 concentration and global ice distribution are treated separately to construct a time-continuous climate forcing.

As a benchmark experiment, we used this new model set-up to simulate the entire last glacial cycle. Computational efficiency is such that this simulation could be performed within roughly 12 h on a consumer-grade system. When compared with the default ANICE version by de Boer et al. (2014), the new model set-up performed better in simulating the volumes of the continental ice sheets and their geographical position and comparably well at simulating global mean deep-water temperature and isotopic content. The improved performance in terms of geographical position is likely a result of the improved dynamically driven changes in precipitation as solved by the GCM. Niu et al. (2017) showed that forcing the PISM ice-sheet model with output from several different GCM simulations of the LGM from PMIP3, all of which were prescribed the same initial ice sheets, resulted in a wide range of ice-sheet sizes at the LGM (50 to 150 m SLE). This illustrates that, even though the ice sheet prescribed to the GCM leaves a clear local “fingerprint” in the resulting climate, especially in the simulated temperature, this is by no means a guarantee that forcing an ice-sheet model with that climate will reproduce an ice sheet that resembles the ice sheet in the boundary conditions.

Modelled temperature anomalies over Greenland and Antarctica agree well with ice-core isotope-based reconstructions. When accounting for uncertainty in the applied forcing and model parameters, the simulated volume of the four major continental ice sheets (excluding contributions from smaller ice caps, glaciers, thermal expansion, and ocean area changes) at the LGM amounted to 97±11 m sea-level equivalent (±2σ from the ensemble of simulations from the sensitivity analysis).

During the first 20 kyr of the inception, the model fails to reproduce the rapid drop in temperature and increase in ice volume visible in both benthic oxygen isotope records and ice-core isotope-based temperature reconstructions, implying that pCO2 forcing alone is not sufficient to explain these observations without including some additional non-linear feedback processes. This is in line with results from other studies; studies like van de Wal et al. (2011) and de Boer et al. (2014) were able to reproduce the rapid cooling by using a forcing, such as a benthic oxygen isotope stack, that already incorporated the rapid decrease during the initial phase of the glacial cycle, whereas Bintanja and van de Wal (2008) and Niu et al. (2017) were unable to reproduce the rapid ice growth with pCO2 forcing alone.

The effects of a growing ice sheet on local and regional temperature are accounted for in the model through the resulting changes in albedo, but non-linear and non-local effects remain difficult to capture. Abe-Ouchi et al. (2013) constructed a model set-up similar to the matrix method presented here, but with more dimensions and corresponding GCM snapshots added to the matrix to decouple the different processes affecting temperature more explicitly: pCO2, albedo, altitude, and atmospheric stationary waves. Although their modelled ice sheets at the LGM do not match geomorphological reconstructions or the results presented here, they do report a stronger increase in ice volume during the inception. Expanding our climate matrix along the lines of their approach to more accurately describe the interplay between ice and climate for smaller ice sheets could therefore potentially solve some of the repeatedly observed discrepancy between sea-level records and benthic δ18O records, on the one hand, and pCO2 and temperature records on the other hand.

Other processes not accounted for in the albedo-based parameterization of our climate matrix include glacial–interglacial changes in sea-ice cover and changes in land albedo caused by changing vegetation. Including these feedback processes in the model could improve model performance in terms of the quantitative relation between pCO2 and ice volume.

Code and data availability

NetCDF files containing output data from the benchmark simulation (ice thickness, bedrock topography, mean annual temperature, annual precipitation, albedo, and surface mass balance) are available online in the Supplement at (Berends et al., 2018a).

The source code of ANICE2.1, including the new matrix method, is available online at (Berends et al., 2018b). Note that the model code can be compiled but cannot be run without input data describing present-day climate and topography, initial ice thickness and topography, and GCM output files constituting the climate matrix. For any questions regarding ANICE, please contact

The output of the HadCM3 experiments which we used to construct the climate matrix can be obtained from Paul Valdes at the University of Bristol (

Appendix A: Mass balance

In the ANICE version used by de Boer et al. (2014), the entire mass balance module is forced by a global temperature offset calculated from a prescribed δ18O value and modelled global ice volume using the inverse routine by de Boer et al. (2013). This temperature offset, combined with a constant lapse rate orography correction to account for changing ice thickness, is used to calculate a new monthly surface temperature field in every model time step:

(A1) T x , y = T ref x , y + d T glob + λ h x , y - h ref x , y .

Thus, the applied temperature T at horizontal location x,y is calculated at every model time step from the ERA-40 reference temperature Tref, the global temperature offset dTglob, and the difference between the model orography h and the reference orography href multiplied by the constant lapse rate λ of −8 K km−1. For Greenland and Antarctica, the applied precipitation P is then calculated by correcting the monthly present-day reference value Pref based on the difference between applied and reference temperature (Jouzel and Merlivat, 1984; Huybrechts, 1992):

(A2) P x , y = P ref x , y 1.0266 T x , y - T ref x , y .

When simulating entire glacial cycles, the changes in ice-sheet geometry over North America and Eurasia are of a much larger scale than those over Greenland and Antarctica. In order to recreate the hypothesized westward growth of those ice sheets during glacial inception caused by orographic forcing of precipitation as moist wind blows up the slope of the ice sheet and releases its moisture content, the precipitation model by Roe and Lindzen (2001) and Roe (2002) is used to calculate monthly precipitation values over these regions:


Here, esat is the saturation vapour pressure at the surface, which is a good proxy for the moisture content of the overlying air column. It is described by the Clausius–Clapeyron in Eq. (A5) using the monthly mean surface temperature T, where e0=6.112 mbar, c1=17.67, and c2=243.5 K. The vertical wind velocity wvv is calculated from the 850 hPa wind and the surface gradient according to Eq. (A7). The precipitation PRoe is related to vertical wind velocity wvv through a probability distribution fwvvdwvv, which is the probability that wvv lies between wvv and wvv+dwvv according to Eq. (A6), where N is a normalization factor and α=1.15 cm s−1 is the measure of variability (Roe, 2002) in the vertical wind velocity. The precipitation PRoe is given by Eq. (A4), where the constants a=2.5×10-11 kg−1 s2 m and b=5.9×10-9 s3 kg were obtained by tuning to observations of Greenland (Roe, 2002). Equation (A4) is solved analytically using error functions (Roe and Lindzen, 2001).

Figure A1Annual mean 2 m temperature for the Northern Hemisphere (a) and Antarctica (b) and the total annual precipitation (c, d) resulting from applying the constant lapse rate temperature change and the Roe precipitation model to the ERA-40 climate fields.


Figure A2Annual mean 2 m temperature for the Northern Hemisphere (a) and Antarctica (b) and the total annual precipitation (c, d) resulting from applying the constant lapse rate temperature change plus global offset and the Roe precipitation model to the ERA-40 climate fields and the ANICE LGM ice sheets.


Both wvv and esat are calculated for both the reference state, using the reference temperature and orography fields, and for the model state, using the values at that model time step. The relative difference between the two modelled precipitation fields resulting from Eq. (A4) is applied as an anomaly to the reference precipitation field to yield the applied precipitation field as described by Eq. (A3).

Figures A1 and A2 show the mean annual temperature and total annual precipitation fields at present day and the LGM, respectively, resulting from applying these two methods to the initial ERA-40 temperature and precipitation fields using the difference between the reference ERA-40 orography and the modelled orography at present day and the LGM.

The monthly surface mass balance is calculated from the applied surface temperature and precipitation fields and the prescribed incoming radiation at the top of the atmosphere following Laskar et al. (2004). Monthly values for accumulation, refreezing, and ablation are calculated separately and added. First, the snow fraction of precipitation is calculated according to the parameterization by Ohmura (1999):

(A8) f snow x , y = 1 - 0.796 tan - 1 T x , y - T 0 3.5 2 ,

where the spatially variable monthly snow fraction fsnow is defined as a function of 2 m air temperature. Monthly accumulation is simply the product of this fraction and monthly precipitation:

(A9) Acc x , y = P x , y f snow x , y .

Local monthly ablation Abl is parameterized as a function of the 2 m air temperature Tano, albedo α, and incoming solar radiation at the top of the atmosphere QTOA following the approach by Bintanja et al. (2002):


with c1=0.0788 m yr−1 K−1, c2=0.004 m3 J−1, and c3 is a tuning parameter different for each individual ice sheet (tuned values listed in Table 1).

The local monthly refreezing Refr is calculated from the available liquid water content Lw (the sum of liquid precipitation and ablation) and the superimposed water content Lsup following the approach by Huybrechts and de Wolde (1999) and Janssens and Huybrechts (2000):


The surface mass balance SMB that will be used by the ice-sheet model is calculated as the sum of the accumulation Acc, the refreezing Refr, and the ablation Abl:

(A14) SMB x , y = Acc x , y + Refr x , y - Abl x , y .
Author contributions

CJB, BdB, and RSWvdW designed the study. CJB created the model set-up and carried out the simulations, with support from BdB and RSWvdW. CJB drafted the paper, and all authors contributed to the final version.

Competing interests

The authors declare that they have no conflict of interest.


The Ministry of Education, Culture and Science (OCW) in the Netherlands provided financial support for this study via the programme of the Netherlands Earth System Science Centre (NESSC). Bas de Boer is funded by NWO Earth and Life Sciences (ALW), project 863.15.019. This work was sponsored by NWO Exact and Natural Sciences for the use of supercomputer facilities. Model runs were performed on the LISA Computer Cluster, and we would like to acknowledge SurfSARA Computing and Networking Services for their support. Special thanks go to Paul Valdes for sharing the data from his HadCM3 simulations with us.

Edited by: Steven Phipps
Reviewed by: Lev Tarasov and Fuyuki Saito


Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., and Blatter, H.: Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume, Nature Letters 500, 190–194, 2013. 

Alley, R. B.: The Younger Dryas cold interval as viewed from central Greenland, Quaternary Sci. Rev., 19, 213–226, 2000. 

Annan, J. D. and Hargreaves, J. C.: A new global reconstruction of temperature changes at the Last Glacial Maximum, Clim. Past, 9, 367–376,, 2013. 

Bamber, J. L., Riva, R. E. M., Vermeersen, B. L. A., and LeBrocq, A. M.: Reassessment of the potential sea-level rise from a collapse of the West Antarctic Ice Sheet, Science, 324, 901–903, 2009. 

Berends, C., de Boer, B., and van de Wal, R.: Berends_etal_2018_GMD_supplement, [Data set],, 2018a. 

Berends, C., de Boer, B., and van de Wal, R.: Berends_etal_2018_GMD_code,, 2018b. 

Bintanja, R. and van de Wal, R. S. W.: North American ice-sheet dynamics and the onset of 100,000-year glacial cycles, Nature, 454, 869–872, 2008. 

Bintanja, R., van de Wal, R. S. W., and Oerlemans, J.: Global ice volume variations through the last glacial cycle simulated by a 3-D ice-dynamical model, Quatern. Int., 95–96, 11–23, 2002. 

Bintanja, R., van de Wal, R., and Oerlemans, J.: Modelled atmospheric temperatures and global sea levels over the past million years, Nature, 437, 125–128, 2005. 

Braconnot, P., Otto-Bliesner, B., Harrison, S., Joussaume, S., Peterchmitt, J.-Y., Abe-Ouchi, A., Crucifix, M., Driesschaert, E., Fichefet, Th., Hewitt, C. D., Kageyama, M., Kitoh, A., Laîné, A., Loutre, M.-F., Marti, O., Merkel, U., Ramstein, G., Valdes, P., Weber, S. L., Yu, Y., and Zhao, Y.: Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum – Part 1: experiments and large-scale features, Clim. Past, 3, 261–277,, 2007. 

Bradley, S. L., Reerink, T. J., van de Wal, R. S. W., and Helsen, M. M.: Simulation of the Greenland Ice Sheet over two glacial–interglacial cycles: investigating a sub-ice-shelf melt parameterization and relative sea level forcing in an ice-sheet–ice-shelf model, Clim. Past, 14, 619–635,, 2018. 

Charbit, S., Ritz, C., and Ramstein, G.: Simulations of Northern Hemisphere ice-sheet retreat: sensitivity to physical mechanisms involved during the Last Deglaciation, Quaternary Sci. Rev., 21, 243–265, 2002. 

Charbit, S., Kageyama, M., Roche, D., Ritz, C., and Ramstein, G.: Investigating the mechanisms leading to the deglaciation of past continental northern hemisphere ice sheets with the CLIMBER–GREMLINS coupled model, Global Planet. Change, 48, 253–273, 2005. 

Charbit, S., Ritz, C., Philippon, G., Peyaud, V., and Kageyama, M.: Numerical reconstructions of the Northern Hemisphere ice sheets through the last glacial-interglacial cycle, Clim. Past, 3, 15–37,, 2007. 

Church, J. A., Clark, P. U., Cazenave, A., Gregory, J. M., Jevrejeva, S., Levermann, A., Merrifield, M. A., Milne, G. A., Nerem, R. S., Nunn, P. D., Payne, A. J., Pfeffer, W. T., Stammer, D., and Unnikrishnan, A. S.: Sea-Level Change, Climate Change 2013: The Physical Science Basis, Contribution of Working Group 1 to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 2013. 

Dahl-Jensen, D., Mosegaard, K., Gundestrup, N., Clow, G. D., Johnsen, S. J., Hansen, A. W., and Balling, N.: Past Temperatures Directly from the Greenland Ice Sheet, Science, 282, 268–272, 1998. 

de Boer, B., van de Wal, R. S. W., Lourens, L. J., Bintanja, R., and Reerink, T. J.: A continuous simulation of global ice volume over the past 1 million years with 3-D ice-sheet models, Clim. Dynam., 41, 1365–1384, 2013. 

de Boer, B., Stocchi, P., and van de Wal, R. S. W.: A fully coupled 3-D ice-sheet–sea-level model: algorithm and applications, Geosci. Model Dev., 7, 2141–2156,, 2014. 

de Boer, B., Dolan, A. M., Bernales, J., Gasson, E., Goelzer, H., Golledge, N. R., Sutter, J., Huybrechts, P., Lohmann, G., Rogozhina, I., Abe-Ouchi, A., Saito, F., and van de Wal, R. S. W.: Simulating the Antarctic ice sheet in the late-Pliocene warm period: PLISMIP-ANT, an ice-sheet model intercomparison project, The Cryosphere, 9, 881–903,, 2015. 

Dolan, A. M., Haywood, A. M., Lunt, D. J., Dowsett, H. J., Hunter, S. J., Lunt, D. J., and Pickering, S. J.: Sensitivity of Pliocene ice sheets to orbital forcing, Palaeogeogr. Palaeocl., 309, 98–110, 2011. 

Dolan, A. M., Haywood, A. M., Hunter, S. J., Tindall, J. C., Dowsett, H. J., Hill, D. J., and Pickering, S. J.: Modelling the enigmatic Late Pliocene Glacial Event – Marine Isotope Stage M2, Global Planet. Change, 128, 47–60, 2015. 

Duplessy, J.-C., Labeyrie, L., and Waelbroeck, C.: Constraints on the ocean oxygen isotopic enrichment between the Last Glacial Maximum and the Holocene: Paleoceanographic implications, Quaternary Sci. Rev., 21, 315–330, 2002. 

Dutton, A., Carlson, A. E., Long, A. J., Milne, G. A., Clark, P. U., DeConto, R. M., Horton, B. P., Rahmstorf, S., and Raymo, M. E.: Sea-level rise due to polar ice-sheet mass loss during past warm periods, Science, 349, 153–163, 2015. 

Ehlers, J. and Gibbard, P. L.: The extent and chronology of Cenozoic Global Glaciation, Quatern. Int., 164, 6–20, 2007. 

Ganopolski, A., Calov, R., and Claussen, M.: Simulation of the last glacial cycle with a coupled climate ice-sheet model of intermediate complexity, Clim. Past, 6, 229–244,, 2010. 

Gordon, C., Cooper, C., Senior, C. A., Banks, H., Gregory, J. M., Johns, T. C., Mitchell, J. F. B., and Wood, R. A.: The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments, Clim. Dynam., 16, 147–168, 2000. 

Gregory, J. M., Browne, O. J. H., Payne, A. J., Ridley, J. K., and Rutt, I. C.: Modelling large-scale ice-sheet–climate interactions following glacial inception, Clim. Past, 8, 1565–1580,, 2012. 

Haywood, A. M. and Valdes, P. J.: Modelling Pliocene warmth: contribution of atmosphere, oceans and cryosphere, Earth Planet. Sc. Lett., 218, 363–377, 2003. 

Haywood, A. M., Hill, D. J., Dolan, A. M., Otto-Bliesner, B. L., Bragg, F., Chan, W.-L., Chandler, M. A., Contoux, C., Dowsett, H. J., Jost, A., Kamae, Y., Lohmann, G., Lunt, D. J., Abe-Ouchi, A., Pickering, S. J., Ramstein, G., Rosenbloom, N. A., Salzmann, U., Sohl, L., Stepanek, C., Ueda, H., Yan, Q., and Zhang, Z.: Large-scale features of Pliocene climate: results from the Pliocene Model Intercomparison Project, Clim. Past, 9, 191–209,, 2013. 

Helsen, M. M., van de Berg, W. J., van de Wal, R. S. W., van den Broeke, M. R., and Oerlemans, J.: Coupled regional climate–ice-sheet simulation shows limited Greenland ice loss during the Eemian, Clim. Past, 9, 1773–1788,, 2013. 

Hewitt, C. D., Broccoli, A. J., Mitchell, J. F. B., and Stouffer, R. J.: A coupled model study of the last glacial maximum: Was part of the North Atlantic relatively warm?, Geophys. Res. Lett., 28, 1571–1574, 2001. 

Herrington, A. and Poulsen, C. J.: Terminating the Last Interglacial: The Role of Ice Sheet-Climate Feedbacks in a GCM Asynchronously Coupled to an Ice Sheet Model, J. Climate, 25, 1871–1882, 2012. 

Hughes, A. L. C., Gyllencreutz, R., Lohne, Ø., Mangerud, J., and Svendsen, J. I.: The last Eurasian ice sheets – a chronological database and time-slice reconstruction, DATED-1, Boreas, 45, 1-45, 2016. 

Huybrechts, P.: The Antarctic ice sheet and environmental change: a three-dimensional modelling stud,, Berichte zur Polarforschung, 99, available at: (last access: 19 November 2018), 1992. 

Huybrechts, P.: Sea-level changes at the LGM from ice-fynamic reconstructions of the Greenland and Antarctic ice sheets during the glacial cycles, Quaternary Sci. Rev., 21, 203–231, 2002. 

Huybrechts, P. and de Wolde, J.: The Dynamic Response of the Greenland and Antarctic Ice Sheets to Multiple-Century Climatic Warming, J. Climate, 12, 2169–2188, 1999. 

Jansen, E., Overpeck, J., Briffa, K. R., Duplessy, J.-C., Joos, F., Masson-Delmotte, V., Olago, D. O., Otto-Bliesner, B., Peltier, W. R., Rahmstorf, S., Ramesh, R., Raynaud, D., Rind, D. H., Solomina, O., Villalba, R., and Zhang, D.: Palaeoclimate, in: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental panel on Climate Change, edited by: Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Averyt, K. B., Tignor, M., Mand iller, H. L., Cambridge University Press, Cambridge, UK, and New York, 433–497, 2007. 

Janssens, I. and Huybrechts, P.: The treatment of meltwater retention in mass-balance parameterizations of the Greenland ice sheet, Ann. Glaciol., 31, 133–140, 2000. 

Jouzel, J. and Merlivat, L.: Deuterium and Oxygen 18 in Precipitation: Modelling of the Isotopic Effects During Snow Formation, J. Geophys. Res., 89, 11749–11757, 1984. 

Jouzel, J., Masson-Delmote, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J. M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, J., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and Millenioal Antarctic Climate Variability over the Past 800,000 Years, Science, 317, 793–797, 2007. 

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902,, 2014. 

Laskar, J., Robutel, P., Gastineau, M., Correia, A. C. M., and Levrard, B.: A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285, 2004. 

Lhomme, N. and Clarke, G. K. C.: Global budget of water isotopes inferred from polar ice sheets, Geophys. Res. Lett., 32, L20502,, 2005. 

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic delta-18-O records, Paleoceanography, 20, L20502,, 2005. 

Loulergue, L., Schilt, A., Spahni, R., Masson-Delmotte, V., Blunier, T., Lemieux, B., Barnola, J.-M., Raynaud, D., Stocker, T. F., and Chappellaz, J.: Orbital and millennial-scale features of atmospheric CH4 over the past 800,000 years, Nature, 453, 383–386, 2008. 

Lüthi, D., Le Floch, M., Bereiter, B., Blunier, T., Barnola, J.-M., Siegenthaler, U., Raynaud, D., Jouzel, J., Fisher, H., Kawamura, K., and Stocker, T. F.: High-resolution carbon dioxide concentration record 650,000-800,000 years before present, Nature Letters, 452, 379–382, 2008. 

Ma, Y., Gagliardini, O., Ritz, C., Gillet-Chaulet, F., Durand, G., and Montagnat, M.: Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model, J. Glaciol., 56, 805–812, 2010. 

Marshall, S. J., Tarasov, L., Clarke, G. K., and Peltier, W. R.: Glaciological reconstruction of the Laurentide Ice Sheet: physical processes and modelling challenges, Can. J. Earth Sci., 37, 769–793, 2000. 

Marshall, S. J., James, T. S., and Clarke, G. K.: North American ice sheet reconstructions at the Last Glacial Maximum, Quaternary Sci. Rev., 21, 175–192, 2002. 

Martin, M. A., Winkelmann, R., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 2: Dynamic equilibrium simulation of the Antarctic ice sheet, The Cryosphere, 5, 727–740,, 2011. 

Masson-Delmote, V., Jouzel, J., Landais, A., Stievenard, M., Johnsen, S. J., White, J. W. C., Werner, M., Sveinbjornsdottir, A., and Fuhrer, K.: GRIP Deuterium Excess Reveals Rapid and Orbital-Scale Changes in Greenland Moisture Origin, Science, 309, 118–122, 2005. 

Niu, L., Lohmann, G., Hinck, S., and Gowan, E. J.: Sensitivity of atmospheric forcing on Northern Hemisphere ice sheets during the last glacial-interglacial cycle using output from PMIP3, Clim. Past Discuss.,, in review, 2017. 

Ohmura, A.: Precipitation, accumulation and mass balance of the Greenland ice sheet, Zeitschrift fur Gletscherkunde und Glazialgeologie, 35, 1–20, 1999. 

Peltier, W. R.: Global Glacial Isostasy and the Surface of the Ice-Age Earth: The ICE-5G (VM2) Model and GRACE, Annu. Rev. Earth Pl. Sc., 32, 111–149, 2004. 

Petit, J. R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chappellaz, J., Davis, J., Delaygue, G., Delmotte, M., Kotlyakov, V. M., Legrand, M., Lipenkov, V., Lorius, C., Pepin, L., Ritz, C., Saltzman, E., and Stievenard, M.: Climate and atmospheric history of the past 420 000 years from the Vostok Ice Core, Antarctica, Nature, 399, 429–436, 1999. 

Pollard, D.: A retrospective look at coupled ice sheet-climate modeiing, Climatic Change, 100, 173–194, 2010. 

Pollard, D. and DeConto, R. M.: Modelling West Antarctic ice sheet growth and collapse through the past five million years, Nature, 458, 329–332, 2009. 

Pollard, D., Kump, L. R., and Zachos, J. C.: Interactions between carbon dioxide, climate, weathering, and the Antarctic ice sheet in the earliest Oligocene, Global Planet. Change, 111, 258–267, 2013. 

Roe, G. H.: Modelling precipitation over ice sheets: an assessment using Greenland, J. Glaciol., 48, 70–80, 2002. 

Roe, G. H. and Lindzen, R. S.: The Mutual Interaction between Continental-Scale Ice Sheets and Atmospheric Stationary Waves, J. Climate, 14, 1450–1465, 2001. 

Rohling, E. J., Grant, K., Bolshaw, M., Roberts, A. P., Siddall, M., Hemleben, Ch., and Kucera, M.: Antarctic temperature and global sea level closely coupled over the past five glacial cycles, Nat. Geosci. Lett., 2, 500–504, 2009. 

Shakun, J. D., Lea, D. W., Lisiecki, L. E., and Raymo, M. E.: An 800-kyr record of global surface ocean d18O and implications for ice volume-temperature coupling, Earth Planet. Sc. Lett., 426, 58–68, 2015. 

Singarayer, J. S. and Valdes, P. J.: High-latitude climate sensitivty to ice-sheet forcing over the last 120 kyr, Quaternary Sci. Rev., 29, 43–55, 2010.  

Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., and Averyt, K. B.: Climate Change 2007: the physical science basis, Contribution Of Working Group I To The Fourth Assessment Report Of The Intergovernmental Panel On Climate Change, 2007. 

Stap, L. B., van de Wal, R. S. W., de Boer, B., Bintanja, R., and Lourens, L. J.: Interaction of ice sheets and climate during the past 800 000 years, Clim. Past, 10, 2135–2152,, 2014. 

Stap, L. B., de Boer, B., Ziegler, M., Bintanja, R., Lourens, L. J., and van de Wal, R. S. W.: CO2 over the past 5 million years: Continuous simulation and new d11B-based proxy data, Earth Planet. Sc. Lett., 439, 1–10, 2016. 

Tarasov, L. and Peltier, W. R.: A geophysically constrained large ensemble analysis of the deglacial history of the North American ice-sheet complex, Quaternary Sci. Rev., 23, 359–388, 2004. 

Thompson, W. G. and Goldstein, Steven L.: A radiometric calibration of the SPECMAP timescale, Quaternary Sci. Rev., 25, 3207–3215, 2006. 

Uppala, S. M., Kallberg, P. W., Simmons, A. J., Andrae, U., da Costa Bechtold, V., Fiorino, M., Gibson, J. K., Haseler, J., Hernandez, A., Kelly, G. A., Li, X., Onogi, K., Saarinen, S., Sokka, N., Allan, R. P., Andersson, E., Arpe, K., Balmaseda, M. A., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Caires, S., Chevallier, F., Dethof, A., Dragosavac, M., Fisher, M., Fuentes, M., Hagemann, S., Holm, E., Hoskins, B. J., Isaksen, L., Janssen, P. A. E. M., Jenne, R., McNally, A. P., Mahfouf, J.-F., Morcrette, J.-J., Rayner, N. A., Saunders, R. W., Simon, P., Sterl, A., Trenberth, K. E., Untch, A., Vasiljevic, D., Viterbo, P., and Woollen, J.: The ERA-40 re-analysis, Q. J. Roy. Meteor. Soc., 131, 2961–3012, 2005. 

Valdes, P. J., Armstrong, E., Badger, M. P. S., Bradshaw, C. D., Bragg, F., Crucifix, M., Davies-Barnard, T., Day, J. J., Farnsworth, A., Gordon, C., Hopcroft, P. O., Kennedy, A. T., Lord, N. S., Lunt, D. J., Marzocchi, A., Parry, L. M., Pope, V., Roberts, W. H. G., Stone, E. J., Tourte, G. J. L., and Williams, J. H. T.: The BRIDGE HadCM3 family of climate models: HadCM3@Bristol v1.0, Geosci. Model Dev., 10, 3715–3743,, 2017. 

Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 1: Model description, The Cryosphere, 5, 715–726,, 2011. 

van de Wal, R. S. W., de Boer, B., Lourens, L. J., Köhler, P., and Bintanja, R.: Reconstruction of a continuous high-resolution CO2 record over the past 20 million years, Clim. Past, 7, 1459–1469,, 2011. 

Zweck, C. and Huybrechts, P.: Modeling of the northern hemisphere ice sheets during the last glacial cycle and glaciological sensitivity, J. Geophys. Res., 110, D07103,, 2005. 

Short summary
We have devised a novel way to couple a climate model to an ice-sheet model. Usually, climate models are too slow to simulate more than a few centuries, whereas our new model set-up can simulate a full 120 000-year ice age in about 12 h. This makes it possible to look at the interactions between global climate and ice sheets on long timescales, something which is relevant for both research into past climate and future projections.