Articles | Volume 19, issue 2
Research article
08 Feb 2023
Research article |  | 08 Feb 2023

Modelling feedbacks between the Northern Hemisphere ice sheets and climate during the last glacial cycle

Meike D. W. Scherrenberg, Constantijn J. Berends, Lennert B. Stap, and Roderik S. W. van de Wal

During the last glacial cycle (LGC), ice sheets covered large parts of Eurasia and North America, which resulted in ∼120 m of sea level change. Ice sheet–climate interactions have considerable influence on temperature and precipitation patterns and therefore need to be included when simulating this time period. Ideally, ice sheet–climate interactions are simulated by a high-resolution Earth system model. While these models are capable of simulating climates at a certain point in time, such as the pre-industrial (PI) or the Last Glacial Maximum (LGM; 21 000 years ago), a full transient glacial cycle is currently computationally unfeasible as it requires a too-large amount of computation time. Nevertheless, ice sheet models require forcing that captures the gradual change in climate over time to calculate the accumulation and melt of ice and its effect on ice sheet extent and volume changes.

Here we simulate the LGC using an ice sheet model forced by LGM and PI climates. The gradual change in climate is modelled by transiently interpolating between pre-calculated results from a climate model for the LGM and the PI. To assess the influence of ice sheet–climate interactions, we use two different interpolation methods: the climate matrix method, which includes a temperature–albedo and precipitation–topography feedback, and the glacial index method, which does not. To investigate the sensitivity of the results to the prescribed climate forcing, we use the output of several models that are part of the Paleoclimate Modelling Intercomparison Project Phase III (PMIP3). In these simulations, ice volume is prescribed, and the climate is reconstructed with a general circulation model (GCM). Here we test those models by using their climate to drive an ice sheet model over the LGC.

We find that the ice volume differences caused by the climate forcing exceed the differences caused by the interpolation method. Some GCMs produced unrealistic LGM volumes, and only four resulted in reasonable ice sheets, with LGM Northern Hemisphere sea level contribution ranging between 74–113 m with respect to the present day. The glacial index and climate matrix methods result in similar ice volumes at the LGM but yield a different ice evolution with different ice domes during the inception phase of the glacial cycle and different sea level rates during the deglaciation phase. The temperature–albedo feedback is the main cause of differences between the glacial index and climate matrix methods.

1 Introduction

Sea level rise due to the melt of the Greenland and Antarctic ice sheets is one of the biggest threats posed by anthropogenic climate change (Fox-Kemper et al., 2021). Ice sheets have a substantial influence on the climate system, as these can amplify changes in temperatures and alter precipitation patterns, which in turn affect ice accumulation and ablation. It is therefore important to make accurate projections of future sea level change that include the interactions between the climate and the ice sheets. These interactions are important on long timescales because ice sheets respond generally slowly to changes in temperature and precipitation (Clark et al., 1999), although brief periods of rapid ice loss have occasionally occurred in the geological past (e.g. Gomez et al., 2015; Brendryen et al., 2020). Direct observations are insufficient to study ice sheet–climate interaction as these only capture the changes in ice sheets over the past century. Instead, the palaeo-record provides information on the climate system, such as ice sheet extent and thickness, atmospheric CO2 concentrations, and eustatic sea level, dating back well before modern observations. These palaeo-reconstructions allow the study of considerable changes in the ice sheet and climate that took place over multiple millennia. The most recent period in the Earth's history with substantial changes in ice sheet extent is the last glacial cycle (LGC; 120 000–8000 years ago). This period is associated with a sea level change of approximately 120 m (e.g. Simms et al., 2019; Gowan et al., 2021) and a decrease in global temperature of 4–6 K (e.g. Annan et al., 2022; Tierney et al., 2020) with respect to present day.

The climate during the LGC is affected by several internal and external processes in the climate system. Over long timescales, insolation changes due to orbital parameters are significant enough to affect climate and ice sheets (Löfverström et al., 2014). Additionally, atmospheric CO2 concentrations changed between 190 and 280 ppm during the LGC, also acting as a forcing to the climate system. Topography and albedo changes result from a change in ice thickness and extent. Changes in albedo and topography affect temperature and precipitation (e.g. Abe-Ouchi et al., 2007; Clark et al., 1999). Especially the temperature–albedo and precipitation–topography interactions induced by changes in ice volume have a substantial impact on ice sheets (Abe-Ouchi et al., 2007; Stap et al., 2014). For example, a decrease in temperature prompts an increase in snow coverage and thereby albedo. This increase in albedo decreases the total thermal energy stored in the climate system, as a larger portion of solar radiation entering the atmosphere is reflected towards space. Therefore, the change in albedo causes a positive feedback, enhancing temperature change on regional and global spatial scales. In addition, local temperatures are also affected by changes in surface topography due to the atmospheric temperature lapse rate.

In addition, topography has a local and regional impact on precipitation. Precipitation is enhanced on the slopes of ice margins, since cooling and condensation take place when air is lifted from the margin to the inland ice. During this transport, air cools, and moisture is removed, resulting in low precipitation on the lee side and on ice plateaus. Changes in surface topography during glacial cycles can also affect atmospheric circulation (e.g, Löfverström et al., 2016), again affecting temperature and precipitation patterns (Pausata et al., 2011; Ullman et al., 2014). These feedbacks, which act over multiple millennia as the ice sheet gradually incepts, grows, and retreats, must be accounted for during transient climate and ice sheet simulations.

An ideal set-up to transiently simulate ice sheet–climate interactions during the LGC would involve a fully coupled general circulation model (GCM) that simulates ice sheets, oceans, and atmosphere. However, these simulations require a large amount of computation time, making them currently unfeasible. One strategy to deal with this excessive computational demand is to decrease grid resolution, use asynchronous coupling while certain physical processes are artificially accelerated (e.g. Smith and Gregory, 2012), or use Earth system models of intermediate complexity which have a reduced computational cost compared to GCM's but can still explicitly simulate ice sheet–climate interactions (Ganopolski et al., 2010). Alternatively, when simulating specific periods in time such as the Last Glacial Maximum (LGM; 21 000 years ago) and pre-industrial (PI), GCMs can be used. Modelling efforts such as the Paleoclimate Modelling Intercomparison Project (PMIP) intercompare a set of GCMs that used similar boundary conditions, such as atmospheric CO2 concentration, orbital configuration, and fixed prescribed ice sheet geometry, to simulate the climate for a specific time slice. The third phase of PMIP, PMIP3 (Braconnot et al., 2011), used nine GCMs to simulate climates during the LGM and PI. Each climate model used prescribed ice sheets by Abe-Ouchi et al. (2015); hence transient changes in ice topography were not simulated. Despite using the same boundary conditions, ice sheet model studies using the PMIP3 models such as Niu et al. (2019) and Alder and Hostetler (2019) found substantial differences in LGM ice volume and extent. In these studies, Niu simulated the entire last glacial cycle, while Adler and Hostetler (2019) used steady-state LGM climate forcing, both finding large differences between the ice sheets resulting from the GCM climates, some showing unrealistic ice sheets.

Simulating a full glacial cycle using a standalone ice sheet model is feasible as these models require only a small amount of runtime compared to GCMs. GCMs model the entire globe and simulate processes that require small time steps. Ice sheet models can run at a much lower temporal resolution and only require a grid that covers a small portion of the Earth, but they need a higher spatial resolution.

In the past different methods have been used to include transiently changing climate forcing to an ice sheet model. These include applying a globally uniform temperature anomaly directly to present-day climate (de Boer et al., 2014) by interpolating between LGM, mid-Holocene, and PI surface mass balances from a climate model (Fyke et al., 2014); by coupling the ice sheet to a simple energy balance model (Stap et al., 2014); using Earth system models of intermediate complexity (Ganopolski et al., 2010); or by using a glacial index method (Niu et al., 2019). With the glacial index method LGM and PI temperature and precipitation are transiently interpolated with respect to one parameter obtained from palaeo-reconstructions, such as CO2. Therefore, as CO2 decreases, the climate cools and dries as it approaches the LGM. This method has been shown to yield LGM ice volume and extent that agrees well with reconstructions (e.g, Charbit et al., 2007; Niu et al., 2019). However, this interpolation does not include any ice sheet–climate interactions, other than (usually) a lapse-rate-based linear elevation correction. To parameterise these feedbacks, we can use a climate matrix method instead (Pollard, 2010). For this we need results from at least two time slices, although more will provide better constraints. Hence with a climate matrix method the ice sheet evolution affects the dynamics of the climate and ice system, whereas with a glacial index the ice sheet only responds passively to the climate forcing. The climate matrix method implicitly resolves the temperature–albedo and precipitation–topography feedbacks. Here, the climate time slices are interpolated with respect to both prescribed forcing, for example atmospheric CO2 concentrations, and the internally calculated albedo and ice sheet geometry. Since it incorporates calculated fields from the ice sheet model into the interpolation, a climate matrix method can implicitly resolve ice sheet–climate feedbacks requiring only a small amount of additional computational resources. When performing realistic scenarios of the LGC, the climate matrix method was shown to be able to successfully replicate ice evolution (Berends et al., 2018).

