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

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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂. Therefore, as CO₂ 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 CO₂ 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 CO₂. 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 1 The 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 1 The 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.

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 CO₂ obtained from Bereiter et al. (2015). Figure 2 depicts the atmospheric CO₂ 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 2 Atmospheric CO₂ 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 CO₂ 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 3 Sea 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 4 The 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 5 Summer (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 CO₂ 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 CO₂ and albedo. Therefore, the same change in CO₂ 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 6 Ice 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 7 LGM 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 8 Ice 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 CO₂. 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 CO₂ 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 9 Sea 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 CO₂, 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 CO₂ 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 CO₂, with timestamps indicated in the figure. Considering an arbitrary threshold of 5 mm yr⁻¹ shows that the glacial-index-method retreat rate accelerates faster at lower CO₂ concentrations (18 ka, 225 ppm) compared to the climate matrix method (16 ka, 240 ppm). This indicates a lower CO₂ 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⁻¹ (11 ka) compared to 19 mm yr⁻¹ (8 ka) with the climate matrix method. This is because when the CO₂ 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 10 Sea level contribution plotted against CO₂ 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 CO₂ is the main driver of temperature changes. However, the CO₂ 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:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\mathrm{Melt}(x,y,m)& =\left(c_{\mathrm{1}}\left(T(x,y,m)-T_{\mathrm{0}}\right)\right.\\ & \left.+c_{\mathrm{2}}(\mathrm{1}-\mathit{\alpha }(x,y,m\left)\right)Q_{\mathrm{TOA}}(x,y,m)-c_{\mathrm{3}}\right)\\ & \left(\mathrm{s}{\mathrm{yr}}^{-\mathrm{1}}\right)/\left(L_{\mathrm{fusion}}\mathrm{12}\mathrm{000}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{surface}}(x,y,m)& ={\mathit{\alpha }}_{\mathrm{snow}}-\left({\mathit{\alpha }}_{\mathrm{snow}}-{\mathit{\alpha }}_{\mathrm{background}}\right)\\ & e^{-\mathrm{15}\mathrm{FirnDepth}(x,y,m-\mathrm{1})}\\ & -\mathrm{0.015}\mathrm{MeltPreviousYr}(x,y).\end{array}\end{array}$$

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).

$$\begin{array}{}\text{(A3)}& \mathrm{SnowFraction}=\mathrm{0.5}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{atan}\frac{\left(T(x,y,m)-T_{\mathrm{0}}\right)}{\mathrm{3.5}}}{\mathrm{1.25664}}}\right)\end{array}$$

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):

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}& \mathrm{SuperImposedWater}(x,y,m)=\\ & \max \left\{\mathrm{0},\mathrm{0.012}\left(T_{\mathrm{0}}-T(x,y,m)\right)\right\},\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}\mathrm{LiquidWater}(x,y,m)& =\mathrm{RainFall}(x,y,m)\\ & +\mathrm{Melt}(x,y,m),\end{array}\text{(A6)}& {\displaystyle }\begin{array}{rl}& \mathrm{Refreezing}(x,y,m)=\min \mathit{\left\{}\mathrm{SuperImposedWater}\right(x,y,m),\\ & \mathrm{LiquidWater}(x,y,m)\mathit{\left\}}P\right(x,y,m).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}\mathrm{SMB}(x,y,m)& =\mathrm{Snowfall}(x,y,m)\\ & +\mathrm{Refreezing}(x,y,m)-\mathrm{Melt}(x,y,m).\end{array}\end{array}$$

Table A1 The 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).

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 B1 LGM 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 B1 Shown 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.

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.

$$\begin{array}{}\text{(C1)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{obs},\mathrm{SL}}(x,y,m)& =T_{\mathrm{obs},\mathrm{PD}}(x,y,m)\\ & +H_{\mathrm{s},\mathrm{obs},\mathrm{PD}}(x,y)\mathit{\lambda }(x,y),\end{array}\text{(C2)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{GCM},\mathrm{SL}}(x,y,m)& =T_{\mathrm{GCM},\mathrm{PD}}(x,y,m)\\ & +H_{\mathrm{s},\mathrm{GCM},\mathrm{PD}}(x,y)\mathit{\lambda }(x,y)\end{array}\end{array}$$

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

