A hybrid GCM paleo icesheet model , ANICE 2 . 1 – HadCM 3 @ Bristolv 1 . 0 : set up and benchmark experiments

Fully coupled ice-sheet-climate modelling over 10,000 – 100,000-year time scales on high spatial and temporal resolution remains beyond the capability of current computational systems. Hybrid GCM-ice-sheet modelling 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 3D thermodynamic ice-sheet-shelf model calculating all ice on Earth, with pre-calculated output from 10 several steady-state simulations with the HadCM3 general circulation model (GCM), using a so-called matrix method of coupling both components, where 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. As a benchmark experiment to assess the 15 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 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 20 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 millions of years in duration, capturing climate model driven variations in the mass balance of the ice sheet.


Introduction 25
Sea-level rise due to 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 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License.
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 5 simulations of 10,000 -100,000 years to be carried out within a reasonable amount of time, climate models are much more computationally demanding, limiting simulation time to decadal or centennial time-scales.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 in the long-term interactions between the climate and the cryosphere despite these computational 10 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(Marshall et al., , 2002;;Charbit et al., 2002Charbit et al., , 2007;;Zweck and Huybrechts, 2005;Niu et al., 2017).Both types of studies share the 15 shortcoming of having no clear physical cause for the prescribed climatological variations, and no explicit feedback from the cryosphere back onto the prescribed climate.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 parametrised.20 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 where 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 pre-calculated climate states, a time-continuous climate history can be constructed, which can be 25 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 40km 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 . 30 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 icesheets, is applied to the mass balance module of the ANICE ice-sheet model, which simulates the evolution of all four major