While the climate matrix method and glacial index method have been used in the past, a comparison with a realistic scenario for the LGC has not yet been explored. Ladant et al. (2014) simulated ice sheets during the Eocene–Oligocene transition, and Abe-Ouchi et al. (2013) investigated Pleistocene glacial cycles by interpolating between snapshots with varying ice sheet sizes, orbital parameters, and atmospheric CO2. Recently, Stap et al. (2022) compared a matrix and index method with schematic experiments of Antarctica during the Miocene. This showed that the temperature–albedo and precipitation–topography feedback operate in opposite ways: the temperature–albedo feedback substantially reduces and the precipitation–topography feedback slightly increases glacial–interglacial variability. They suggested that the temperature–albedo feedback is stronger than the precipitation–topography feedback.

Here we aim to build upon the work by Stap et al. (2022) by intercomparing a glacial index and climate matrix method using a realistic experiment by simulating the Northern Hemisphere ice sheets during the LGC. We force the ice sheet model using different available climate model output from the PMIP3 ensemble. This study uses a climate matrix method and therefore also builds on Niu et al. (2019), who simulated the LGC using PMIP3 climate forcing interpolated with a glacial index method. Our ice sheet model and climate forcing set-up are described in Sect. 2. Section 3 introduces our simulations of the LGC and shows the differences due to climate forcing and due to the interpolation methods. These findings are discussed in Sect. 4.

2 Methods

2.1 Ice sheet model

In this study, we use the three-dimensional thermodynamically coupled ice sheet model IMAU-ICE version 2.0 (Berends et al., 2022). This model uses the depth-integrated viscosity approximation (DIVA; Goldberg, 2011) to calculate the dynamics of floating and grounded ice. This vertically integrated approximation to the stress balance is similar to the hybrid shallow ice–shallow shelf approximation (Bueler and Brown, 2009) but has improved physics, is more efficient, and has improved numerical stability (Robinson et al., 2022). A regularised Coulomb sliding law is used to calculate basal friction (Bueler and van Pelt, 2015). Proper grounding-line migration is achieved using a sub-grid friction-scaling scheme, which is based on the approach used in the Parallel Ice Sheet Model (PISM; Feldmann et al., 2014) and the Community Ice Sheet Model (CISM; Leguy et al., 2021). Sub-shelf melt rates are calculated using a temperature and depth-dependent sub-shelf melt parameterisation (Martin et al., 2011) in combination with parameterised, globally uniform ocean temperature changes (de Boer et al., 2013). As a result, while computationally fast, this parameterisation does not include ocean temperature fields from the GCM simulations and also does not capture the spatial pattern in ocean temperature changes. Bedrock adjustment to changes in ice load is modelled using an elastic lithosphere–relaxing asthenosphere model (Le Meur and Huybrechts, 1996). Calving is not included in this version of the model.

We simulate the North American, Eurasian, and Greenland ice sheets concurrently in three separate domains with a 40×40, 40×40, and 20×20 km resolution respectively (see Fig. 1). We use a higher resolution for Greenland as the ice sheet is smaller, allowing us to capture the effect of small topographical changes with a comparative number of grid cells compared to North America and Eurasia. To prevent double-counting of ice, ice growth is prevented in the Greenlandic parts of the North American and Eurasian domains and Ellesmere Island in the Greenland domain. Antarctica is not included in these simulations as the feedbacks from topography and albedo throughout the LGC are small compared to Eurasia and North America, due to the relatively small change in extent and elevation. In reality, Antarctica may have a substantial impact on ocean circulation and the carbon cycle (Adkins, 2013) and thereby climate evolution, but these feedbacks can only be captured in fully coupled models and not with the parameterised forcing applied here.

Figure 1The North American (red), Greenland (green), and Eurasian (blue) domains used in the ice sheet model. The LGM extent from Abe-Ouchi et al. (2015) is shown in black. This ice sheet reconstruction is used as a boundary condition in the PMIP3 climate model simulations.

2.2 Surface mass balance

The monthly surface mass balance (SMB) is calculated from the monthly temperature and precipitation resulting from the climate interpolation (see Sect. 2.4) using an insolation–temperature model (IMAU-ITM; Berends et al., 2018), which has been part of the Greenland Surface Mass Balance Model Intercomparison Project (GrSMBMIP; Fettweis et al., 2020). In this parameterised SMB scheme, snow accumulation is calculated using a precipitation- and temperature-dependent snow–rain partitioning by Ohmura et al. (1999). This snow–rain partitioning was tuned towards the regional climate model RACMO as part of GrSMBMIP. Annual refreezing is calculated using the approach by Huybrechts and de Wolde (1999) and Janssens and Huybrechts (2000). Ablation is parameterised and depends on temperature, insolation, and surface albedo (Bintanja et al., 2002). The insolation at the top of the atmosphere is prescribed and obtained from Laskar et al. (2004). Surface albedo is calculated internally: first a snow-free albedo is applied, with 0.1 for ocean, 0.2 for land, and 0.5 for bare ice. A firn layer is added on top, which, depending on depth, can increase the albedo to a maximum of 0.85. The planetary albedo changes due to cloud cover changes are not taken into account. The equations governing IMAU-ITM are presented in Appendix A.

2.3 PMIP3 climate time slices

We obtained climate forcing for near-surface air temperature, precipitation, and topography from the nine available LGM and PI GCM simulations that are part of PMIP3. While recently several LGM and PI simulations for PMIP4 have become available, we nevertheless decided to use PMIP3. PMIP4 is not fundamentally different to PMIP3 (Kageyama et al., 2021). However, PMIP3 has been more widely studied, allowing for comparisons with studies that have been conducted in the past such as Niu et al. (2019) and Alder and Hostetler (2019). The PMIP3 simulations are listed in Table 1. Each of these models used identical boundary conditions following the PMIP3 protocol, such as orbital parameters, trace gases, and ice sheets. The PMIP3 protocol contains a prescribed ice sheet by Abe-Ouchi et al. (2015), which is a composite of the ICE-6G v2.0 (Argus and Peltier, 2010), GLAC-1a (Tarasov et al., 2012), and ANU (Lambeck et al., 2010). We selected a subset of the nine available PMIP3 LGM and PI simulations to obtain a climate forcing that results in good agreement with reconstructions. This step is described in Appendix B. The selection of acceptable models consists of COSMOS, IPSL, MIROC, and MPI. Our parameterised SMB scheme was tuned to the mean climate of this sub-selection, which is hereafter referred to as PMIP3-Ensemble.

Table 1The climate forcing from PMIP3 model output that was used in this study. PMIP3-Ensemble represents the mean of COSMOS, IPSL, MIROC, and MPI. The global annual temperature (T) and precipitation (P) difference between the LGM and PI is shown for each climate model.

Download Print Version | Download XLSX

To apply the GCM fields to the ice sheet model, corrections need to be applied to account for the difference in resolution and topography. First a correction is applied to deal with biases in the GCM models, which is described in Appendix C. Secondly, to account for the large difference in resolution, the GCM fields are bilinearly interpolated to the ice sheet model grid. Thirdly, since the ice sheet model topography evolves during the LGC, which is not accounted for in the GCM climate time slices, a topographic correction needs to be applied. Here we use the approach by Berends et al. (2018), which implies that temperature is downscaled using a dynamic lapse rate correction. Precipitation is downscaled using the Roe and Lindzen (2001) model. This model accounts for the precipitation changes as a result of large-scale changes in topography and takes into account the surface slope, wind direction, and changes in surface height. Throughout the LGC, Greenland's topography change is less substantial compared to North America and Eurasia. Therefore, instead of the Roe and Lindzen model, we apply a Clausius–Clapeyron relation that adapts the applied precipitation based on the change in surface temperature.