$$\begin{array}{}\text{(C3)}& T_{\mathrm{GCM},\mathrm{bias}}(x,y,m)=T_{\mathrm{GCM},\mathrm{SL}}(x,y,m)-T_{\mathrm{obs},\mathrm{SL}}(x,y,m)\end{array}$$

T_GCM,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.

$$\begin{array}{}\text{(C4)}& T_{\mathrm{GCM},\mathrm{corr}}(x,y,m)=T_{\mathrm{GCM}}(x,y,m)-T_{\mathrm{GCM},\mathrm{bias}}(x,y,m)\end{array}$$

Here, T_GCM is the GCM output, either the PI or LGM, while T_GCM,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.

$$\begin{array}{}\text{(C5)}& P_{\mathrm{GCM},\mathrm{bias}}(x,y,m)=P_{\mathrm{GCM},\mathrm{PI}}(x,y,m)/P_{\mathrm{obs},\mathrm{PI}}(x,y,m)\end{array}$$

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

$$\begin{array}{}\text{(C6)}& P_{\mathrm{GCM},\mathrm{corr}}(x,y,m)=P_{\mathrm{GCM}}(x,y,m)/P_{\mathrm{GCM},\mathrm{bias}}(x,y,m)\end{array}$$

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:

$$\begin{array}{}\text{(D1)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{ref}}(x,y,m)& =w_{\mathrm{tot}}(x,y)T_{\mathrm{PI},\mathrm{corr}}(x,y,m)\\ & +\left(\mathrm{1}-w_{\mathrm{tot}}(x,y)\right)T_{\mathrm{LGM},\mathrm{corr}}(x,y,m)\end{array}\text{(D2)}& {\displaystyle }\begin{array}{rl}{P}_{\mathrm{ref}}& =e\left(w_{\mathrm{PI}}(x,y)\log \left(P_{\mathrm{PI},\mathrm{corr}}(x,y,m)\right)\right.\\ & \left.+w_{\mathrm{LGM}}(x,y)\log \left(P_{\mathrm{LGM},\mathrm{corr}}(x,y,m)\right)\right).\end{array}\end{array}$$

Here T_ref and P_ref are the temperature and precipitation climate forcing respectively. The weights for the interpolation are represented by w_tot, w_LGM, and w_PI. The calculation for these interpolation weights determines the difference between the climate matrix and glacial index method.

With the glacial index method w_tot, w_LGM, and w_PI are calculated with only prescribed CO₂.

$$\begin{array}{}\text{(D3)}& w_{{\mathrm{CO}}_{\mathrm{2}}}=\left({\mathrm{CO}}_{\mathrm{2}}-{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{LGM}}\right)/\left({{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{PI}}-{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{LGM}}\right)\end{array}$$

Here, CO₂ represents the prescribed CO₂ concentration obtained from Bereiter et al. (2015) at the current model time step. LGM and PI CO₂ concentrations are 190 and 280 ppm respectively. Note here that $w_{{\mathrm{CO}}_{\mathrm{2}}}$ is a scalar. With the glacial index method, $w_{{\mathrm{CO}}_{\mathrm{2}}}$ is equal to w_tot, w_PI, and (1−w_LGM). Since CO₂ 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, w_tot is calculated using both prescribed CO₂ and absorbed insolation. The annual absorbed insolation is calculated as follows:

$$\begin{array}{}\text{(D4)}& I_{\mathrm{abs}}(x,y)=\sum _{m=\mathrm{1}}^{\mathrm{12}}{Q}_{\mathrm{TOA}}(x,y,m)\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{surface}}(x,y,m)\right).\end{array}$$

Here Q_TOA 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 I_abs, we use the sum of the monthly I_abs. We can use this to calculate the interpolation weights for the absorbed insolation.