Climate model
HadCM3 is a coupled atmosphere-ocean general circulation model (Gordon et al., 2000;Valdes et al., 2017).It has been shown 20 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 palaeoclimate reconstructions such as PMIP2 (Braconnot et al., 2007) and PlioMIP (Haywood and Valdes, 2003;Dolan et al., 2011Dolan et al., , 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 northsouth resolution of about 278 km, whereas east-west resolution varies from about 70 km over northern Greenland (80 ° latitude) 25 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 Last Glacial Maximum (LGM), pre-industrial (PI) and several time slices in between.Orbital parameters representative of the era are used according to Laskar 30 et al. (2004), atmospheric CO2 concentration is prescribed according to the Vostok ice-core record (190 ppmv at 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 Fig. 2 and Fig. 3.
The modelled glacial-interglacial global mean temperature difference is 4.3 °C, 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 5 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.10

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., 2013de Boer et al., , 2014de Boer et al., , 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 15 andvan de Wal (2005, 2008), which uses present-day monthly precipitation values, where 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).
Surface temperature is calculated from present-day monthly values, including a global temperature offset calculated based on 20 some external forcing, and a constant lapse-rate orographic correction.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.
Horizontal resolution is 20 km for Greenland and 40 km for the other three regions.
In their 2014 study, de Boer et al. simulated global ice distribution and sea level variation over the last 400 kyr, forcing ANICE 25 with the LR04 benthic d 18 O record using the inverse routine by de Boer et al. (2013).Their simulated LGM ice-sheets are shown in Fig. 4.They showed that modelled variations in sea-level, global mean temperature and benthic d 18 O are in good agreement with existing independent literature, although ice-sheet location and extent do not agree well with evidence from geomorphology.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 -30 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 strongly parameterized climate forcing and resulting computational efficiency of ANICE allow these transient simulations of multiple glacial cycles to be carried out within 10 -100h on single-core systems, making ensemble simulations feasible.

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 5 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 pre-calculated 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 span 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 10 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: 15 with pCO2,PI = 280 ppmv and pCO2,LGM = 190 ppmv.To determine the position of the model state along the pCO2 dimension of the climate matrix, we use the EPICA ice core record by Luthi 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 20 represent this process in the model by making the ice-sheet coverage dimension of the climate matrix a spatially variable field  012 (, ), calculated by scaling between the local absorbed insolation at present-day and at LGM.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 a, the latter being calculated internally by ANICE: 25 (2) The weighting field is calculated by scaling between the PI and LGM reference fields: 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: with the weights of the respective unsmoothed, smoothed and average values determined experimentally.For all four icesheets, these spatially variable ice weighting fields are combined with the scalar pCO2 weight wCO2 to yield the final weighting parameter wtot: 10 which is used to linearly interpolate between the states in the climate matrix and calculate the reference temperature, precipitation and orography.Precipitation is customarily interpolated logarithmically to accurately reflect relative changes and to prevent the occurrence of negative values: 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.

10
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 LGM.Here, results from the LGM experiment with HadCM3 (Singarayer and Valdes, 2010) indicate mean annual surface temperatures of around 235 K, or -38 °C.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 to the ice-sheet surface (with 15 an ice thickness of up to 5000 m, as indicated by ICE-5G), the resulting value is about 220 K, or -53 °C, about 15 degrees colder than calculated with the GCM, as shown in Fig. 5.A problem occurs during the inception and the subsequent build-up towards 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.20 A solution to this is to slightly adapt the constant lapse-rate approximation.Assuming the GCM-generated temperature field at 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: calculated temperature between the LGM and PI fields over the ice-free area in the region at LGM.For North America, this results in a value of ∆ +,-= −14.9 K.This methodology ensures that when the modelled ice-sheet is identical to the ICE-5G ice-sheet at LGM and the CO2 concentration is at the LGM values (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 the forcing parameters is at the LGM values, but no ice is present in the model, mean annual surface temperatures are uniformly lower 5 than present-day by ∆ +,-.
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, that the applied lapse-10 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.

Precipitation 15
Whereas a continental-sized ice-sheet influences temperature mainly through albedo, the effect 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.The different character of this process calls for a different representation in the model.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 20 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.The interpolation parameter for the "ice" dimension of the climate matrix wice is expressed as: 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, 25 only the total ice volume limitation is applied and the interpolation weight is calculated as: LGM ice-sheets that were used to calculate the corresponding GCM states are, 5 in many places, thinner at LGM than at present-day, even though the total volume of the ice-sheet is larger.This would mean that an increase of 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 adjusted to account for the difference in 10 temperature between the model state Tmod and the reference GCM state TGCM, as shown in Eq. 12, according to a Clausius-Clapeyron relationship.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.

Last glacial cycle benchmark 20
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 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 25 about 12 hours 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.
Shown in Fig. 6 are the results of this experiment in terms of the global mean sea-level contributions of the four separate icesheets over time, as well as the total global mean sea-level, together with the same values from a simulation of the same period which is in better agreement with the reconstruction.The southern margin lies a little too far to the north, especially in the mid-west.The Antarctic ice-sheet now shows a much stronger increase in ice volume around LGM, matching the 16 m of eustatic sea-level contribution postulated by ICE-5G (Peltier, 2004).Most of the growth ice mass increase takes place happens 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 10 studies of Greenland using the ANICE model (de Boer et al., 2013(de Boer et al., , 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 15 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 LGM shown in Fig. 8, the modelled ice-sheet does not extend as far south over 20 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.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 25 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 sealevel equivalent that matches the DATED-1 horizontal extent very well, as shown in Fig. 9.However, since it was the ICE-5G reconstruction that was used as input for the HadCM3 simulation by Singarayer and Valdes (2004), we aim to maintain consistency and reproduce that particular ice-sheet with our model rather than the DATED-1 LGM ice sheet.We will therefore 30 not use this new version as our benchmark.

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 (Luthi 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 where these parameters had been either 5 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. 10.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.10

Benthic oxygen isotope abundance
Included in ANICE is a module that tracks the oxygen isotope abundances of the ocean (d 18 Osw), precipitation and the icesheets.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 15 anomaly from the climate matrix is used to determine a deep-water temperature anomaly (DTdw), and hence a modelled value for benthic d 18 O.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 4,000-year running average and a scaling factor of 0.25.As opposed to the approach by de Boer et al. (2014), where an inverse method was used to match modelled benthic d 18 O to an externally prescribed record, modelled d 18 O can now be independently compared to such a record in order 20 to test the performance of the matrix method.
We compared our modelled benthic oxygen isotope abundance and the relative contributions to this signal by sea-water heavy oxygen enrichment and deep-water temperature change to data by Shakun et al. (2015), who analysed 49 ODP drilling locations where both surface-dwelling planktonic and benthic foraminiferal oxygen isotope abundance data were available, thereby 25 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

Ice core temperature reconstructions
Shown in Fig. 12 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 5 reconstruction by Jouzel et al. (2007), 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 GISP2 and NGRIP isotope-based reconstructions (Alley, 2000;Kindler et al., 2014) regarding Greenland surface temperature anomalies when neglecting the strong negative excursions during 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 10 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, where the model fails to reproduce the observed rapid cooling.

Conclusions
We have presented and evaluated a hybrid ice-sheet-climate model set-up that combines results from pre-calculated GCM 15 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 20 is such that this simulation could be performed within roughly 12 hours 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.Modelled temperature anomalies over Greenland and 25 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 LGM amounted to 97 ± 6 m sea-level equivalent.
During the first 20 kyr of the inception, the model fails to reproduce the rapid drop in temperature and increase in ice volume 30 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) andde 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.5 The effects of a growing ice-sheet on local and regional temperature are accounted for in the model through the resulting changes in albedo.However, several processes that would influence this effect are currently not included.Glacial-interglacial changes in sea ice cover, though smaller than the changes in land ice cover, would have a strong influence on this process.
Since ANICE does not model sea ice, this effect is not accounted for in the model.Changes in land albedo caused by changing 10 vegetation, although of a smaller magnitude, are likewise not accounted for.Including these 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 as online supplementary material at doi: 15 10.5194/gmd-2018-145supplement, or directly at https://zenodo.org/record/1288386.The source code of ANICE2.1, including the new matrix method, is available online at doi: 10.5194/gmd-2018-145code, or directly at https://zenodo.org/record/1299522.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 c.j.berends@uu.nl.20 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 (P.J.Valdes@bristol.ac.uk).

Acknowledgements
The Ministry of Education, Culture and Science (OCW), in the Netherlands, provided financial support for this study via the Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License.continental ice-sheets (North America, Eurasia, Greenland and Antarctica) simultaneously.Difficulties in bridging the differences in model resolution, as well as other inconsistencies between model states, are addressed and solved.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 LGM that agree very well with geomorphology-based reconstructions.5 Previous work with the ANICE ice-sheet model (de Boer et al., 2013, 2014) used an inverse coupling method, where 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 d 18 O record.This approach essentially determines how global climate should have behaved in order to produce the observed d 18 O record -regardless of what, if anything, could have caused the 10 resulting strong, rapid climatic variations.The new model set-up presented here differs from this approach in that we use the climate matrix to determine what global climate (both temperature and precipitation according to a GCM) would have looked like for a certain prescribed pCO2 and modelled ice-sheet configuration, where the latter is calculated by the ice-sheet model, which is forced by this prescribed climate state.This ensures the constructed climate history is in agreement with the observed pCO2 record and the modelled ice-sheet configuration, thereby capturing the major feedback process between global climate 15 and the cryosphere, where any change in ice-sheet configuration has an immediate impact on local climate through changes in albedo and orographic forcing of precipitation, and any change in climate changes the mass balance of the ice-sheet.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License. 0D: (, ) =  89:,EFG (, ) −  89:,+,-(, )  89:,./(, ) −  89:,+,-(, ) , ,"-(, ) =  OFO • ℎ ./(, ) + (1 −  OFO ) • ℎ +,-(, ).(8) 15Being 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. 20 Since the relative changes in ice-sheet size for Greenland and Antarctica are much smaller than those for North America and Eurasia, the changes in absorbed insolation in those regions should have less impact on local climate.This is reflected in the model by giving more weight to the pCO2 parameter: Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License. OFO = 3 •  "#$ +  012 (, ) 4 .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.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 5

)
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License.The first term in Eq. 10 describes the local ice weighting method byPollard (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.The reason that the local ice thickness term is absent in the calculation for Greenland and Antarctica shown in Eq. 11 is that the ICE-5G Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License. of time with the default ANICE model forced with the LR04 benthic d 18 O 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 LGM are shown in Fig. 7.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. 4.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, 5 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License.
Fig. 11.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 d 18 O at 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 30 model set-up, though still in good agreement with the calculated global mean offset of 2 to 3 K at LGM.The new model set-Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License. up fails to reproduce the strong drop in benthic d 18 O during the inception of the glacial cycle, "catching up" only around 75 kyr.
program of the Netherlands Earth System Science Centre (NESSC).B. de Boer is funded by NWO Earth and Life Sciences 25 (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, 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.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-145Manuscript under review for journal Geosci.Model Dev. Discussion started: 9 July 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 3 :
Figure 3: Annual mean 2m temperature for the Northern Hemisphere (A) and Antarctica (B) and total annual precipitation (C and D), calculated with HadCM3 in the LGM experiment (Singarayer and Valdes, 2010).

Figure 4 :
Figure 4: Ice-sheets (white) and shelves (light blue) at LGM over A) the Northern Hemisphere and B) Antarctica, as simulated with 5

Figure 5 :
Figure 5: Mean annual surface temperature at 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).

Figure 6 :
Figure 6: Global 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 LGM (dashed line).

Figure 7 :
Figure 7: Ice-sheets (white) and shelves (light blue) at LGM over A) the Northern Hemisphere and B) Antarctica, as simulated with the new model set-up.10

Figure 8 : 5 Figure 9 :
Figure 8: Comparison of the simulated Eurasian ice-sheet at LGM with the DATED-1 reconstruction (Hughes et al., 2016).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.

Figure 10 :
Figure 10: Modelled sea level contribution over time for all four individual ice-sheets, and the total sum.Dotted lines indicate upper and lower limits for prescribed CO2 forcing (blue), ablation tuning parameter (green) and SIA/SSA flow enhancement factor ratio (purple).

Figure 11 :
Figure 11: 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) d 18 O 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.

Figure A1 : 5
Figure A1: Annual mean 2m temperature for the Northern Hemisphere (A) and Antarctica (B) and total annual precipitation (C and D), resulting from applying the constant lapse-rate temperature change and the Roe precipitation model to the ERA-40 climate fields.5

Figure A2 :
Figure A2: Annual mean 2m temperature for the Northern Hemisphere (A) and Antarctica (B) and total annual precipitation (C and 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.