2.4 Transiently changing forcing

In this study we use two methods to transiently interpolate between climate time slices: a glacial index method and a climate matrix method. With our glacial index method, precipitation and temperature fields are interpolated with respect to prescribed CO2 obtained from Bereiter et al. (2015). Figure 2 depicts the atmospheric CO2 concentration during the LGC and the corresponding values for the glacial index. A glacial index of 1 (0) represents full glacial (interglacial) conditions. Hence, the climate forcing is equal to the LGM (PI). Temperature and precipitation are respectively interpolated linearly and logarithmically. The logarithmical interpolation for precipitation prevents negative values and is used to describe the relative changes during the LGC.

Figure 2Atmospheric CO2 concentrations during the LGC from Bereiter et al. (2015) and the corresponding glacial index values. A glacial index value of 1 represents “full glacial conditions” corresponding to an LGM climate. A glacial index value of 0, or “full interglacial conditions”, represents a PI climate.

The climate matrix method expands upon the glacial index method by implicitly resolving the temperature–albedo and precipitation–topography feedbacks. This is achieved by making the interpolation parameter spatially variable. In this study, we applied the climate matrix method using the approach by Berends et al. (2018). Temperature is interpolated with respect to CO2 and to the absorbed insolation at the surface. The absorbed insolation is computed using the internally calculated albedo and the insolation solution by Laskar et al. (2004). Therefore, any modelled advance of the ice will result in an increased albedo, an increase in the interpolation weight, and therefore a more glacial climate. Precipitation is interpolated with respect to the change in topography between the LGM and PI. Regions without LGM ice cover do not undergo substantial topographical changes. In these regions precipitation is interpolated with respect to the total changes in topography throughout the domain. Equations governing the glacial index and climate matrix methods are presented in Appendix D. Neither the climate matrix nor the glacial index method includes interactions between the ocean and ice or climate. However, changes in the ocean circulation had substantial influence on the climate and the carbon cycle (Adkins, 2013; Toucanne et al., 2021) and affected sub-shelf melt rates. Atmospheric circulation changes due to ice sheet topography changes are also not explicitly modelled. Nevertheless, the climate matrix method is a computationally fast method to implicitly resolve ice–temperature and topography–precipitation feedbacks without transient GCM simulations. This method should therefore be seen as an alternative to the glacial index method, not an alternative to climate models.

3 Results

3.1 Climate forcing

In this section, we compare simulations of the LGC resulting from different climate forcings using the climate matrix method. Figure 3 shows the sea level contribution during the LGC with the climate forcing from COSMOS, IPSL, MIROC, MPI, and the PMIP3-Ensemble. The sea level contribution at the LGM ranges substantially between the five simulations. As shown in Fig. 3a, the LGM sea level contribution of the Northern Hemisphere ice sheets ranges between 74 m (COSMOS) and 119 m (MIROC) of sea level equivalent (m s.l.e.). When including the contribution for Antarctica of 10 m s.l.e. (Simms et al., 2019), the total LGM global sea level contribution is within the range of Simms et al. (2019) for the PMIP3-Ensemble, MPI, and IPSL.

Figure 3Sea level contribution of the Northern Hemisphere ice sheets using the climate matrix method and the uncertainty range of LGM ice volume from Simms et al. (2019). Each simulation was forced with a climate obtained from a member of the PMIP3 ensemble.

Figure 4The time at which a region was first covered by ice during the LGC. Darker colours represent earlier inception. Lighter-colour regions were covered by ice later. The black contour depicts the ice extent reconstruction Abe-Ouchi et al. (2015) used in the GCM simulations. Each simulation was forced by GCM output from PMIP3 interpolated using the climate matrix method.

Figure 5Summer (JJA) temperature differences between the LGM and PI. IPSL and MPI have lower temperatures in Arctic Canada compared to the other GCM climates, which results in larger rates of ice growth at the onset of the LGC. Black contours represent the LGM extent of the Abe-Ouchi et al. (2015) ice sheets. The LGM extent of the ice sheet model simulations is shown in green.

The North American ice sheet is the largest contributor to the differences in ice volume, with a minimum LGM contribution of 56 m s.l.e. (COSMOS) and a maximum of 87 m s.l.e. (MIROC; Fig. 3b). For the Eurasian ice sheet, LGM sea level contribution varies between 12 m (IPSL) and 29 m (MIROC) (Fig. 3c). This substantial difference between ice sheet model runs forced by individual PMIP3 members is in line with findings by Niu et al. (2019) and Alder and Hostetler (2019).

Figure 4 shows when each region was covered by ice for the first time during the LGC, thereby indicating the timing of inception and the gradual increase in extent of the ice sheets. For example, dark-purple areas over Ellesmere Island and Baffin Island show the immediate inception of the ice sheets in those regions. They are currently covered by ice caps, and therefore near-pre-industrial temperatures are sufficient for ice growth. The light-orange colours in southern Scandinavia indicate that ice covered these regions from ∼30 ka onwards. In North America, the ice sheets incept from the North American Cordillera and eastern Canada, forming two large ice domes. In Eurasia, ice sheets incept at the islands surrounding the Barents Sea. Little evidence exists on the inception phase of the LGC, since later glaciations tend to have removed most geological evidence. Modelling studies based on geological evidence, such as Bahadory et al. (2021) and Dalton et al. (2022), suggest that the North American ice sheet incepted at Ellesmere Island and the Canadian Cordillera, which agrees reasonably well for most GCM forcings except for the COSMOS forcing.

Between 120 and 60 ka, simulations forced with MPI and IPSL output have a comparatively large sea level contribution (Fig. 3a). This can also be seen in the ice evolution (Fig. 4) for MPI and IPSL forcing, leading to a comparatively large extent of the Baffin and Innuitian ice sheets during the inception phase of the ice sheets. Similarly, MPI forcing resulted in a large extent of the Barents–Kara ice sheets at the onset of the LGC. Since each simulation was forced by the same prescribed CO2 and insolation, this is a result of the climate forcing and internal feedback processes.

Figure 5 shows the temperature difference between LGM and PI. Shown here is that regions with large ice extent at the early parts of the simulations correspond generally well to large LGM–PI temperature differences promoting ice growth. The LGM and PI temperatures are linearly interpolated with respect to prescribed CO2 and albedo. Therefore, the same change in CO2 and albedo results in a larger temperature change with increasing LGM–PI temperature difference provided by the GCMs.

Figure 6 depicts ice thickness at the LGM, showing that the ice extent varies considerably with climate forcing. In North America, the differences in extent are mostly located along the southern margins. As the ice sheet forced with MIROC has an unrealistically large extent compared to reconstructions, the LGM ice sheet exceeds the southern domain boundary. In the Eurasian ice sheet, the differences are found in western Europe. Although in reality the British Isles were covered by glacial ice at the LGM (e.g. Abe-Ouchi et al., 2015; Batchelor et al., 2019), this only occurs in the simulation forced by the MIROC climate, though this is accompanied with a large ice volume compared to the reconstructions by Simms et al. (2019) (see Fig. 3d). IPSL and COSMOS forcings lead to limited LGM ice extent south of the Scandinavian mountains, which does not match palaeo-reconstructions (e.g. Abe-Ouchi et al., 2015; Batchelor et al., 2019). Figure 7 shows LGM temperature and compares it to the extent of the ice sheet simulations. This figure indicates that low LGM summer temperatures tend to match ice extent well, which is in line with Niu et al. (2019). This is to be expected considering that temperature is important for the SMB by affecting the amount of melt, refreezing, and snowfall.

Figure 6Ice thickness using forcing obtained from members of the PMIP3 ensemble. The ice extent by Abe-Ouchi et al. (2015) is shown as a black contour.

Figure 7LGM summer temperature as it is applied to the ice sheet model. The green contours indicate the extent of the ice sheets that resulted from the GCM forcing. The black contours show the extent of the reconstruction by Abe-Ouchi et al. (2015).

Figure 8Ice thickness (a, b) and timing of first ice (c, d) of the climate matrix (a, c) and glacial index (b, d) methods. The black contour represents the extent of the ice reconstruction by Abe-Ouchi et al. (2015).