$$\begin{array}{}\text{(D5)}& \begin{array}{rl}{w}_{\mathrm{ins}}(x,y)=& \left(I_{\mathrm{abs}}(x,y)-I_{\mathrm{abs},\mathrm{LGM}}(x,y)\right)/\\ & \left(I_{\mathrm{abs},\mathrm{PI}}\left(xy\right)-I_{\mathrm{abs},\mathrm{LGM}}\left(xy\right)\right)\end{array}\end{array}$$

I_abs,PI and I_abs,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 Q_TOA 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, w_ins is not uniform. To obtain the regional effect of albedo on temperature, we calculate the mean w_ins (w_ins,av) and w_ins with a Gaussian smoothing of 200 km (w_ins,smooth). Using these insolation weights, we can calculate the interpolation weight (w_ice) resulting from ice interactions.

In the North America and Eurasia domain,

$$\begin{array}{}\text{(D6)}& {\displaystyle }\begin{array}{rl}& w_{\mathrm{ice}}(x,y)=\\ & {\displaystyle \frac{w_{\mathrm{ins}}(x,y)+\mathrm{3}{w}_{\mathrm{ins},\mathrm{smooth}}(x,y)+\mathrm{3}{w}_{\mathrm{ins},\mathrm{av}}(x,y)}{\mathrm{7}}},\end{array}\text{(D7)}& {\displaystyle }{w}_{\mathrm{tot},\mathrm{T}}(x,y)={\displaystyle \frac{\left(w_{{\mathrm{CO}}_{\mathrm{2}}}(x,y)+w_{\mathrm{ice}}(x,y)\right)}{\mathrm{2}}}.\end{array}$$

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.

$$\begin{array}{}\text{(D8)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{ice}}(x,y)={\displaystyle \frac{w_{\mathrm{ins},\mathrm{smooth}}(x,y)+\mathrm{6}{w}_{\mathrm{ins},\mathrm{av}}}{\mathrm{7}}},\text{(D9)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{tot},\mathrm{T}}(x,y)={\displaystyle \frac{\left(\mathrm{3}{w}_{{\mathrm{CO}}_{\mathrm{2}}}(x,y)+w_{\mathrm{ice}}(x,y)\right)}{\mathrm{4}}}\end{array}$$

The field w_tot 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 (H_s) between the LGM and PI.

$$\begin{array}{}\text{(D10)}& \begin{array}{rl}{w}_{\mathrm{tot},\mathrm{P}}=& \left(\sum H_{\mathrm{s}}(x,y)-\sum H_{\mathrm{s},\mathrm{PI}}(x,y)\right)/\\ & \left(\sum H_{\mathrm{s},\mathrm{LGM}}(x,y)-\sum H_{\mathrm{s},\mathrm{PI}}(x,y)\right)\end{array}\end{array}$$

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

$$\begin{array}{}\text{(D11)}& \begin{array}{rl}{w}_{\mathrm{LGM},H_{\mathrm{s}}}(x,y)=& {\displaystyle \frac{H_{\mathrm{s},\mathrm{ISM}}(x,y)-H_{\mathrm{s},\mathrm{GCM},\mathrm{PI}}(x,y)}{H_{\mathrm{s},\mathrm{GCM},\mathrm{LGM}}(x,y)-H_{\mathrm{s},\mathrm{GCM},\mathrm{PI}}(x,y)}}\\ & w_{\mathrm{tot},\mathrm{P}}(x,y)\end{array}\end{array}$$

In regions that are outside the bounds of this ice extent, $w_{\mathrm{LGM},H_{\mathrm{s}}}$ is equal to w_tot,P. We multiply w_LGM by w_tot to account for the regional effect on precipitation due to a change in surface topography.

$$\begin{array}{}\text{(D12)}& w_{\mathrm{LGM}}(x,y)=w_{\mathrm{LGM},H_{\mathrm{s}}}(x,y)w_{\mathrm{tot},\mathrm{P}}(x,y)\end{array}$$

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