We are able to capture the LGM sea level (Simms et al., 2019) and extent (Abe-Ouchi et al., 2015) well using the climate matrix method. While we tuned ice volume towards Simms et al. (2019), the extent of the LGM ice sheet is smaller compared to Abe-Ouchi et al. (2015). However, the rate of ice growth during Marine Isotope Stage (MIS) 5, as well as peak ice volumes at 60 ka, is not captured well (e.g. Gowan et al., 2021). This may be a result of a too-small effect of insolation on temperature and SMB. The interior of the Laurentide ice sheet in our simulations has a thickness over 4000 m, exceeding the thickness of the ICE6G-C reconstruction by Peltier et al. (2015). This suggests that not enough ice is transported from the interior towards the margin, resulting in a small extent and large thickness in the interior, possibly a result of shortcomings in the basal sliding and specifically a too-high basal friction.

3.2 Ice sheet–climate feedbacks

Here we investigate the effect of the albedo–temperature and topography–precipitation feedbacks. We present two simulations using the PMIP3-Ensemble forcing, which differ in their use of either the glacial index or the climate matrix method.

Figure 8c and d compare when the first ice accumulates in regions for the glacial index and climate matrix methods. In the Eurasian ice sheet, which incepts mostly in the Barents–Kara Sea region, more domes are formed with the glacial index method. This difference in the number of domes is even more pronounced in North America. Using the climate matrix method, only the Laurentide and Cordilleran ice domes develop, which merge around 40 ka. With the glacial index method, many smaller domes are formed in the Keewatin, Baffin, and Cordilleran regions. The ice domes in the glacial-index-method simulations have much more irregular shapes compared to the smoother margins of the climate matrix ice domes. This North American ice sheet with few inception regions with the climate matrix method agrees better with studies conducted with geological constraints (e.g. Batchelor et al., 2019; Dalton et al., 2022).

This difference in inception between the glacial index method and the climate matrix method is due to the feedback processes. An ice sheet has a high albedo and thereby creates a regionally cold climate, enhancing ice growth. In regions without snow or ice, albedo is low, leading to higher temperatures, inhibiting ice inception. This feedback process is included with the climate matrix method, which separates temperature change caused by albedo and insolation and CO2. Consequently, the albedo in the ice sheet model has a pronounced influence on temperature and thereby on the ice sheet evolution. With the glacial index method temperature is only affected by CO2 and does not account for albedo changes. Hence, the glacial index method underestimates cooling in high-albedo regions and overestimates cooling in low-albedo regions.

Figure 9Sea level contribution during the LGC for simulations forced by PMIP3-Ensemble using either the glacial index or climate matrix method. The error bars indicate the sea level contribution by Simms et al. (2019).


Figure 9 shows the sea level contribution of the Northern Hemisphere ice sheets over time. Ice volume at the LGM is slightly larger with the glacial index method (110 m s.l.e.) compared to the climate matrix method (96 m s.l.e.). Accordingly, the glacial index method has larger volumes for both the North American (Fig. 9b) and Eurasian (Fig. 9c) ice sheets. Ice thickness at the LGM is shown in Fig. 8a and b. The Eurasian ice sheet with the glacial index method has more ice covering the British Islands, but less in the Barents Sea region. The North American ice sheet margin extends more towards the south with the glacial index method. The smaller volume with the climate matrix method is mostly caused by the temperature–albedo interaction. With the climate matrix method, temperature and precipitation are only equal to the LGM time slice when CO2, albedo and topography are the same as the ice sheet reconstruction. However, the extent of the modelled LGM ice sheet is smaller compared to the reconstruction by Abe-Ouchi et al. (2015). As a result, the modelled albedo is too low in the region, and consequentially the temperature is too high. Therefore, at the LGM, the temperature in the ice sheet model is higher compared to the LGM climate forcing. With the glacial index method, when CO2 levels are equal to the LGM the climate forcing is equal to the LGM as well. Therefore, in these simulations, the LGM volume is higher with the glacial index method.

After the LGM, the Eurasian and North American ice sheets retreat. Figure 10b shows the modelled sea level contribution rate during the LGC. While sea level contribution rates reach negative values at approximately the same time, the climate matrix method finishes retreating later (5 ka) compared to the glacial index method (8 ka). Furthermore, Fig. 10a compares sea level contribution to CO2, with timestamps indicated in the figure. Considering an arbitrary threshold of 5 mm yr−1 shows that the glacial-index-method retreat rate accelerates faster at lower CO2 concentrations (18 ka, 225 ppm) compared to the climate matrix method (16 ka, 240 ppm). This indicates a lower CO2 threshold for retreat with the glacial index method. The peak sea level rate is higher and earlier with the glacial index method, with a decrease of 32 mm yr−1 (11 ka) compared to 19 mm yr−1 (8 ka) with the climate matrix method. This is because when the CO2 increases rapidly, low temperatures persist longer with the climate matrix method due to the high albedo of the large ice sheet, while the response of the glacial index method is instantaneous. The climate matrix has the tendency to maintain the ice sheet as it is and does not warm as quickly, resulting in lower retreat rates.

Figure 10Sea level contribution plotted against CO2 concentrations (a). The sea level contribution rate during the LGC (b). LI indicates last interglacial; PD indicates present day.


4 Discussion and conclusion

In this study we simulated the Northern Hemisphere ice sheets of the LGC using an ice sheet model. Our aim was to investigate the sensitivity to palaeoclimate forcing and to assess the effects of the albedo–temperature and precipitation–topography interactions.

Climate forcing is obtained from the PMIP3 ensemble to investigate the sensitivity of a climate matrix and glacial index method. We find that the differences in ice volume due to the climate forcing exceed the differences caused by the transient climate interpolation method (glacial index versus climate matrix). Several PMIP3 models yield unrealistic ice sheet configurations. Despite choosing only a subset of the PMIP3 simulations, the sea level contribution of the Northern Hemisphere ice sheets still shows a considerable range of 74–113 m. These large differences are in line with findings by Niu et al. (2019) and Alder and Hostetler (2019) and are not exclusively found with PMIP3 forcing, but also for the first PMIP ensemble by Charbit et al. (2007). Our study additionally used the climate matrix method instead of only the glacial index method used by Charbit et al. (2007) and Niu et al. (2019), showing that even when including atmospheric feedback processes, the range in sea level contribution is still large. As originally suggested by Niu et al. (2019), cold LGM summer temperatures generally correspond well with areas of LGM ice cover when using PMIP3 climate forcing to simulate the LGC. This is caused by the fact that high ablation rates at the margin inhibit ice advance, and ablation rates depend strongly on temperature.

In this study, we compared the glacial index and climate matrix method with a realistic experiment for the first time. The LGC simulations using the glacial index and climate matrix method were conducted with the same climate forcing, namely the PMIP3 ensemble mean. The glacial index method and climate matrix methods yield similar LGM ice volumes but show a much larger difference in ice evolution. Generally, the climate matrix method incepts at fewer domes, with the few domes that form being larger and more regularly shaped compared to the glacial index method. With the glacial index method, many more smaller domes form, often far away from the main ice sheet centres. With both methods, the domes gradually merge to form one big Eurasian or Laurentide ice sheet. The difference in inception is due to the albedo–temperature feedback, where temperatures with the climate matrix method are low along ice margins and high in regions with low albedo. With the climate matrix method this enhances ice growth close to ice margins and inhibits ice inception of new domes in regions with low albedo. We find that the ice sheets, especially the Laurentide ice sheet, deglaciate faster with the glacial index method compared to the climate matrix method. This is attributed to the albedo–temperature feedback, as the high albedo of the ice sheet can maintain colder temperatures at the onset of deglaciation.

LGM extent matches reconstructions well, but there are discrepancies for the British Islands. This may have been a result of our simplistic approach to sub-shelf melt rates. Furthermore, Greenland and North America are simulated in separate domains, preventing ice dynamical interactions between these ice sheets. Additionally, temperatures in the climate forcing may have been too high for ice inception in the British ice sheet, except for the simulation forced by climate from MIROC. Neither the glacial index nor climate matrix method matches global eustatic sea level evolution well during MIS5 compared to reconstructions (e.g. Gowan et al., 2021). With either method ice volume increases gradually, while studies that are based on geological constraints such as Batchelor et al. (2019) and Dalton et al. (2022) suggest a more dynamic change in ice volume and area for the North American and Eurasian ice sheets. Our method is unable to capture the maximum volume and extent of the Eurasian ice sheets at 60 ka (e.g. Gowan et al., 2021; Kleman et al., 2013). This indicates that we are missing forcing and/or feedback processes that would improve ice evolution. One of the main reasons could be a too-weak effect of insolation on temperature and SMB. With both the climate matrix and glacial index methods CO2 is the main driver of temperature changes. However, the CO2 record alone is not sufficient to capture the variability in ice evolution throughout the LGC. Insolation was also found to have a substantial impact on glacial cycles (Abe-Ouchi et al., 2013).

In addition, the climate matrix method can be expanded by including snapshots that were generated with different ice sheets and orbital parameters. With more snapshots, a large range of orbital parameters and ice sheets can be studied (e.g. Abe-Ouchi et al., 2013; Ladant et al., 2014). Here we have chosen to use a wide array of GCM simulations, limiting the availability of different snapshots such as insolation minima and maxima. As a consequence of using only LGM and PI time slices we do not capture threshold behaviour such as the closure of the Canadian Arctic Archipelago gateways (Löfverström et al., 2022), abrupt changes in atmospheric circulations in the North Atlantic during deglaciation (Löfverström and Lora, 2017), or processes that are not captured in LGM and PI time slices such as Dansgaard/Oeschger and Heinrich events and their impact on climate (Claussen et al., 2003). While we included parameterised albedo–temperature and precipitation–topography interactions, there are many more ice sheet–climate interactions that are not included. Ice sheets can affect atmospheric circulation and vice versa, which is not well resolved in the glacial index or climate matrix methods, especially given the limited number of GCM snapshots. For example, the North American ice sheet affects atmospheric circulations, which can alter the location and size of the Eurasian ice sheet (Liakka et al., 2016), with a large North American ice sheet preventing growth of ice in Siberia. Currently, this feedback is not resolved, which may have affected the evolution and shape of the Eurasian ice sheet throughout the LGC. Freshwater influx into the ocean as the Laurentide ice sheet melted may have had a large impact on the Atlantic meridional overturning circulation (Otto-Bliesner and Brady, 2010). With our method, the ocean is not explicitly modelled. Our basal melt method is simplistic and parameterised and does not use ocean data from a GCM. Therefore, this method is able to capture neither any ice–ocean interactions nor the effect of ocean circulation changes. With the climate matrix method, albedo is annually averaged, ignoring seasonal fluctuations. To explicitly simulate the above-mentioned processes, transient GCM simulations are required, which perform adequately, rather than any temporal interpolation method as used here. This requires a major investment of computational resources. The climate matrix method should therefore be viewed as an alternative to the glacial index method, not an alternative to climate models.

Appendix A: Surface mass balance model

We use IMAU-ITM to calculate surface mass balance (Berends et al., 2018). This SMB scheme requires temperature and precipitation fields that we obtain from the GCM fields after applying the methods described in Appendix C and D. Monthly melt in IMAU-ITM is parameterised using Bintanja et al. (2002) and calculated as follows:

(A1) Melt ( x , y , m ) = c 1 T ( x , y , m ) - T 0 + c 2 ( 1 - α ( x , y , m ) ) Q TOA ( x , y , m ) - c 3 s yr - 1 / L fusion 12 000 .

Here, the melting temperature of ice (T) is 273.16 K, QTOA is the insolation at the top of the atmosphere from Laskar et al. (2004) in watts per square metre, and Lfusion is the latent heat of fusion in joules per kilogramme. Values for c1, c2, and c3 are tuning parameters shown in Table A1. Albedo (α) is calculated internally, based on the approach by Bintanja et al. (2002):

(A2) α surface ( x , y , m ) = α snow - α snow - α background e - 15 FirnDepth ( x , y , m - 1 ) - 0.015 MeltPreviousYr ( x , y ) .

Here, αbackground is 0.1 for water, 0.2 for land, and 0.5 for bare ice; αsnow is 0.85; and “MeltPreviousYr” is the melt of the previous model year. The albedo is capped between the background and snow albedo. Firn depth is calculated using the amount of snowfall that has been added without melting and is capped at 10 m. For snowfall we use a temperature-based snow–rain partitioning from Ohmura et al. (1999).

(A3) SnowFraction = 0.5 1 - atan T ( x , y , m ) - T 0 3.5 1.25664

The “SnowFraction” is the fraction of precipitation (P) that falls as snow with the remainder falling as rain. To calculate the total accumulation, we need refreezing. This is calculated using the approach by Huybrechts and de Wolde (1999):


Using the previously calculated snowfall, rainfall, refreezing, and melt, the SMB can be obtained:

(A7) SMB ( x , y , m ) = Snowfall ( x , y , m ) + Refreezing ( x , y , m ) - Melt ( x , y , m ) .

Table A1The parameters used to calculate ablation in the approach by Bintanja et al. (2002). The surface mass balance parameters were tuned towards the PMIP3-Ensemble simulation to obtain ice sheets that agree well with Simms et al. (2019).

Download Print Version | Download XLSX

Appendix B: Climate forcing selection and SMB tuning

Niu et al. (2019) and Alder and Hostetler (2019) have shown that there are large differences in the modelled ice sheets between the different PMIP3 climates. For our ice sheet simulations, we would ideally use climate forcing with which we can obtain a good representation of the ice sheets at the LGM. First, we conducted simulations forced with all nine available PMIP3 simulations using an untuned ablation (see Table A1); these are shown in Fig. B1. Since we obtained large differences between the simulations forced by the individual PMIP3 ensemble members, some deviating substantially from reconstructions, we used a smaller selection of models.

To make a selection of the PMIP3 climates, we compared the simulated ice sheet extent to the reconstruction from Abe-Ouchi et al. (2015). This reconstruction was used as a prescribed LGM ice extent in the climate model simulations. We compare the ice sheet model and reconstruction at the LGM to compute the percentage of too-large and too-small ice extent. For example, the simulation with MRI forcing resulted in a small Eurasian ice sheet. Ice in the simulation does not cover 91 % of the reconstructed extent. Some ice is also found outside the extent of the reconstruction. Therefore, MRI climate forcing also resulted in 1 % too-large ice extent outside the reconstruction. These percentages for too-large and too-small ice sheets are listed in Table B1. We added these values together and applied a threshold of 40 % to select the climate forcing. MIROC, IPSL, COSMOS, and MPI stayed below this threshold for both the Eurasian and North American ice sheets. The other models are above this threshold and are excluded from further study. A climate forcing was made from the mean of this selection and called PMIP3-Ensemble (four GCMs). Using this ensemble climate, the SMB model was tuned towards sea level by Simms et al. (2019).

Figure B1LGM ice thickness of an ice sheet model forced by nine different PMIP3 climates, interpolated using the climate matrix method. Black contours represent the LGM ice sheet reconstruction by Abe-Ouchi et al. (2015).

Table B1Shown here are the percentages of ice extent that is too large or too small compared to the reconstruction by Abe-Ouchi et al. (2015). These percentages are added to obtain a quality assessment for the preliminary simulations. Lower percentages indicate a better match to Abe-Ouchi et al. (2015). Simulations that stayed below the 40 % threshold are depicted in bold.

Download Print Version | Download XLSX

Appendix C: Bias correction

We apply a correction to the climate forcing to correct for biases in the climate models. First, we calculate the difference between the PI climate from the general circulation models (GCMs) and observed climate from ERA40 (Uppala et al., 2005). These differences are then applied to the LGM and PI climate forcing.

To apply the bias correction for temperature, we first need to account for differences in topography between the model and observation data. Here we apply a lapse rate (λ) correction to calculate temperature at sea level.


Using the temperature at sea level from observations (Tobs,SL) and the climate model (TGCM,SL), we can calculate the difference between observed and modelled temperature.

(C3) T GCM , bias ( x , y , m ) = T GCM , SL ( x , y , m ) - T obs , SL ( x , y , m )

TGCM,bias represents the bias between modelled and observed data. We can apply this to the LGM and PI time slices to obtain a bias-corrected temperature.

(C4) T GCM , corr ( x , y , m ) = T GCM ( x , y , m ) - T GCM , bias ( x , y , m )

Here, TGCM is the GCM output, either the PI or LGM, while TGCM,corr is the bias-corrected GCM data.

For precipitation (P), we use the ratio between model and observed data instead. First, we calculate the bias between model and observational data directly.

(C5) P GCM , bias ( x , y , m ) = P GCM , PI ( x , y , m ) / P obs , PI ( x , y , m )

Secondly, we apply this bias correction to the GCM data.

(C6) P GCM , corr ( x , y , m ) = P GCM ( x , y , m ) / P GCM , bias ( x , y , m )

After applying these corrections, we have obtained bias-corrected temperature and precipitation fields for the LGM and PI.

Appendix D: Climate interpolation

The glacial index and climate matrix methods are used to interpolate between the climate model LGM and PI time slices. As a result, both interpolation methods produce a transiently changing climate throughout the LGC. Here we apply the method by Berends et al. (2018), for which we use the bias-corrected fields for temperature and precipitation. To compute the interpolated temperature and precipitation forcing with both the glacial index and climate matrix method, we use the following two equations:


Here Tref and Pref are the temperature and precipitation climate forcing respectively. The weights for the interpolation are represented by wtot, wLGM, and wPI. The calculation for these interpolation weights determines the difference between the climate matrix and glacial index method.

With the glacial index method wtot, wLGM, and wPI are calculated with only prescribed CO2.

(D3) w CO 2 = CO 2 - CO 2 LGM / CO 2 PI - CO 2 LGM

Here, CO2 represents the prescribed CO2 concentration obtained from Bereiter et al. (2015) at the current model time step. LGM and PI CO2 concentrations are 190 and 280 ppm respectively. Note here that wCO2 is a scalar. With the glacial index method, wCO2 is equal to wtot, wPI, and (1−wLGM). Since CO2 is the same across the entire globe, not only are these values uniform, but they are also the same for each domain. There is no interaction between the ice sheet model and the interpolation.

The climate matrix method expands on the glacial index method by including a parameterised temperature–albedo and precipitation–topography feedback. For temperature, wtot is calculated using both prescribed CO2 and absorbed insolation. The annual absorbed insolation is calculated as follows:

(D4) I abs ( x , y ) = m = 1 12 Q TOA ( x , y , m ) 1 - α surface ( x , y , m ) .

Here QTOA is the insolation at the top of the atmosphere, for which we use prescribed forcing from Laskar et al. (2004). Albedo (α) is calculated internally; more details on this are described in Appendix A. To calculate annual Iabs, we use the sum of the monthly Iabs. We can use this to calculate the interpolation weights for the absorbed insolation.

(D5) w ins ( x , y ) = I abs ( x , y ) - I abs , LGM ( x , y ) / I abs , PI ( x y ) - I abs , LGM ( x y )

Iabs,PI and Iabs,LGM are the reference PI and LGM absorbed insolation respectively. This reference albedo is calculated using Eq. (A2), where we use the PI and LGM for QTOA and αsurface. To calculate albedo at the LGM, we apply the topography, ice mask, and land mask from Abe-Ouchi et al. (2015) and run the SMB model to obtain the albedo change due to snow coverage. Since albedo and insolation vary across the domain, wins is not uniform. To obtain the regional effect of albedo on temperature, we calculate the mean wins (wins,av) and wins with a Gaussian smoothing of 200 km (wins,smooth). Using these insolation weights, we can calculate the interpolation weight (wice) resulting from ice interactions.

In the North America and Eurasia domain,


In the Greenland domain, we employed a slightly different method. In Eurasia and North America, the change in albedo is mainly due to a change in snow coverage and the increase in ice extent. The Greenland ice sheet is covered by ice during both the PI and LGM. Albedo changes due to the expansion of ice shelves. Therefore, we interpolate with respect to the total change in albedo for the Greenland domain instead.


The field wtot is used in Eq. (14) to calculate temperature.

Precipitation with the climate matrix method is only dependent on surface topography. First, we calculate the total change in surface topography (Hs) between the LGM and PI.

(D10) w tot , P = H s ( x , y ) - H s , PI ( x , y ) / H s , LGM ( x , y ) - H s , PI ( x , y )

For regions within the ice extent of Abe-Ouchi et al. (2015), we calculate the regional change in topography.

(D11) w LGM , H s ( x , y ) = H s , ISM ( x , y ) - H s , GCM , PI ( x , y ) H s , GCM , LGM ( x , y ) - H s , GCM , PI ( x , y ) w tot , P ( x , y )

In regions that are outside the bounds of this ice extent, wLGM,Hs is equal to wtot,P. We multiply wLGM by wtot to account for the regional effect on precipitation due to a change in surface topography.

(D12) w LGM ( x , y ) = w LGM , H s ( x , y ) w tot , P ( x , y )

This wLGM is used in Eq. (D2) to calculate precipitation. Note here that precipitation is not uniform. Only as the ice sheet increases in thickness does Hs increase, which raises wLGM.

Code availability

The source code of IMAU-ICE is maintained on GitHub at (Scherrenberg, 2023). The exact version used in this study (including makefiles, compiling scripts, run scripts, config files for all the simulations presented here, and MATLAB scripts for creating the figures) is available at Zenodo (; Scherrenberg et al., 2022a). Please note that model simulations cannot be conducted without input files for CO2, climate and initial topography. For more information, contact the corresponding author.

Data availability

The simulation output shown in this study is available for a 10 kyr output frequency at Zenodo (; Scherrenberg et al., 2022b). Output every 1 kyr as well as additional fields can be requested by contacting the corresponding author.


The supplement related to this article is available online at:

Author contributions

MDWS conducted the simulations and wrote the manuscript. Experimental set-up was created by RSWvdW, CJB, and MDWS. Model support was provided by CJB and LBS. All authors helped analyse the results and contributed to the text.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The use of supercomputer facilities was sponsored by the Dutch Research Council (NWO) Exact and Natural Sciences. Model runs were performed on the Dutch National Supercomputer Snellius. We would like to acknowledge SurfSARA Computing and Networking Services for their support. We would finally like to thank the editor and two anonymous reviewers.

Financial support

Meike D. W. Scherrenberg is funded by the Netherlands Earth System Science Centre (NESSC), which is supported financially by the Ministry of Education, Culture and Science through OCW grant no. 024.002.001. Constantijn J. Berends was supported by PROTECT. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 869304. Lennert B. Stap is funded by the NWO through VENI grant VI.Veni.202.031.

Review statement

This paper was edited by Marisa Montoya and reviewed by two anonymous referees.


Abe-Ouchi, A., Segawa, T., and Saito, F.: Climatic Conditions for modelling the Northern Hemisphere ice sheets throughout the ice age cycle, Clim. Past, 3, 423–438,, 2007. 

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, 500, 190–193,, 2013. 

Abe-Ouchi, A., Saito, F., Kageyama, M., Braconnot, P., Harrison, S. P., Lambeck, K., Otto-Bliesner, B. L., Peltier, W. R., Tarasov, L., Peterschmitt, J.-Y., and Takahashi, K.: Ice-sheet configuration in the CMIP5/PMIP3 Last Glacial Maximum experiments, Geosci. Model Dev., 8, 3621–3637,, 2015. 

Adkins, J. F.: The role of deep ocean circulation in setting glacial climates, Paleoceanography, 28, 539–561,, 2013. 

Alder, J. R. and Hostetler, S. W.: Applying the Community Ice Sheet Model to evaluate PMIP3 LGM climatologies over the North American ice sheets, Clim. Dynam., 53, 2807–2824,, 2019. 

Annan, J. D., Hargreaves, J. C., and Mauritsen, T.: A new global surface temperature reconstruction for the Last Glacial Maximum, Clim. Past, 18, 1883–1896,, 2022. 

Argus, D. F. and Peltier, W. R.: Constraining models of postglacial rebound using space geodesy: a detailed assessment of model ICE-5G (VM2) and its relatives, Geophys. J. Int., 181, 697–723,, 2010. 

Bahadory, T., Tarasov, L., and Andres, H.: Last glacial inception trajectories for the Northern Hemisphere from coupled ice and climate modelling, Clim. Past, 17, 397–418,, 2021. 

Batchelor, C. L., Margold, M., Krapp, M., Murton, D. K., Dalton, A. S., Gibbard, P. L., Stokes, C. R., Murton, J. B., and Manica, A.: The configuration of Northern Hemisphere ice sheets through the Quaternary, Nat. Commun., 10, 1–10, 2019. 

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J.: Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present, Geophys. Res. Lett., 42, 542–549,, 2015. 

Berends, C. J., de Boer, B., and van de Wal, R. S. W.: Application of HadCM3@Bristolv1.0 simulations of paleoclimate as forcing for an ice-sheet model, ANICE2.1: set-up and benchmark experiments, Geosci. Model Dev., 11, 4657–4675,, 2018. 

Berends, C. J., Goelzer, H., Reerink, T. J., Stap, L. B., and van de Wal, R. S. W.: Benchmarking the vertically integrated ice-sheet model IMAU-ICE (version 2.0), Geosci. Model Dev., 15, 5667–5688,, 2022. 

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. 

Braconnot, P., Harrison, S. P., Otto-Bliesner, B. L., Abe-Ouchi, A., Jungclaus, J. H., and Peterschmitt, J.-Y.: The Paleoclimate Modeling Intercomparison Project contribution to CMIP5, CLIVAR Exchanges, 56, 15–19, 2011. 

Brady, E. C., Otto-Bliesner, B. L., Kay, J. E., and Rosenbloom, N.: Sensitivity to Glacial Forcing in the CCSM4, J. Climate, 26, 1901–1925,, 2013. 

Brendryen, J., Haflidason, H., Yokoyama, Y., Haaga, K. A., and Hannisdal, B.: Eurasian Ice Sheet collapse was a major source of Meltwater Pulse 1A 14,600 years ago, Nat. Geosci., 13, 363–368,, 2020. 

Budich, R., Gioretta, M., Jungclaus, J., Redler, R., and Reick, C.: The MPI-M Millennium Earth System Model: An assembling guide for the COSMOS configuration, Tech. rep., Max-Planck Institute for Meteorology, Hamburg, Germany, (last access: 2 February 2023), 2010. 

Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res., 114, F03008,, 2009. 

Bueler, E. and van Pelt, W.: Mass-conserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635,, 2015. 

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. 

Clark, P. U., Alley, R. B., and Pollard, D.: Northern Hemisphere Ice-Sheet Influences on Global Climate Change, Science, 286, 1104–1111,, 1999. 

Claussen, M., Ganopolski, A., Brovkin, V., Gerstengarbe, F.-W., and, Werner, P.: Simulated global-scale response of the climate system to Dansgaard/Oeschger and Heinrich events, Clim. Dynam., 21, 361–370, 2003. 

Dalton, A. S., Stokes, C. R., and Batchelor, C. L.: Evolution of the Laurentide and Innuitian ice sheets prior to the Last Glacial Maximum (115 ka to 25 ka), Earth Sci. Rev., 224, 103875,, 2022. 

de Boer, B., van de Wal, R., 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. 

Dufresne, J.-L., Foujols, M.-A., Denvil, S., Caubel, A., Marti, O., Aumont, O., Balkanski, Y., Bekki, S., Bellenger, H., Benshila, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Cheruy, F., Codron, F. F., Cozic, A., Cugnet, D., de Noblet, N., Duvel, J.-P., Ethé, C., Fairhead, L., Fichefet, T., Flavoni, S., Friedlingstein, P., Grandpeix, J.-Y., Guez, L., Guilyardi, E., Hauglustaine, D., Hourdin, F., Idelkadi, A., Ghattas, J., Joussaume, S., Kageyama, M., Krinner, G., Labetoulle, S., Lahellec, A., Lefèbvre, M.-P., Lefèvre, F., Lévy, C., Li, Z. X., Lloyd, J., Lott, F., Madec, G., Mancip, M., Marchand, M., Masson, S., Meurdesoif, Y., Mignot, J., Musat, I., Parouty, S., Polcher, J., Rio, C., Schulz, M., Swingedouw, D., Szopa, S., Talandier, C., Terray, P., and Viovy, N.: Climate change projections using the IPSL-CM5 Earth System Model: from CMIP3 to CMIP5, Clim. Dynam., 40, 2123–2165,, 2013. 

Feldmann, J., Albrecht, T., Khroulev, C., Pattyn, F., and Levermann, A.: Resolution-dependent performance of grounding line motion in a shallow model compared to a full-Stokes model according to the MISMIP3d intercomparison, J. Glaciol., 60, 353–360, 2014. 

Fettweis, X., Hofer, S., Krebs-Kanzow, U., Amory, C., Aoki, T., Berends, C. J., Born, A., Box, J. E., Delhasse, A., Fujita, K., Gierz, P., Goelzer, H., Hanna, E., Hashimoto, A., Huybrechts, P., Kapsch, M.-L., King, M. D., Kittel, C., Lang, C., Langen, P. L., Lenaerts, J. T. M., Liston, G. E., Lohmann, G., Mernild, S. H., Mikolajewicz, U., Modali, K., Mottram, R. H., Niwano, M., Noël, B. P. Y., Ryan, J. C., Smith, A., Streffing, J., Tedesco, M., van de Berg, W. J., van den Broeke, M. R., van de Wal, R. S. W., van Kampenhout, L., Wilton, D., Wouters, B., Ziemen, F., and Zolles, T.: GrSMBMIP: intercomparison of the modelled 1980–2012 surface mass balance over the Greenland Ice Sheet, The Cryosphere 14, 3935–3958,, 2020. 

Fox-Kemper, B., Hewitt, H. T., Xiao, C., Adalgeirsdottir, G., Drijfhout, S. S., Edwards, T. L., Golledge, N. R., Hemer, M., Kopp, R. E., Krinner, G., Mix, A., Notz, D., Nowicki, S., Nurhati, I. S., Ruiz, L., Sallée, J.-B., Slangen, A. B. A., and Yu, Y.: Ocean, Cryosphere and Sea Level Change, in: Climate Change 2021: The Physical Science Basis, Contribution of Working Group 1 to the Sixth Assessment Report of the Intergovernmental Panel on Climate change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthes, J. B. R., Maycock, T. K., Waterfield, T., Yelekci, O., Yu, R., and Zhou, B., Cambridge University Press, Cambridge, UK and New York, NY, USA, 1211–1362, (last access: 2 February 2023), 2021. 

Fyke, J. G., Sacks, W. J., and Lipscomb, W. H.: A technique for generating consistent ice sheet initial conditions for coupled ice sheet/climate models, Geosci. Model Dev., 7, 1183–1195,, 2014. 

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. 

Goldberg, D. N.: A variationally derived, depth-integrated approximation to a higher-order glaciological flow model, J. Glaciol., 57, 157–170, 2011. 

Gomez, N., Gregoire, L., Mitrovica, J., and Payne, A.: Laurentide-Cordilleran Ice Sheet saddle collapse as a contribution to meltwater pulse 1A, Geophys. Res. Lett., 42, 3954–3962,, 2015. 

Gowan, E. J., Zhang, X., Khosravi, S., Rovere, A., Stocchi, P., Hughes, A. L. C., Gyllencreutz, R., Mangerud, J., Svendsen, J.-I., and Lohmann, G.: A new global ice sheet reconstruction for the past 80 000 years, Nat. Commun., 12, 1199,, 2021. 

Huybrechts, P. and de Wolde, J.: The dynamic response of the Greenland and Antarctic ice sheets to multiple-century climatic warming, J. Climate, 1, 2169–2188, 1999. 

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. 

Jungclaus, J., Giorgetta, M., Reick, C., Legutke, S., Brovkin, V., Crueger, T., Esch, M., Fieg, K., Fischer, N., Glushak, K., Gayler, V., Haak, H., Hollweg, H.-D., Kinne, S., Kornblueh, L., Matei, D., Mauritsen, T., Mikolajewicz, U., Müller, W., Notz, D., Pohlmann, T., Raddatz, T., Rast, S., Roeckner, E., Salzmann, M., Schmidt, H., Schnur, R., Segschneider, J., Six, K., Stockhause, M., Wegner, J., Widmann, H., Wieners, K.-H., Claussen, M., Marotzke, J., and Stevens, B.: CMIP5 simulations of the Max Planck Institute for Meteorology (MPI-M) based on the MPI-ESM-P model: The lgm experiment, served by ESGF, WDCC at DKRZ,, 2012. 

Kageyama, M., Harrison, S. P., Kapsch, M.-L., Lofverstrom, M., Lora, J. M., Mikolajewicz, U., Sherriff-Tadano, S., Vadsaria, T., Abe-Ouchi, A., Bouttes, N., Chandan, D., Gregoire, L. J., Ivanovic, R. F., Izumi, K., LeGrande, A. N., Lhardy, F., Lohmann, G., Morozova, P. A., Ohgaito, R., Paul, A., Peltier, W. R., Poulsen, C. J., Quiquet, A., Roche, D. M., Shi, X., Tierney, J. E., Valdes, P. J., Volodin, E., and Zhu, J.: The PMIP4 Last Glacial Maximum experiments: preliminary results and comparison with the PMIP3 simulations, Clim. Past, 17, 1065–1089,, 2021. 

Kleman, J., Fastook, J., Ebert, K., Nilsson, J., and Caballero, R.: Pre-LGM Northern Hemisphere ice sheet topography, Clim. Past, 9, 2365–2378,, 2013. 

Ladant, J.-B., Donnadieu, Y., Lefebvre, V., and Dumas, C.: The respective role of atmospheric carbon dioxide and orbital parameters on ice sheet evolution at the Eocene-Oligocene transition, Paleoceanography, 29, 810–823,, 2014. 

Lambeck, K., Purcell, A., Zhao, J., and Svensson, N.-O.: The Scandinavian Ice Sheet: from MIS 4 to the end of the Last Glacial Maximum, Boreas, 39, 410–435,, 2010. 

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. 

Leguy, G. R., Lipscomb, W. H., and Asay-Davis, X. S.: Marine ice sheet experiments with the Community Ice Sheet Model, The Cryosphere, 15, 3229–3253,, 2021. 

Le Meur, E. and Huybrechts, P.: A comparison of different ways of dealing with isostasy: examples from modelling the Antarctic ice sheet during the last glacial cycle, Ann. Glaciol., 23, 309–317, 1996. 

Liakka, J., Löfverström, M., and Colleoni, F.: The impact of the North American glacial topography on the evolution of the Eurasian ice sheet over the last glacial cycle, Clim. Past, 12, 1225–1241,, 2016. 

Löfverström, M. and Lora, J. M., Abrupt regime shifts in the North Atlantic atmospheric circulation over the last deglaciation, Geophys. Res. Lett., 44, 8047–8055,, 2017. 

Löfverström, M., Caballero, R., Nilsson, J., and Kleman, J.: Evolution of the large-scale atmospheric circulation in response to changing ice sheets over the last glacial cycle, Clim. Past, 10, 1453–1471,, 2014. 

Löfverström, M., Caballero, R., Nilsson, J., and Messori, G.: Stationary Wave Reflection as a Mechanism for Zonalizing the Atlantic Winter Jet at the LGM, J. Atmos. Sci., 73, 3329–3342,, 2016. 

Löfverström, M., Thompson, D. M., Otto-Bliesner, B. L., and Brady, E. C.: The importance of Canadian Arctic Archipelago gateways for glacial expansion in Scandinavia, Nat. Geosci., 15, 482-488,, 2022. 

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. 

Niu, L., Lohmann, G., Hinck, S., Gowan, E. J., and Krebs-Kanzow, U.: The sensitivity of Northern Hemisphere ice sheets to atmospheric forcing during the last glacial cycle using PMIP3 models, J. Glaciol., 65, 645–661,, 2019. 

Ohmura, A., Calanca, P., Wild, M., and Anklin M.: Precipitation, accumulation and mass balance of the Greenland Ice sheet, Z. Gletscherkd. Glazialgeol., 35, 1–20, 1999. 

Otto-Bliesner, B. L. and Brady, E. C.: The sensitivity of the climate response to the magnitude and location of freshwater forcing: last glacial maximum experiments, Quaternary Sci. Rev., 29, 56–73,, 2010. 

Pausata, F. S. R., Li, C., Wettstein, J. J., Kageyama, M., and Nisancioglu, K. H.: The key role of topography in altering North Atlantic atmospheric circulation during the last glacial period, Clim. Past, 7, 1089–1101,, 2011. 

Peltier, W. R., Argus, D. F., and Drummond, R.: Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C(VM5a) model, J. Geophys. Res.-Solid, 120, 450–487,, 2015. 

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

Robinson, A., Goldberg, D., and Lipscomb, W. H.: A comparison of the stability and performance of depth-integrated ice-dynamics solvers, The Cryosphere, 16, 689–709,, 2022. 

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. 

Scherrenberg, M. D. W.: IMAU-paleo/IMAU-ICE, GitHub [code],, last access: 2 February 2023. 

Scherrenberg, M. D. W., Berends, C. J., Stap, L. B., and van de Wal, R. S. W.: Scherrenberg et al., 2023 IMAU-ICE version 2.0 model code and scripts for making figures, Zenodo [code],, 2022a. 

Scherrenberg, M. D. W., Berends, C. J., van de Wal, R. S. W., and Stap, L. B.: Scherrenberg et al., 2023 last glacial cycle ice sheet model output [Data set], Zenodo [data set],, 2022b. 

Simms, A. R., Lisiecki, L., Gebbie, G., Whitehouse, P. L., and Clark, J. F.: Balancing the last glacial maximum (LGM) sea-level budget, Quaternary Sci. Rev., 205, 143–153,, 2019. 

Smith, R. S. and Gregory, J.: The last glacial cycle: transient simulations with an AOGCM, Clim. Dynam., 38, 1545–1559,, 2012. 

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., Berends, C. J., Scherrenberg, M. D., Van De Wal, R. S., and Gasson, E. G.: Net effect of ice-sheet–atmosphere interactions reduces simulated transient Miocene Antarctic ice-sheet variability, The Cryosphere, 16, 1315–1332,, 2022. 

Sueyoshi, T., Ohgaito, R., Yamamoto, A., Chikamoto, M. O., Hajima, T., Okajima, H., Yoshimori, M., Abe, M., O'ishi, R., Saito, F., Watanabe, S., Kawamiya, M., and Abe-Ouchi, A.: Set-up of the PMIP3 paleoclimate experiments conducted using an Earth system model, MIROC-ESM, Geosci. Model Dev., 6, 819–836,, 2013. 

Tarasov, L., Dyke, A. S., Neal, R. M., and Peltier, W. R.: A data-calibrated distribution of deglacial chronologies for the North American ice complex from glaciological modeling, Earth Planet. Sc. Lett., 315, 30–40,, 2012. 

Tierney, J. E., Zhu, J., King, J., Malevich, S. B., Hakim, G. J., and Poulsen, C. J.: Glacial cooling and climate sensitivity revisited, Nature, 584, 569–573,, 2020. 

Toucanne, S., Soulet, G., Riveiros, N. V., Boswell, S. M., Dennielou, B., Waelbroeck, C., Bayon, G., Mojtahid, M., Bosq, M., Sabine, M., Zaragosi, S., Bourillet, J., and Mercier, H.: The North Atlantic Glacial Eastern Boundary Current as a key driver for ice-sheet – AMOC interactions and climate instability, Paleoceanogr. Paleocl., 36, e2020PA004068., 2021. 

Ullman, D. J., LeGrande, A. N., Carlson, A. E., Anslow, F. S., and Licciardi, J. M.: Assessing the impact of Laurentide Ice Sheet topography on glacial climate, Clim. Past, 10, 487–507,, 2014. 

Uppala, S. M., Kållberg, 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., Hólm, 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. Meteorol. Soc., 131, 2961–3012, 2005. 

Voldoire, A., Sanchez-Gomez, E., Salas y Mélia, D., Decharme, B., Cassou, C., Sénési, S., Valcke, S., Beau, I., Alias, A. Chevallier, M., Déqué, M., Deshayes, J., Douville, H., Fernandez, E., Madec, G., Maisonnave, E., Moine, M.-P., Planton, S., Saint-Martin, D., Szopa, S., Tyteca, S., Alkama, R., Belamari, S., Braun, A., Coquart, L., and Chauvin, F.: The CNRM-CM5.1 global climate model: description and basic evaluation, Clim. Dynam., 40, 2091–2121,, 2013. 

Yukimoto, S., Adachi, Y., Hosaka, M., Sakami, T., Yoshimura, H., Hirabara, M., Tanaka, T. Y., Shindo, E., Tsujino, H., Deushi, M., Mizuta, R., Yabu, S., Obata, A., Nakano, H., Ose, T., and Kitoh, A.: A new global climate model of Meteorological Research Institute: MRI-CGCM3 – Model description and basic performance, J. Meteorol. Soc. Jpn., 90a, 23–64, 2012.  

Zheng, W. and Yu, Y.: Paleoclimate simulations of the mid-Holocene and Last Glacial Maximum by FGOALS, Adv. Atmos. Sci., 30, 684–698, 2013. 

Short summary
Ice sheets have a large effect on climate and vice versa. Here we use an ice sheet computer model to simulate the last glacial cycle and compare two methods, one that implicitly includes these feedbacks and one that does not. We found that when including simple climate feedbacks, the North American ice sheet develops from two domes instead of many small domes. Each ice sheet melts slower when including feedbacks. We attribute this difference mostly to air temperature–ice sheet interactions.