# Late Pleistocene glacial terminations accelerated by proglacial lakes

Meike D W Scherrenberg, Constantijn J Berends, Roderik S W van de Wal

► To cite this version:

Meike D W Scherrenberg, Constantijn J Berends, Roderik S W van de Wal. Late Pleistocene glacial terminations accelerated by proglacial lakes. Climate of the Past, 2024, 20 (8), pp.1761-1784. 10.5194/cp-20-1761-2024. hal-04967376

HAL Id: hal-04967376

https://hal.science/hal-04967376v1

Submitted on 26 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Late Pleistocene glacial terminations accelerated by proglacial lakes

Meike D. W. Scherrenberg<fnref n="1" />, Constantijn J. Berends<fnref n="1" />, and Roderik S. W. van de Wal<fnref n="1" />,<fnref n="2" />

**Correspondence:** Meike D. W. Scherrenberg (m.d.w.scherrenberg@uu.nl)

Received: 9 June 2023 – Discussion started: 3 July 2023

Revised: 24 May 2024 – Accepted: 5 June 2024 – Published: 12 August 2024

**Abstract.** During the glacial cycles of the past 800 000 years, Eurasia and North America were periodically covered by large ice sheets, causing up to 100 m of sea-level change. While Late Pleistocene glacial cycles typically lasted 80 000–120 000 years, the termination phases were completed in only 10 000 years. During these glacial terminations, the North American and Eurasian ice sheets retreated, which created large proglacial lakes in front of the ice-sheet margin. Proglacial lakes accelerate deglaciation as they facilitate the formation of ice shelves at the southern margins of the North American and Eurasian ice sheets. These ice shelves are characterized by basal melting, low surface elevations, and negligible friction at the base. Here, we use an ice-sheet model to quantify the (combined) effects of proglacial lakes on Late Pleistocene glacial terminations by examining their interplay with glacial isostatic adjustment (GIA) and basal sliding. We find that proglacial lakes accelerate the deglaciation of ice sheets mainly because there is an absence of basal friction underneath ice shelves. If friction underneath grounded ice is applied to floating ice, full deglaciation is postponed by a few millennia, resulting in more ice remaining during interglacial periods and no extensive ice shelves forming. Additionally, the large uncertainty in melt rates underneath lacustrine ice shelves translates to an uncertainty in the timing of the termination of up to a millennium.

Proglacial lakes are created by depressions in the landscape that remain after an ice sheet has retreated. The depth, size, and timing of proglacial lakes depend on the rate of bedrock rebound. We find that if bedrock rebounds within a few centuries (rather than a few millennia), the mass loss rate of the ice sheet is substantially reduced. This is because

fast bedrock rebound prevents the formation of extensive proglacial lakes. Additionally, a decrease in ice thickness is partly compensated for by faster bedrock rebound, resulting in a higher surface elevation; lower temperatures; and a higher surface mass balance, which delays deglaciation. We find that a very long bedrock relaxation time does not substantially affect terminations, but it may lead to a delayed onset of the next glacial period. This is because inception regions, such as northwestern Canada, remain below sea level throughout the preceding interglacial period.

## 1 Introduction

From palaeoglaciology, we can learn about which processes are important for the evolution of ice sheets. This can improve our understanding of how the Antarctic and Greenland ice sheets might respond to future warming. During the Late Pleistocene (800 000–10 000 years (800–10 kyr) ago), the North American and Eurasian continents were frequently covered by large ice sheets (e.g. Hughes et al., 2015; Batchelor et al., 2019). While a single glacial cycle lasted on average 80–120 kyr, its decay phase only took 10 kyr. The climate experienced global-scale changes during these glacial terminations, and sea levels rose by up to 130 m (Lambeck et al., 2014; Clark and Tarasov, 2014; Simms et al., 2019), primarily due to the mass loss of major ice sheets. As a consequence, the planetary albedo decreased due to the smaller extent of snow, sea ice, and ice sheets (Abe-Ouchi et al., 2013; Stap et al., 2014). Large volumes of carbon stored in the deep ocean were released (e.g. Denton et al., 2010; Menviel, 2019; Hasenfratz et al., 2019; Sigman et al., 2021),

<footnote n="1">¹Institute for Marine and Atmospheric research Utrecht, Utrecht University, Utrecht, the Netherlands</footnote>

<footnote n="2">²Faculty of Geosciences, Department of Physical Geography, Utrecht University, Utrecht, the Netherlands</footnote>
which contributed to an increase in CO₂ concentrations by 80–100 ppm (parts per million) during interglacial periods (Bereiter et al., 2015). These processes, along with changes in insolation (an important driver of glacial cycles (Miankovitch, 1941)), caused global temperatures to increase by roughly 4–5 °C (Annan et al., 2022). While each of these processes enhanced the mass loss of the ice sheets, glaciated regions also influenced the climate and deglaciation. For instance, it has been suggested that Late Pleistocene deglaciation phases only occurred if the ice sheets were large enough (e.g. Abe-Ouchi et al., 2013; Berends et al., 2021b; Parrenin and Paillard, 2003). Additionally, Berends et al. (2021b) suggested that this may have led to three ice regimes: a small ice sheet that melts at every orbital maximum (e.g. Eurasia and Early Pleistocene North America), a medium-sized ice sheet that survives orbital maxima through elevation–temperature and albedo–temperature feedbacks (e.g. interstadial North America), and a large ice sheet which, due to bedrock–ice feedbacks, becomes sensitive to small increases in summer insolation (e.g. North America during glacial terminations). North America may have experienced a change in size regime during the Mid-Pleistocene Transition (MPT; 1.2–0.8 Ma), when the periodicity of glacial cycles shifted from 40 kyr to an average of 100 kyr. This shift may have been caused by the removal of regolith during consecutive glacial cycles (e.g. Clark and Pollard, 1998; Tabor and Poulsen, 2016; Willeit et al., 2019), exposing the underlying bedrock and reducing sliding. This resulted in thicker ice sheets that were less sensitive to insolation maxima and could survive some insolation optima.

Various processes control the interactions between ice sheets and the climate (and vice versa). Regions with ice and snow have a high albedo, which increases the amount of solar irradiance that is reflected and decreases global and local temperatures. An ice–albedo feedback, where albedo decreases due to the retreat of the ice sheets, amplifies temperature increase and ablation. Additionally, the elevation of an ice sheet influences the surface mass balance, inducing a positive surface-mass-balance–height feedback. A decrease in surface elevation increases temperatures, which enhances ablation and causes a further reduction in ice thickness. However, the ice-sheet mass balance also depends on the amount of accumulation. The amount of precipitation may decline with elevation as decreasing temperatures lower vapour pressure and limit available moisture content (i.e. a negative feedback). At the same time, orographic forcing of precipitation can result in windward and leeward effects, influenced by ice-sheet geometry and prevailing winds. Ice-sheet topography can also influence large-scale atmospheric circulation (Beghin et al., 2015; Löfverström et al., 2016), affecting both temperature and accumulation patterns (Pausata et al., 2011; Ullman et al., 2014; Liakka et al., 2016) on a global scale. In addition to surface melting, ice sheets can also lose mass at the base. Specifically, water underneath ice shelves can facilitate melting, which can thin the ice shelves, thereby reduc-

ing buttressing and accelerating ice mass loss. Additionally, ocean- or lake-terminating ice sheets can lose mass at their margins due to calving.

Besides these changes in forcing, dynamical processes within the ice sheet can also influence mass loss rates. Marine ice sheets with grounding lines resting on a retrograde slope may exhibit instability, where a small perturbation can cause a self-sustained advance or retreat (Weertman, 1974; Schoof, 2012). This process is referred to as marine-ice-sheet instability (MISI), and it is thought to be especially important in the current mass loss in the West Antarctic Ice Sheet (Pattyn, 2018), which has substantial parts of its grounding line resting on retrograde slopes. In the past decade, many improvements have been made in capturing MISI in ice-sheet models (e.g. Pattyn et al., 2012, 2013; Schoof, 2007, 2012; Sun et al., 2020).

The North American and Eurasian ice sheets may have experienced the lacustrine equivalent of MISI during deglacial phases, known as proglacial-lake ice-sheet instability (PLISI; Quiquet et al., 2021; Hinck et al., 2022). Proglacial lakes are created by a combination of glacial isostatic adjustment (GIA; Peltier, 1974) and runoff. The large mass of an ice sheet prompts bedrock deformation, creating a depression in the landscape. As the ice sheet starts to retreat, the rebound lags behind, creating an ice-free depression in front of the ice margin. This depression can fill with meltwater, creating a proglacial lake. Evidence for the existence of large proglacial lakes during the last deglaciation has been found in North America (Lake Agassiz; Upham, 1881; Lepper et al., 2013) and Eurasia (Baltic Ice Lake; Patton et al., 2017). Besides PLISI, proglacial lakes facilitate similar mass loss processes to those of marine shelves, though the rates may be different. Lake calving rates are typically at least 1 order of magnitude lower than those of tide water glaciers (e.g. Warren et al., 1995; Warren and Kirkbride, 2003; Benn et al., 2007). Subshelf melting is also different, which is primarily due to the lack of salinity gradients driving subshelf circulation under ocean conditions (Sugiyama et al., 2016).

The interaction between North America and proglacial lakes has been studied using numerical models (e.g. Tarasov and Peltier, 2006; Hinck et al., 2022; Quiquet et al., 2021; Austermann et al., 2022). Hinck et al. (2022) and Quiquet et al. (2021) studied the effect of PLISI on the deglaciation of North America. They showed that proglacial lakes significantly accelerate ice-sheet melting, with PLISI-induced mass loss being accelerated by increased surface melt rates over low-lying lacustrine shelves. Both studies found that the enhanced retreat caused by proglacial lakes is not due to calving or basal melting but rather to PLISI and the negative surface mass balance. This is because ice shelves have a low surface elevation, high temperatures, and strong ablation.

Here, we build on the work of Hinck et al. (2022) and Quiquet et al. (2021) by considering a wider range of processes related to the presence of proglacial lakes during deglaciation. We do this by employing an ice-sheet model that in-
cludes both the North American and Eurasian ice sheets over
consecutive glacial cycles. We use an ice-sheet model forced
by general-circulation-model climate time slices to study the
transient effect of proglacial lakes on glacial terminations
throughout the Late Pleistocene. Rather than focusing on de-
tailed sea-level projections, our main goal is to investigate ice
dynamical processes that may have contributed to the melt-
ing of the North American and Eurasian ice sheets. We focus
on the effects of basal sliding, shelf formation, and subshelf
melting on glacial terminations. In addition, we consider the
effects of different timescales for GIA response on proglacial
lakes and glacial terminations.

## 2 Methods

### 2.1 Ice-sheet model

We simulated the Northern Hemisphere ice sheets us-
ing a vertically integrated ice-sheet model (version 2.0 of
IMAU-ICE; Berends et al., 2022). The hybrid shallow-
ice–shallow-shelf approximation is used to calculate the
flow of ice (Bueler and Brown, 2009). To model GIA, we
used an elastic-lithosphere–relaxing-asthenosphere (ELRA)
model (Le Meur and Huybrechts, 1996). Basal friction
was calculated using a Budd-type sliding law (Bueler and
van Pelt, 2015; see Appendix A). We included a simple basal-
hydrology scheme from Martin et al. (2011), in which pore
water pressure depends on ice thickness and bedrock height
relative to sea level. The till friction angle in North Amer-
ica was determined using the geology map from Gowan et
al. (2019), which was specifically created for ice-sheet model
applications. For the Eurasian ice sheet, we generated a till
friction angle map based on the sediment thickness map from
Laske and Masters (1997). Here, we used till friction angles
of 30° for sediment thicknesses below 100 m and 10° for
thicknesses exceeding this threshold. The 100 m threshold is
high, but due to the coarse resolution of the map from Laske
and Masters (1997), it does not impact the simulations sub-
stantially. Basal friction at the grounding line is treated us-
ing a subgrid friction-scaling scheme (Berends et al., 2022)
and is based on the approach used in the Community Ice
Sheet Model (CISM; Leguy et al., 2021) and the Parallel
Ice Sheet Model (PISM; Feldmann et al., 2014). In this ap-
proach, the basal-friction coefficient for each model grid cell
is multiplied by the grounded fraction of that cell. Therefore,
friction decreases where ice is partially floating (grounding
line) and is negligible where ice is fully floating (shelves).
This grounded fraction is calculated using the approach from
Leguy et al. (2021), in which the thickness above flotation is
bilinearly interpolated. Berends et al. (2022) showed that the
model performs well in the MISMIP and MISMIP+ bench-
mark experiments, resolving the (migrating) grounding line
within a single grid cell across a range of resolutions.

Calving is parameterized using a simple thickness thresh-
old scheme with a threshold thickness of 200 m. To calculate

the subshelf melt, we use a depth-dependent subshelf-melt
parameterization (Martin et al., 2011) in which basal melting
is linearly related to temperature anomalies. Ocean temper-
ature are parameterized using globally uniform ocean tem-
perature changes (de Boer et al., 2013) that do not account
for regional variations in ocean temperatures. We apply the
same subshelf-melt method to oceans and proglacial lakes
unless stated otherwise. This is a simplification as lacustrine
and marine environments have different thermal structures,
salinities and subshelf circulation regimes (Sugiyama et al.,
2016). Lakes and ocean are simulated when the bedrock is
below the modelled sea level. Sea level is calculated while
maintaining a constant ocean area and only accounting for
ice volume above flotation. As the ice sheet grows, the
bedrock subsides due to GIAs. At glacial maxima, large areas
of the ice sheet can thus become grounded below sea level.
During deglaciation, bedrock uplift lags behind changes in
ice load induced by the thinning and melting of ice. If the
bedrock is below sea level, the cell is classified as one of
three different types of surfaces. (1) If the ice is thick enough
to be grounded, the cell is classified as grounded ice. (2) If
the ice is thin enough to float, it is classified an ice shelf.
(3) If there is no ice, the cell is classified as an ocean. There
is no distinction between lakes and oceans in the model, and
processes such as calving and basal melting are treated the
same for both, indicating that these processes may be over-
estimated with respect to lakes. A further simplification is the
constant surface level of the lakes; the consequences of this
are explored in Sect. 3.5.

The North American, Eurasian, and Greenland ice sheets
are simulated in three separate domains. North America and
Eurasia have a spatial resolution of 40 × 40 km, and Green-
land has a spatial resolution of 20 × 20 km. The boundaries
of these domains are shown in Fig. 1. The higher resolution
used for the Greenland ice sheet results in a similar num-
ber of grid cells compared to the other two domains, while
smaller topographic features are also captured. As shown in
Fig. 1, the domains have some overlapping regions. There-
fore, regions that appear in more than one model domain are
allowed to have ice in only one domain – for example, ice
on Ellesmere Island is only simulated in the North Ameri-
can domain and not in the Greenland domain, while ice on
Greenland itself is not simulated in the North American do-
main. We simulate Greenland and North America in separate
domains; however, they are thought to have merged during
glacial periods.

### 2.2 Climate forcing

To calculate the melting and accumulation of ice, our sur-
face mass balance model (see Sect. 2.4) requires informa-
tion on precipitation and temperature as functions of time
and space. To efficiently compute climate forcing, we in-
terpolate between pre-calculated pre-industrial (PI) and Last
Glacial Maximum (LGM; 21 kyr ago) time slices using a ma-
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The extent of the North American (red), Greenland (green), and Eurasian (blue) domains. The present-day coastline is shown using black lines, while the land and ocean during the LGM are shown in grey. The extent of the LGM ice sheets described by Abe-Ouchi et al. (2015) is shown in white.</figcaption></figure>

trix method (Pollard, 2010). This enables us to implicitly include climate-ice-sheet interactions at lower computational costs compared to fully coupled ice-climate set-ups. Details of this method, which is based on Berends et al. (2018) and Scherrenberg et al. (2023), are described in Appendix B.

The main drivers of temperature change are the external forcing and albedo, each contributing 50% to the final temperature interpolation. External forcing is shown in Fig. 2a and is a combination of CO₂ (Fig. 2b; Bereiter et al., 2015) and insolation changes (Fig. 2c; Laskar et al., 2004). To derive an interpolation weight from this forcing, we first determine an index for CO₂, where 0 represents the LGM climate (190 ppm) and 1 represents the PI climate (280 ppm). We then modify this index using a 65° N summer insolation to capture changes in temperature caused by orbital cycles. Therefore, when summer insolation decreases, climate forcing approaches the LGM leading to cooling. The forcing index remains unchanged if the insolation is 440 W m⁻² (see Fig. 2a and c). As a result, for LGM CO₂ concentrations, the forcing index can still be relatively high for strong insolation, and for PI CO₂ levels, the forcing index can be relatively low for weak insolation values. We tuned the effect of insolation on the forcing index to capture the periodicity of glacial cycles over the past 800 000 years. To calculate the albedo feedback, we multiply the modelled surface albedo by the insolation to obtain the annual amount of insolation absorbed by the surface. We then calculate an interpolation weight from the concurrent amount of absorbed insolation and compare it to the fields obtained using LGM and PI climates and masks. The matrix method includes a precipitation-topography feedback as precipitation is inter-

polated with respect to local and domain-wide changes in the topography relative to the LGM and PI topography.

## 2.3 Climate time slices and downscaling

By simulating the last glacial cycle using an ice-sheet model, it has been shown that the LGM extent and volume were strongly dependent on climate forcing (Charbit et al., 2007; Niu et al., 2019; Alder and Hostetler, 2019; Scherrenberg et al., 2023). Not all general-circulation-model (GCM) simulations can be used to model an LGM ice-sheet extent that agrees well with reconstructions (Scherrenberg et al., 2023; Niu et al., 2019). Here, we use the mean of the MIROC (Sueyoshi et al., 2013), IPSL (Dufresne et al., 2013), COSMOS (Budich et al., 2010), and MPI (Jungclaus et al., 2012) members of phase 3 of the Paleoclimate Modelling Intercomparison Project (PMIP3; Braconnot et al., 2011). This ensemble has been shown to yield a good LGM extent in combination with IMAU-ICE (Scherrenberg et al., 2023). To correct biases in the GCM data, we calculate the difference between the PI time slice and the reanalysis from ERA-40 (Uppala et al., 2005). This bias is then applied to both the PI and LGM time slices. As a consequence, the resulting PI time slice may contain some of the anthropogenic warming present in ERA-40.

The topography and spatial resolution differ between the climate forcing and the ice-sheet model. Therefore, some corrections need to be applied before the climate forcing can be used in IMAU-ICE. First, we bilinearly interpolate the climate forcing onto the finer ice-sheet model grid. As the climate forcing has a lower resolution and therefore a smoother topography, some topographic corrections need to be applied to the temperature and precipitation fields. For temperature, we apply a lapse-rate-based correction. For precipitation, we use the Roe and Lindzen (2001) model to capture the orographic forcing of precipitation on the sloping ice margin and the plateau desert effect in the ice-sheet interior. A more detailed description of the bias correction and downscaling methods can be found in Appendix C.

## 2.4 Surface mass balance model

The surface mass balance (SMB) is calculated monthly using the insolation temperature model IMAU-ITM (Berends et al., 2018). For the present-day climate, this model provides an adequate SMB distribution, as shown in the Greenland Surface Mass Balance Model Intercomparison Project (GrSMB-MIP; Fettweis et al., 2020). Using this model, the accumulation of snow is calculated using the large-scale snow-rain partitioning method proposed by Ohmura et al. (1999). Refreezing is calculated following schemes by Huybrechts and de Wolde (1999) and Janssens and Huybrechts (2000). Ablation is calculated based on Bintanja et al. (2002) and depends on temperature, insolation, and albedo. The equations
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Panel (a) shows the forcing index, which combined with an albedo feedback, drives temperature changes in the ice-sheet model. The forcing index depends on the prescribed CO₂ (b; Bereiter et al., 2015) and insolation (c; Laskar et al., 2004). The pathway of the forcing index over a 100 kyr period is shown in red. A forcing index of 0 (1) represents the LGM (PI) temperature contribution from external forcing.</figcaption></figure>

describing this scheme (IMAU-ITM) are discussed in more detail in Appendix D.

## 3 Results

We conducted a “baseline” simulation by simulating North Hemisphere ice sheets from 782 ka to the present day. Figure 3a shows the extent of the ice sheets during the LGM compared to reconstructions of the North American ice sheet by Dalton et al. (2020) and the Eurasian ice sheet by Hughes et al. (2015; orange contours). The modelled LGM extent matches the reconstructions reasonably well, although ice coverage is lacking in the British Isles. A small proglacial lake formed in North America around 14 kyr ago (Fig. 3d). The North American ice sheet retreats faster than reconstructions suggest for the period starting 12 kyr ago (Fig. 3e–h), with full retreat occurring 10 kyr ago rather than 3 to 4 mil- lennia later. Among other reasons (e.g. uncertainty in the climate forcing, atmospheric circulation, and ice-sheet model parameterizations), this discrepancy could be partially due to the absence of feedback between meltwater, the ocean, and the climate. For example, prior to the Younger Dryas (12.9– 11.7 kyr ago), large amounts of freshwater flowed into the North Atlantic (Teller et al., 2002; Tarasov and Peltier, 2005, 2006; Condron and Winsor, 2012), which could have caused a collapse in the Atlantic Meridional Overturning Circula- tion (AMOC), resulting in cooling (McManus et al., 2004; Velay-Vitow et al., 2024) and a stagnation in sea-level rise (Lambeck et al., 2014). At the present day (see Fig. 3i), the North American and Eurasian ice sheets have fully melted, with the exception of some small ice caps in regions that are currently partly glaciated (e.g. the Arctic Archipelago, Sval- bard, and Iceland).

The total sea-level contribution of this baseline simulation is shown in Fig. 4, where it is compared to ice volume

reconstructions by Spratt and Lisiecki (2016) and Grant et al. (2014). Since we only simulate Northern Hemisphere ice sheets, we added 30% to the ice-sheet contribution to account for sea-level changes caused by phenomena other than the Northern Hemisphere ice sheets, such as the Antarctic (~10 m), Greenland, and Patagonian ice sheets (Simms et al., 2019). While Antarctic ice volume has been shown to be correlated with Northern Hemisphere ice volume as sea-level changes prompt the advance and retreat of grounding lines (Gomez et al., 2020), this approach neglects out-of-phase behaviour of ice mass loss between the Northern and Southern hemispheres.

We find that the modelled sea level captures all major melting events except Marine Isotope Stage (MIS) 13 (530– 480 kyr ago), which had relatively low atmospheric CO₂ concentrations (~240–250 ppm) compared to the interglacial periods that succeeded it (Siegenthaler et al., 2005; Bereiter et al., 2015). The modelled inception periods are longer compared to reconstructions, particularly for the warm interglacial MIS 11 (425–375 kyr ago). This results from the bias correction based on observations in which our PI time slice shows some anthropogenic warming. Consequently, modelled ice inception requires relatively low CO₂ concentrations and weak insolation. The ability of the model to capture the overall pattern of glacial terminations and fully melting North American and Eurasian ice sheets allows us to study the importance of ice dynamical processes that may have contributed to the decay of the ice sheets.

Figure 5 shows the ice volume time series for North America (Fig. 5b), Eurasia (Fig. 5c), and the Northern Hemisphere (Fig. 5a). Here, colours denote net melt (red) and net accumulation (blue), and black squares indicate the onset of deglaciation. In Fig. 5d–f, ice volume is plotted against climate forcing, with a glacial climate on the left and an interglacial climate on the right. As the climate becomes colder, the ice sheet tends to have a net positive mass balance and
<figure><img src="image_5.png" /></figure>

**Figure 3.** Ice thickness and bedrock topography from the baseline simulation during the last deglaciation. Reconstructions of the North American ice sheet by Dalton et al. (2020) and reconstructions of the Eurasian ice sheet by Hughes et al. (2015) are indicated by orange contours, while the present-day coastline is indicated by black contours. Please note that w.r.t stands for "with respect to".

<figure><img src="image_6.png" /></figure>

**Figure 4.** The simulated global mean sea level for the North American (dotted line) and Eurasian (dashed line) ice sheets, derived from the baseline experiment. The global mean sea-level change (solid line) represents the volume change from the Northern Hemisphere ice sheet plus an additional 30% to account for other ice sheets (e.g. Patagonian, Himalayan, and Antarctic ice sheets). These results are compared with data from Grant et al. (2014) and Spratt and Lisiecki (2016).
thus grows. However, as the ice sheet becomes larger and extends further towards the warmer south, CO₂ concentrations and insolation need to be lower to maintain net growth. This aligns with observations from Abe-Ouchi et al. (2013) and Parrenin and Paillard (2003), showing that a larger ice sheet is more vulnerable to increases in insolation and that deglaciation tends to take place when the ice volume is large enough.

The Eurasian ice sheet is more sensitive to collapse compared to the North American ice sheet. Occasionally, the Eurasian ice sheet melts completely (e.g. 168 kyr ago and 50 kyr ago), while the North American ice sheet only undergoes partial melting. This suggest that the Eurasian ice sheet is more vulnerable to complete melting during climate optima and that CO₂ concentrations and insolation can facilitate the decay of the Eurasian ice sheet while remaining favourable enough to allow the North American ice sheet to survive the climatic optimum. On average, the Eurasia ice sheet achieves full deglaciation roughly 3 millennia earlier than the North American ice sheet.

This is in line with the findings of Bonelli et al. (2009), Abe-Ouchi et al. (2013), and Tarasov and Peltier (1997), who found that the Eurasian ice sheet requires lower CO₂ concentrations or insolation levels compared to the North American ice sheet to survive climatic optima. The higher sensitivity of the Eurasian ice sheet is also supported by ice reconstructions, such as those by Gowan et al. (2021) and Mangerud et al. (2023), which show that the Eurasian ice sheet lost most of its volume during the MIS 3 interstadial period (60–25 kyr ago), while the North American ice sheet continued to survive until the LGM. This higher sensitivity of the Eurasian ice sheet to melting during a climate optimum may partly be due to its smaller size (the Eurasian sea-level contribution during the LGM is 19.1 m, compared to 80.1 m for the North American sea-level contribution) and higher LGM temperatures in the climate forcing, resulting largely from the thinner Eurasian ice sheet.

## 3.1 Design of the perturbed experiments

To investigate the effect of proglacial lakes and GIA on Late Pleistocene terminations, we carry out a set of experiments that are similar to the baseline experiment but with one process changed at a time. In the baseline set-up, our model reproduces the basic features of glacial terminations throughout the Late Pleistocene. For the sensitivity experiments, we modify the baseline simulation to investigate the effects of subshelf melting, basal friction of grounded ice, subshelf basal friction, and GIAs on the deglaciation of the Northern Hemisphere ice sheets. Each simulation starts at a point in time corresponding to 782 kyr ago, during an interglacial period when the North American and Eurasian continents were mostly ice-free. In the next few sections, we introduce these perturbation experiments to investigate which processes are

important for the decay of the ice sheets. These experiments are described in Sect. 3.2–3.6 and summarized in Table 1.

## 3.2 The effect of basal friction of grounded ice

The till friction angle is prescribed from and based on a geological map (Gowan et al., 2019; North America) and a sediment map (Laske and Masters, 1997; Eurasia). To assess the impact of basal friction on deglaciation, we conducted two sensitivity experiments. The first involves a low-friction simulation, with till friction angles set to 10° across the North American and Eurasian domains, and is a rough representation of continents fully covered by easily deformable sediments. The second involves a high-friction simulation that roughly represents full bedrock coverage, with till friction angles set to 30°. The basal hydrology is applied uniformly across all till friction angle maps and does not distinguish between sediment and bedrock.

Figure 6 shows a time series of the global mean sea-level contribution for the low-friction, baseline, and high-friction experiments. Figure 6b–k zoom in on individual glacial terminations. Basal friction has a substantial impact on deglaciation. As shown in Fig. 6, decreasing the till friction angle always results in an earlier completion of deglaciation. In the low-friction simulation, the peak melt rate in North America is achieved, on average, 4 centuries earlier compared to the baseline experiment, while in the high-friction simulation, the peak mass loss rate is achieved a millennium later (see Figs. S1 and S2 in the Supplement). In the low-friction scenario, the ice sheet fully melts at MIS 13, while in the high-friction scenario, full deglaciation is not achieved during MIS 13 or MIS 17 (712–676 kyr ago).

Increasing friction also leads to increased ice volume at glacial maxima, although ice extent is not impacted as much. During the LGM, the North American global mean sea-level contribution in the high-friction simulation is 93.6 m, 17% larger than the baseline simulation (80.1 m), while the low-friction simulation has a contribution of 75.3 m (6% smaller). Similarly, the Eurasian sea-level contribution during the LGM is 19.1 m in the baseline simulation, 18.8 m in the low-friction simulation, and 23.7 m in the high-friction simulation. Figure 7 shows maps illustrating ice thickness in the baseline, low-friction, and high-friction simulations. Ice area during the LGM (Fig. 7a, f, and k) deviates by less than 1% for North America and 3% for Eurasia between simulations. Additional ice volume therefore mostly results from an increase in thickness (see Fig. 7a, f, and k). A time series comparing ice volume and area throughout the Late Pleistocene is shown in Fig. S3.

These results suggest that basal friction has a large influence on melt rates during climate optima. Lower friction results in thinner ice sheets with gentler slopes, a large ablation area, and a lower SMB. Combined with increased ice velocities transporting more ice towards the ablation area, these processes can explain the increased sensitivity of the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time series of the simulated ice volume in the baseline experiment are shown in panels (a)–(c). Panels (d)–(f) show the ice volume along with the corresponding (glacial to interglacial) climate forcing, based on the external forcing (prescribed CO₂ and insolation forcing). Panels (a) and (c) show the combined ice volume of Eurasia, Greenland, and North America. Colours indicate net melt (red), net accumulation (blue), and the onset of deglaciation (black squares). The onset of deglaciation is only marked if the volume is at least 30 % of the simulation's maximum at the onset and falls below 30 % by the end of the termination.</figcaption></figure>

<table><caption>**Table 1.** Descriptions of the experiments. Each perturbed experiment is similar to the baseline experiment, except for the described features.</caption><thead><tr><td>Experiment</td><td>Description</td><td>Section</td></tr></thead><tbody><tr><td>Low friction</td><td>Till friction angle is set to 10°, representing full sediment coverage</td><td>3.2</td></tr><tr><td>High friction</td><td>Till friction angle is set to 30°, representing full bedrock coverage</td><td>3.2</td></tr><tr><td>Zero BMB</td><td>Basal mass balance is set to 0 everywhere</td><td>3.3</td></tr><tr><td>Rough water</td><td>The basal friction of floating ice is the same as that of land</td><td>3.4</td></tr><tr><td>Lake 100 m</td><td>Increased lake levels in North America by 100 m</td><td>3.5</td></tr><tr><td>Lake 200 m</td><td>Increased lake levels in North America by 200 m</td><td>3.5</td></tr><tr><td>Lake 300 m</td><td>Increased lake levels in North America by 300 m</td><td>3.5</td></tr><tr><td>Lake 400 m</td><td>Increased lake levels in North America by 400 m</td><td>3.5</td></tr><tr><td>Fast GIA</td><td>GIA relaxation time of 300 years</td><td>3.6</td></tr><tr><td>Slow GIA</td><td>GIA relaxation time of 10 000 years</td><td>3.6</td></tr></tbody></table>

ice sheets during deglaciation and interstadial periods. As friction decreases, the ice sheets may become more sensitive to collapse during climate optima, which forms the basis of the regolith hypothesis, which suggests that the MPT may be explained by an increase in basal drag due to sediment being gradually removed during glacial cycles.

To study the robustness of our results, we also conducted the low-friction and high-friction experiments with the set-up used for the rough-water and fast-GIA simulations, discussed in Sect. 3.4 and 3.6 (see Figs. S4 and S5).

### 3.3 The effect of basal melting on glacial cycles

Proglacial lakes facilitate the formation of ice shelves, which can undergo subshelf melting. Subshelf melting is considered an important process in the current mass loss of Antarctica (Pritchard et al., 2012; Shepherd et al., 2018), where it is dominated by temperature and salinity gradients. Lake Agassiz, a freshwater lake created by the melting ice sheet, could therefore have substantially reduced basal-melt rates compared to parameterizations made for ocean shelves. Here, we conduct a sensitivity test to investigate whether subshelf
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Ice volume time series (metres of sea level equivalent (s.l.e)) for the baseline, high-friction, and low-friction scenarios are shown in panel (a). Grey regions in panel (a) indicate the regions of the zoomed-in time series shown in panels (b)–(k).</figcaption></figure>

melting has a significant impact on the retreat of the North American and Eurasian ice sheets.

The zero BMB experiment deviates from the baseline experiment as the subshelf-melt rate in North America and Eurasia is set to 0. Figure 8 shows the ice volume time series calculated in the zero BMB experiment and compares it to the baseline experiment. Figure 8b–k show the terminations in more detail. The zero BMB experiment is similar to the baseline experiment, although the ice sheets during glacial periods are slightly bigger and shelves can extend far into the ocean. During most terminations, removing subshelf melting only has a small effect, delaying full deglaciation by up to a millennium. MIS 13 and MIS 17 are exceptions. In MIS 17, substantially less ice is melted in the zero BMB experiment compared to the baseline experiment, with the lowest Northern Hemisphere sea-level contribution reaching 15.9 m in the zero BMB experiment and 4.0 m in the baseline experiment. Conversely, more ice is lost during MIS 13 in the zero BMB simulation, despite the lack of subshelf melting. This is caused by the increased sensitivity of the ice sheets when the ice volume is higher, as shown in the relationship between ice-sheet decay and climate forcing in Fig. 5. The zero BMB experiment has a substantially larger ice volume during the glacial period preceding MIS 13, and the retreat produces a substantial proglacial lake, which does not happen in the baseline experiment (see Fig. S6).

### 3.4 The effect of basal friction on floating ice

The ice shelves floating on proglacial lakes or seas experience negligible basal friction, resulting in relatively high flow velocities. To study the impact of this lack of friction, we conduct the rough-water experiment. In this experiment, the subgrid friction coefficient is not scaled with the grounded fraction but is instead calculated as if all grid cells were grounded. Therefore, friction beneath shelves is not negligible anymore. This essentially prevents the formation of shelves, meaning that a migration of the grounding line does not cause a change in friction, which prevents PLISI/MISI. While this is a very unrealistic scenario, the grounding line does not migrate far into the ocean to cover the entire ice domain due to strong ablation, as shown in the ice thickness and bedrock topography maps in Fig. 9i–l.

The sea-level contribution in the rough-water experiment, shown in Fig. 8, can be compared to that in the zero BMB and baseline experiments. During the onset of the termination, the rough-water experiment loses mass at roughly the same pace as the baseline experiment. However, once more than half of the ice volume is lost, the mass loss rate in the rough-water experiment slows down with respect to the baseline experiment. This is because proglacial lakes are only created once the ice sheet has already partly retreated. Therefore, while the baseline and rough-water experiments have a similar retreat rate at the onset of the termination, the retreat in the baseline simulation accelerates once the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Ice thickness and bedrock topography for the low-friction (**a**–**e**), baseline (**f**–**i**), and high-friction (**k**–**o**) experiments. Time slices corresponding to 21 000 years (21 kyr) ago (**a**, **f**, **k**), 12 kyr ago (**c**–**m**), and 10 kyr ago (**d**–**n**) are shown. The LGM extent (**a**, **f**, **k**) from Abe-Ouchi et al. (2015) is indicated by grey contours, whereas black contours indicate the present-day coastlines.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Time series of the North American and Eurasian ice sheets during various deglacial periods. The full 800 kyr time series is shown in panel (a). The grey regions in panel (a) correspond to the time series shown in panels (b)–(k).</figcaption></figure>

proglacial lake has formed. Further differences between the rough-water and baseline simulations can be seen in the transect shown in Fig. 10. While the baseline simulation includes large ice shelves in North America, almost the entire ice sheet is grounded in the rough-water simulation. The shelves in the baseline simulation exhibit large ice velocities due to negligible subshelf friction. Consequently, the surfaces of the ice shelves in the baseline simulation are flat and close to sea level and therefore experience high temperatures and a strongly negative SMB. In the rough-water simulation, these shelves are very small. Without extensive shelves, higher elevation and steeper slopes result in a smaller ablation area.

Interglacial ice volume is generally larger in the rough-water simulation compared to the baseline simulation. The Eurasian ice sheet only melts completely during glacial terminations and not during interstadial periods. Additionally, more ice tends to survive interglacial periods in both North America and Eurasia, which can be observed in the Barents–Kara Sea region and the Arctic Archipelago (see Fig. 9l). However, these results are sensitive to basal friction, and an increase in friction can lead to more ice surviving interstadial periods. In the high-friction equivalent of the rough-water experiment, neither the Eurasian nor the North American ice sheet fully deglaciate, with the smallest sea-level contributions corresponding to 4.5 and 11 m respectively (see Fig. S5).

## 3.5 Lake depth

Lakes are assumed to exist in the model wherever (glacial-isostatic-adjusted) bedrock is below the modelled sea level and where the modelled ice is thin enough to float. In the baseline experiment, the water level of the modelled lakes is assumed to be equal to sea level, while in reality, it can be significantly higher. To assess the effect of this simplifying assumption on our results, we simulate the last deglaciation (21 kyr ago to the present day) using a set of four different constant sea levels (relative to present-day levels) applied across the entire ice-sheet model grid (lake 100 m, lake 200 m, lake 300 m, and lake 400 m experiments). Therefore, lakes are simulated when bedrock is below these fixed sea levels, resulting in increased lake levels. As these sea levels are applied to the entire ice-sheet model grid, the experiments shown here overestimate the melt rates compared to a high-resolution lake model which allows for more drainage. Obviously, very high sea levels (e.g. 400 m) create an unrealistic inland sea that can hinder glacial inception. Therefore, we only apply these fixed levels for the last deglaciation and focus exclusively on the North American ice sheet. Despite this, even our most extreme scenario (400 m) does not induce an immediate collapse after the start of the simulation.

Time series of ice volume for the four lake-level experiments and the baseline experiment are shown in Fig. 11a. The evolution of the ice extent is shown in Fig. 11b–f, with colours indicating when each region of North America was
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Bedrock topography and ice thickness in the baseline, rough-water, slow-GIA, and fast-GIA simulations (columns) corresponding to 21, 11, 9, and 0 ka (rows). The reconstruction from Abe-Ouchi et al. (2015) is shown in grey (**a**, **e**, **i**, **m**, **q**), and the present-day coastlines are shown in black.</figcaption></figure>

last covered by ice. Here, we find that the retreat of the North American ice sheet occurs much faster with increasing sea levels, translating to a lead of roughly 1 millennium in terms of ice volume. These results indicate that lake levels significantly affect the rate of deglaciation for the North American ice sheet and that higher lake levels can substantially accelerate deglaciation.

## 3.6 Glacial isostatic adjustment

Proglacial lakes form through the interaction between GIA, ice sheets, and meltwater. Consequently, the GIA relaxation time controls how quickly bedrock fully recovers from changes in ice load. This relaxation time, as well as the thickness of the ice sheet and the retreat rate, controls the size and shape of the proglacial lake. Here, we assess the effect of GIA response time on multiple glacial cycles while mod-

elling PLISI. We compare three simulations: the previously shown baseline simulation (3000 years), the slow-GIA simulation (10 000 years), and the fast-GIA simulation (300-year relaxation time).

Figure 12 shows the ice volume time series for these three simulations, with smaller panels zooming in on individual terminations (Fig. 12b–k). The retreat in the slow-GIA simulation is generally slower and can lag up to a millennium behind the baseline simulation in terms of ice volume. However, some deglaciations (e.g. the last and penultimate deglaciations) exhibit minimal differences compared to the baseline simulation. This can be explained by differences in the bedrock topography at the onset of the deglaciations and the presence of proglacial lakes. The slower subsidence rates in the slow-GIA simulations may result in higher bedrock topographies. This reduces the modelled saturation at the base as saturation increases with a decreasing bedrock topography
<figure><img src="image_12.png" /></figure>

**Figure 10.** Transects (a, c, e) from the baseline (a, b), rough-water (c, d), and fast-GIA simulations (e, f) corresponding to 11 kyr ago. Present-day bedrock is indicated by dashed lines in panels (a), (c), and (e). The 0 km distance represents the northernmost point of the transect, as shown in panels (b), (d), and (f).

<figure><img src="image_13.png" /></figure>

**Figure 11.** (a) Ice volume time series illustrating varying lake levels compared to the baseline simulation. Panels (b)–(f) show the timing of deglaciation. Background colours are coded from dark to light as follows: ocean in the LGM and PI time slices (dark grey), land in the LGM time slice but not in the PI time slice (grey), land in both the LGM and PI time slices (light grey), and present-day ice coverage (white).

relative to the concurrent sea level. Therefore, in the slow-GIA simulation, basal friction may be lower at the onset of some terminations, causing reduced retreat rates. However, if the bedrock topography during a glacial maximum is similar to that in the baseline simulation, the retreat will also be similar (see Fig. 12c, j, and k). An exception is MIS 13, where the North American ice sheet fully melts in the slow-GIA

simulation but not in the baseline simulation. Retreat rates and ice-sheet volume are not sufficient to induce a proglacial lake in the baseline simulation. However, in the slow-GIA simulation, slower bedrock uplift allows for the creation and sustainment of a proglacial lake throughout this termination event, enabling the North American ice sheet to fully melt (see Fig. S6).
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Time series illustrating ice volume with different GIA relaxation times. The time series for the past 800 kyr is shown in panel (a). Panels (b)–(k) each show one termination, indicated by the grey regions in panel (a).</figcaption></figure>

The ice thickness and bedrock topography in the slow-GIA simulations corresponding to 21, 12, 10, and 0 kyr ago are shown in Fig. 9m–p. During the retreat, the proglacial lakes are significantly larger compared to those in the baseline simulation as bedrock uplift is slower. The slow-GIA simulation also has a delayed onset phase, as seen in Fig. 12. Due to the slow bedrock uplift, large parts of North America still remain below sea level millennia after the ice sheet fully receded (see Fig. 9p). As a result, regions such as northeastern Canada, a key location for ice inception, remained below sea level throughout the entire preceding interglacial period.

The fast-GIA simulation exhibits substantially slower melt rates compared to the baseline, slow-GIA, and rough-water simulations, as shown in the transects in Fig. 10 and the ice thickness maps in Fig. 9q-t. The delayed deglaciation is due to two processes. Firstly, as the ice sheet retreats, rapid bedrock rebound quickly eliminates the depression left by the ice sheet, preventing the formation of proglacial lakes (see Fig. 9q-t). Secondly, the reduction in ice thickness is more efficiently compensated for by the bedrock rebound, which reduces the elevation-temperature feedback, thereby reducing surface melt rates and slowing deglaciation. Therefore, while the fast-GIA and rough-water simulations eliminate the effects of proglacial lakes, only the fast-GIA simulation reduces melt rates from the onset of the termination. During many interglacial periods, including the present, an ice dome persists or the termination is skipped. However, this is again

strongly dependent on basal friction: the low-friction equivalent of the fast-GIA experiment induces full melting during most interglacial periods, while the high-friction equivalent never melts completely (see Figs. S4 and S5). The lack of lakes, combined with the SMB-elevation feedback, makes the fast-GIA simulation the one with the slowest deglaciation.

## 4 Discussion

In this study, we investigate the effect of proglacial lakes on the deglaciation of the Eurasian and North American ice sheets throughout the past 800 kyr. In the baseline configuration, the modelled ice volume over time generally agrees well with different sea-level reconstructions, meaning that all major deglaciations throughout the Late Pleistocene are captured. However, shortcomings include the lack of ice coverage over the British Isles and the fact that our simulations tend to produce interglacial periods that are too long compared to reconstructions. This is especially evident for the MIS 11 interglacial period, which exhibits higher global temperatures and sea levels compared to the present day (Hearty et al., 1999; Raymo and Mitrovica, 2012). Warmer-than-pre-industrial temperatures are also not captured as our climate forcing was interpolated from PI and LGM time slices. Adding additional time slices to our matrix method, similar to the approach used by Abe-Ouchi et al. (2013), especially
for warmer conditions than those of the present day, may im-
prove our representation of various ice volume, CO2 concen-
tration, and orbital parameters. Results nevertheless depend
strongly on the quality of the climate forcing (Scherrenberg
et al., 2023). Here, we chose a forcing that results in an LGM
extent that agrees well with reconstructions, rather than in-
cluding a large number of time slices.

The size and shape of proglacial lakes result from the in-
teraction between GIA and ice thickness. Here, we used uni-
form GIA response times, but in reality, GIA varies spatially
(e.g. Forte et al., 2010; van Calcar et al., 2023) and has large
uncertainties. GIA can increase the size of proglacial lakes
compared to scenarios without GIA (e.g. Austermann et al.,
2022), and in this study, we demonstrated that increasing the
GIA response time leads to proglacial lakes being larger. A
faster GIA response decreases the size of proglacial lakes,
and increased bedrock uplift can partly compensate for thick-
ness loss, resulting in reduced melt rates.

Here, we simulated lakes where bedrock is situated below
sea level. This is a simplification as lakes can exist well above
sea level, and, as shown in the lake 100-400 m experiments,
increasing lake levels results in substantially faster deglacia-
tion. Lake models and drainage algorithms incur high com-
putational costs when applied at a fine topographic resolu-
tion that captures all valleys and drainage channels. While
this may not be feasible, computational costs could be re-
duced by asynchronously coupling with a flood-fill algorithm
(Berends and van de Wal, 2016) or using a drainage algo-
rithm (Tarasov and Peltier, 2006) at a reduced spatial resolu-
tion. These methods do not capture every drainage channel;
however, they may give a reasonable estimate of lake levels
and thus may improve the current method.

In this study, we treated lakes as if they were oceans. How-
ever, as lakes have lower salinity levels and are smaller in
size, they can exhibit substantially different types of sub-
shelf melting (Sugiyama et al., 2016) and lower calving rates
(e.g. Warren et al., 1995; Warren and Kirkbride, 2003; Benn
et al., 2007), and the clogging of icebergs can create a po-
tential buttressing effect (Geirsdóttir et al., 2008). These lim-
itations suggest that the effects of basal melting and calving
may be overestimated in the proglacial lakes, and they are
further compounded by uncertainties in the subshelf-melting
and calving schemes used here. For example, applying a
more sophisticated subshelf-melting scheme which simulates
high melt rates near the grounding line (Rignot and Jacobs,
2002) may result in a greater impact from subshelf melting.
Additionally, we have not explored the effect of calving rates
in detail; however, as shown by Cutler et al. (2001), lake calv-
ing can have an effect on the ice flow of the Laurentide ice
sheet. Nevertheless, our results suggest that while subshelf
melting can enhance the melting of the North American and
Eurasian ice sheets, it is overshadowed by large surface melt
rates and PLISI. Due to these limitations, we cannot give an
exact estimate of how much longer a deglaciation phase will
take when subshelf melting is removed; however, our find-

ings serve as an indication that full melting can still take
place without subshelf melting.

Another limitation is the treatment of basal friction. For
basal hydrology, we used a parameterization from Martin
et al. (2011), and the till friction angle was based on ge-
ology and sediment masks. We performed experiments by
decreasing the till friction angle to 10 or 30°, representing
full sediment or full hard-bed coverage. We found that these
changes significantly influence ice volume at glacial max-
ima, with a much smaller effect on ice extent. Deglaciation is
also slower with increasing friction. While our main conclu-
sions are consistent across various basal-friction scenarios,
the timing of deglaciation and the ice volume that persists
through interglacial periods differ substantially. This high-
lights that the volume of the glacial cycle and melt rates are
sensitive to basal friction. Here, we used a static till fric-
tion angle mask. However, friction can change over time as
bedrock is eroded and sediment is transported. Basal friction
could therefore be improved by including a sediment trans-
port model (e.g. Hildes et al., 2004; Melanson et al., 2013)
and a more sophisticated basal-hydrology method (e.g. Hoff-
man and Price, 2014; Flowers, 2015).

The matrix method used to interpolate between the LGM and PI time slices implicitly includes temperature–albedo and precipitation–topography feedbacks. However, ice-sheet–climate interactions can exhibit threshold behaviours which cannot be simulated using our method. These behaviours include the opening and closing of straits, such as the Canadian Arctic Archipelago (Löfverström et al., 2022); the response of Heinrich and Dansgaard–Oeschger events (Claussen et al., 2003); and rapid changes in ocean circulation and sea ice due to the influx of meltwater into the ocean (Otto-Bliesner and Brady, 2010). Since we do not model the ocean, there are no interactions between meltwater, ocean circulation, and the climate. This may partially explain why the baseline simulation shows excessively high retreat rates from 12 kyr ago, coinciding with the Younger Dryas. The Younger Dryas exhibits reduced rates of sea-level rise (e.g. Lambeck et al., 2014) and lower Northern Hemisphere temperatures due to a stagnation of the AMOC (McManus et al., 2004; Velay-Vitow et al., 2024), likely caused by an influx of meltwater into the North Atlantic (Teller et al., 2002; Tarasov and Peltier, 2005, 2006; Condron and Winsor, 2012). This process is not included in our set-up as temperature is only affected by CO2, insolation, elevation and albedo; therefore, it fails to capture the stagnation in melt rates observed during the Younger Dryas.

Including many of these behaviours would require a model
that more explicitly simulates the climate system. GCMs
may be able to simulate these interactions, but simulating
glacial cycles with a reasonably high resolution would re-
quire an excessive number of computational resources. Alter-
natively, ocean–atmosphere circulation models can be used
to simulate individual glacial terminations (Obase et al.,
2021), or intermediate-complexity models (Ganopolski and
Calov, 2011) can simulate multiple glacial cycles with more explicit feedbacks in the climate system; however, these approaches still require higher computational costs than ice-sheet models and more parameterizations than full GCMs.

## 5 Conclusion

We studied the relative importance of different ice dynamical processes for glacial terminations. The onset of terminations is dominated by a decrease in the SMB, which causes the retreat of the ice sheet. We found that the Eurasian ice sheet requires lower CO₂ concentrations and/or insolation compared to the North American ice sheet in order to survive a climate optimum. Once the ice sheets have retreated significantly, proglacial lakes are created at the margin of the ice sheets. Our results show that these proglacial lakes can significantly accelerate the collapse of the North American and Eurasian ice sheets. If certain processes facilitated by the lakes are removed, the North American and Eurasian ice sheets deglaciate at a reduced pace and may remain partially covered in ice during interglacial periods.

The largest impact of proglacial lakes on deglaciation is their facilitation of low friction under floating ice. If the basal friction of shelves is the same as that of grounded ice, removing the effects of PLISI and MISI, the Eurasian ice sheet only fully melts during interglacial periods rather than during interstadial periods. Additionally, more ice may persist in the Northern Hemisphere during interglacial periods. The high ice velocities caused by negligible subshelf friction create large shelves with low surface elevations, leading to high temperatures and surface melt. However, this process is also sensitive to basal friction. Lowering basal friction decreases ice thickness at glacial maxima, and, consequently, surface temperatures and melting increase. Therefore, reduced friction results in faster deglaciation. This shows that modelling sediments and hydrology is important for simulating glacial cycles. Furthermore, we found that subshelf melting is only a secondary effect with respect to the mass loss of the ice sheets. Applying a zero subshelf-melt rate still results in full deglaciation, although it may take an additional millennium for it to be completed.

We also investigated the effect of the GIA response time on consecutive glacial cycles and deglaciations. If the GIA response is slower than in our baseline simulation, the termination can be slower by up to a millennium, and the subsequent inception phase is delayed. Since inception sites are typically the last regions to deglaciate, land may remain below sea level at the onset of the next glacial period if the bedrock rebound is too slow. We find that a GIA response that is substantially faster than the baseline simulation results in slower deglaciation and larger interglacial ice volumes as the North American and Eurasian ice sheets may not fully deglaciate during some interglacial periods. This is because proglacial lakes are not created when bedrock uplift is too

fast. Additionally, surface melt is reduced as bedrock uplift compensates for ice thickness loss more efficiently.

When projecting the future mass loss of the Greenland and Antarctic ice sheets, marine-ice-sheet dynamics and ice-sheet-climate interactions are considered very important. For Antarctica, the marine equivalent of PLISI, MISI, may induce a runaway collapse of parts of the West Antarctic Ice Sheet (e.g. Ritz et al., 2015). However, this is not a perfect analogy as, for example, the high sensitivity of the Eurasian ice sheet to temperature does not directly translate to future melting in West Antarctica (van Aalderen et al., 2024). The Greenland ice sheet is currently mostly land-terminating, but proglacial lakes are created during the retreat – a phenomenon that is already being observed (e.g. Carrivick and Quincey, 2014). These lakes may accelerate the retreat of the Greenland ice sheet (Carrivick et al., 2022). Therefore, our results underline the relevance of these processes in understanding past ice-sheet evolution and highlight the importance of reproducing this evolution to help constrain these processes in the context of current changes in Antarctica and Greenland.

## Appendix A: Basal sliding, friction, and hydrology

In this study, we used a Budd-type sliding law to simulate the sliding at the base of the ice sheets. To calculate the friction underneath the base of the ice sheet, we applied a parameterization based on Martin et al. (2011). Table A1 lists the units and corresponding values of the constants named here.

The pore water pressure ($\psi$) is parameterized based on Martin et al. (2011). The pore water pressure scaling factor ($\lambda$) determines the saturation of the base and depends on the local bedrock height ($b$) and the concurrent sea level (SL). It is expressed as

$$ \lambda = 1 - \frac{b(x, y) - \text{SL} - b_{\psi, \text{min}}}{b_{\psi, \text{max}} - b_{\psi, \text{min}}}, \quad (A1) $$

where $x$ and $y$ represent the horizontal grid indices. The pore water pressure scaling factor is limited between 1 (fully saturated) and 0 (fully unsaturated). Here, $b_{\psi, \text{max}}$ is the bedrock elevation at which the beds become fully unsaturated (1000 m), relative to sea level, while $b_{\psi, \text{max}}$ is the elevation at which the beds become fully saturated (0 m). This pore water pressure scaling factor is used in combination with the overburden pressure to calculate the pore water pressure ($\psi$), which is given by

$$ \psi = P_w \rho g H(x, y) \lambda. \quad (A2) $$

Here, $H$ represents the thickness of the ice and $P_w$ is the pore water pressure scaling fraction. Here, we use a value of 0.99, which yields a reasonable modelled LGM ice-sheet volume and extent in combination with the prescribed till friction map. To determine the till friction angle ($\phi$), we use different data sets for North America and Eurasia. For North
<table><caption>**Table A1.** Constants describing till friction angle and basal hydrology.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>b<sub>ψ,min</sub></td><td>Elevation with fully saturated beds</td><td>m</td><td>0</td></tr><tr><td>b<sub>ψ,max</sub></td><td>Elevation with fully unsaturated beds</td><td>m</td><td>1000</td></tr><tr><td>ρ</td><td>Density of ice</td><td>kg m<sup>-3</sup></td><td>910</td></tr><tr><td>g</td><td>Gravitational acceleration</td><td>m s<sup>-2</sup></td><td>9.81</td></tr><tr><td>P<sub>w</sub></td><td>Pore water pressure scaling</td><td></td><td>0.99</td></tr><tr><td>q</td><td>Exponent in Budd-type sliding law</td><td></td><td>0.3</td></tr><tr><td>u<sub>0</sub></td><td>Threshold velocity in Budd-type sliding law</td><td>m yr<sup>-1</sup></td><td>100</td></tr></tbody></table>

America, we use a geology reconstruction from Gowan et al. (2019), which provides full coverage of our ice domain and indicates where bedrock could be more easily eroded. The map from Gowan et al. (2019) has been specifically created for ice-sheet modelling. For Eurasia, we created a till friction angle map based on sediment thicknesses from Laske and Masters (1997). Here, till friction angles are 10° when sediment thicknesses exceed 100 m, while for other regions, till friction angles of 30° are used. In the high-friction simulations, till friction angles are 30° in North America and Eurasia, while the low-friction simulations have till friction angles of 10°. For Greenland, which is not the focus of our experiments, we always use a till friction angle map with values of 30° for present-day land and 10° for oceans. The pore water pressure and till friction angle can then be used to calculate the till yield stress (τ), expressed as

$$ \tau(x, y) = \tan\left(\frac{\pi}{180}\right) \phi(\rho g H(x, y) - \psi). \quad (A3) $$

To calculate the basal-friction coefficient (β), we use the Budd-type sliding law (Bueler and van Pelt, 2015), which is given by

$$ \beta(xy) = \tau(x, y) \frac{u^{(q-1)}}{u_0^q}. \quad (A4) $$

Here, *u* represents the ice velocity. The sliding term, is then used to calculate the basal sliding. The basal-friction coefficient is multiplied by the grounded fraction, which is calculated based on Leguy et al. (2021). Therefore, regions that are fully floating (grounded fraction of 0) receive negligible friction, the grounding line (grounded fraction between 0 and 1) receives reduced friction, and fully grounded ice receives full friction (grounded fraction of 1). In the rough-water simulation, the basal-friction coefficient is always multiplied by 1, which means friction is applied as if all ice were grounded.

## Appendix B: Climate time slice interpolation

To provide the ice-sheet model with transiently changing forcing using minimal computational resources, we interpolate between pre-calculated LGM and PI climate time slices. To interpolate the time slices, we use a matrix method. Our

approach is based on Berends et al. (2018) and Scherrenberg et al. (2023) and uses different methods for temperature and precipitation.

To calculate the temperature forcing, we use linear interpolation:

$$ T(x, y, \text{mnth}) = w_T(x, y)T_{\text{PI}}(x, y, \text{mnth}) + (1 - w_T(x, y))T_{\text{LGM}}(x, y, \text{mnth}). \quad (B1) $$

Here, *T* is the monthly (mnth) temperature forcing in the ice-sheet model. *T*<sub>PI</sub> and *T*<sub>LGM</sub> are the climate model temperatures for the PI and LGM time slices respectively. Moreover, *w*<sub>T</sub> is the interpolation weight and depends on two processes: external forcing (*w*<sub>e</sub>) and albedo feedback (*w*<sub>a</sub>). To prevent excessive extrapolation, *w*<sub>T</sub> is capped between -0.25 and 1.25. For North America and Eurasia, *w*<sub>T</sub> is calculated as follows (see Berends et al., 2018):

$$ w_T(x, y) = \frac{w_e(x, y) + w_a(x, y)}{2}. \quad (B2) $$

In Greenland, albedo changes are almost exclusively due to changes in the ice-sheet extent. Our model does not include sea ice, and the Greenland domain does not contain extensive tundra areas. Therefore, we apply a smaller albedo contribution, expressed as

$$ w_T(x, y) = \frac{3w_e(x, y) + w_a(x, y)}{4}. \quad (B3) $$

To calculate *w*<sub>e</sub>, we combine the effects of CO<sub>2</sub> concentration (ppm; Bereiter et al., 2015) and insolation at 65° N (Q<sub>65°N</sub>; W m<sup>-2</sup>; Laskar et al., 2004):

$$ w_e = \frac{C - C_{\text{LGM}}}{C_{\text{PI}} - C_{\text{LGM}}} + \frac{Q_{65^{\circ}\text{N}} - 440 \text{ W m}^{-2}}{70 \text{ W m}^{-2}}. \quad (B4) $$

In this equation, *C*<sub>PI</sub> and *C*<sub>LGM</sub> correspond to 190 and 280 ppm respectively. By including summer (June–July–August) insolation at 65° N (Q<sub>65°N</sub>), the climate can become colder or warmer even if the CO<sub>2</sub> concentration is constant (see Fig. 2). We obtained the ratio between CO<sub>2</sub> and insolation by first conducting a preliminary simulation based on the methods from de Boer et al. (2013) and Berends et al. (2021a), where forcing was modified to match the benthic δ<sup>18</sup>O record from Ahn et al. (2017). This approach essentially reproduces the forcing needed to match the benthic δ<sup>18</sup>O record. We then fitted CO<sub>2</sub> and summer insolation to this forcing to obtain Eq. (B4).

To calculate *w*<sub>a</sub>, which represents albedo feedback, we calculate the annual insolation absorbed by the surface. The absorbed insolation (*I*) depends on the monthly internally calculated surface albedo (*α*<sub>s</sub>) and insolation at the top of the atmosphere (*Q*). *I*(*x*, *y*) is given by

$$ I(x, y) = \sum_{m=1}^{12} Q(x, y, \text{mnth})(1 - \alpha_s(x, y, \text{mnth})). \quad (B5) $$
The albedo is calculated in the ice-sheet model using Eq. (D3). To calculate an interpolation weight for the absorbed insolation ($w_i$), we need to calculate reference fields for the LGM and PI time slices. To calculate the albedo for the LGM time slice, we use land and ice masks from the ice-sheet reconstruction by Abe-Ouchi et al. (2015) as these were also used in the climate model simulations. We integrate the SMB model forward through time with a fixed climate and ice-sheet geometry until the firn layer reaches a steady state (typically after ~30 years). We can then use Eq. (B5) to calculate the reference fields for the absorbed insolation. These absorbed-insolation fields can then be used to calculate $w_i$.

$$w_i(x, y) = (I(x, y) - I_{\text{LGM}}(x, y)) / (I_{\text{PI}}(xy) - I_{\text{LGM}}(xy)) \quad (B6)$$

To account for both local and domain-wide changes in albedo and insolation, we use the following equation based on Berends et al. (2018):

$$w_a(x, y) = \frac{w_i(x, y) + 3w_{i, \text{smooth}}(x, y) + 3w_{i, \text{domain}}(x, y)}{7} \quad (B7)$$

Here, $w_{i, domain}$ is the domain-wide-averaged interpolation weight and $w_i$ is the local interpolation weight. Moreover, $w_{i, smooth}$ represents the regional temperature effect and is calculated by applying a 200 km Gaussian smoothing filter to $w_i$. Once again, we use a different method for Greenland due to the lack of tundra regions; this method is expressed as

$$w_a(x, y) = \frac{w_{i, \text{smooth}}(x, y) + 6w_{i, \text{domain}}}{7} \quad (B8)$$

The interpolation weight used to calculate temperature ($w_T$) can now be derived from $w_a$ and $w_e$ using Eqs. (B2) or (B3). This interpolation weight changes depending on albedo, insolation, and CO$_2$.

For precipitation, we apply a different method as precipitation does not change linearly when the climate cools down and the topography changes. We use the following equation to interpolate precipitation from the climate time slices:

$$P = \exp((1 - w_P(x, y))\log(P_{\text{PI}}(x, y, \text{mnth})) + w_P(x, y)\log(P_{\text{LGM}}(x, y, \text{mnth}))). \quad (B9)$$

Here, $w_P$ is the interpolation weight and depends on local and domain-wide topography changes. First, we compare the domain-wide topography in the model to the climate time slices using the following equation:

$$w_{\text{s, domain}} = \frac{\left(\sum s(x, y) - \sum s_{\text{PI}}(x, y)\right)}{\left(\sum s_{\text{LGM}}(x, y) - \sum s_{\text{PI}}(x, y)\right)} \quad (B10)$$

The surface topography is represented by $s$. Moreover, $w_{s, domain}$ represents the interpolation weight from the domain-wide change in the topography. If a grid cell was

covered with ice during the LGM, we also interpolate with respect to local changes in the topography as follows:

$$w_{\text{s, local}}(x, y) = \frac{S(x, y) - S_{\text{PI}}(x, y)}{S_{\text{LGM}}(x, y) - S_{\text{PI}}(x, y)} w_{\text{s}}(x, y). \quad (B11)$$

If a grid cell did not have ice during the LGM, $w_{s, local}$ is equal to $w_{s, domain}$. In the last step, we multiply the local and regional precipitation effects to obtain the interpolation weight for precipitation, given by

$$w_P(x, y) = w_{\text{s, local}}(x, y)w_{\text{s, domain}}(x, y). \quad (B12)$$

The resulting $w_P$ from Eq. (B12) is used in Eq. (B9) to calculate the interpolated precipitation forcing.

## Appendix C: Downscaling and bias correction

To account for differences between the general-circulation-model (GCM) simulations and the observed climate (ERA-40; Uppala et al., 2005), we apply a bias correction to the LGM and PI snapshots.

To account for temperature biases, we first have to scale the temperature to sea level using a lapse rate correction. This is to account for topographical differences between the GCM and ERA-40 data.

$$T_{\text{obs,SL}}(x, y, \text{mnth}) = T_{\text{obs}}(x, y, \text{mnth}) + s_{\text{obs}}(x, y)\lambda(x, y) \quad (C1)$$

$$T_{\text{GCM,SL}}(x, y, \text{mnth}) = T_{\text{GCM,PI}}(x, y, \text{mnth}) + s_{\text{GCM,PI}}(x, y)\lambda(x, y) \quad (C2)$$

Here, $T$ represents the temperature from ERA-40 (obs) and the pre-industrial (PI) time slices of the climate model (GCM). Surface height is denoted by $s$. The temperature lapse rate ($\lambda$) is equal to $0.008 \text{ K m}^{-1}$. Once the temperature is scaled to sea level (SL), we calculate the temperature difference between the climate model and observed climate as follows:

$$T_{\text{GCM,bias}}(x, y, \text{mnth}) = T_{\text{GCM,SL}}(x, y, \text{mnth}) - T_{\text{obs,SL}}(x, y, \text{mnth}). \quad (C3)$$

This bias correction is then subtracted from the PI and LGM snapshots. As a result, the PI snapshot is equal to the ERA-40 data, which contain some anthropogenic warming.

For precipitation ($P$), biases are applied as ratios rather than absolute differences to ensure that the bias-corrected values are always positive. Therefore, we use the ratio between the model and the observed fields:

$$P_{\text{GCM,bias}}(x, y, \text{mnth}) = P_{\text{GCM}}(x, y, \text{mnth}) / P_{\text{obs}}(x, y, \text{mnth}). \quad (C4)$$

This ratio is used to calculate the bias-corrected precipitation for the PI and LGM time slices and is expressed as

$$P(x, y, \text{mnth}) = P_{\text{GCM}}(x, y, \text{mnth}) / P_{\text{GCM,bias}}(x, y, \text{mnth}). \quad (C5)$$
## Appendix D: Surface mass balance model

The surface mass balance (SMB) is calculated using an insolation temperature model (IMAU-ITM; Berends et al., 2018). To calculate the SMB, ice is added through snow and refreezing and is removed through melting. To calculate the accumulation and ablation of ice, the model requires temperature and precipitation fields, obtained from downscaled and bias-corrected GCM output (see Appendix B and C). To calculate the amount of snowfall, we use a temperature-based snow-rain partitioning method with respect to the melting point ($T_0$) given by Ohmura et al. (1999).

$$f = 0.5 \left( 1 - \frac{\operatorname{atan}\left(\frac{T(x, y, \text{mnth}) - T_0}{3.5}\right)}{1.25664} \right) \quad (\text{D1})$$

The snow fraction ($f$) determines the amount of precipitation that falls as snow; the remainder falls as rain. The snow fraction is always limited between 0 (100% rain) and 1 (100% snow). Furthermore, $x$ and $y$ correspond to the horizontal grid, while $m$ indicates the month. To calculate the ablation of ice, we use the parameterized scheme from Bintanja et al. (2002), which accounts for ablation resulting from temperature and insolation:

$$M(x, y, m) = (c_1(T(x, y, \text{mnth}) - T_0) + c_2(1 - \alpha(x, y, \text{mnth}))Q(x, y, \text{mnth}) - c_3). \quad (\text{D2})$$

Here, $T$ is the 2 m air temperature, $T_0$ is the melting temperature of ice (273.16 K), and $Q_{TOA}$ is the insolation at the top of the atmosphere (Laskar et al., 2004). The parameters $c_1$ and $c_2$ are set to $0.079 \text{ m yr}^{-1} \text{ K}^{-1}$ and $7.9 \times 10^{-4} \text{ m J}^{-1} \text{ yr}^{-1}$ respectively. The parameter $c_3$ is used for tuning. Here, we tuned the model to obtain realistic LGM ice volumes, with $c_3$ values provided for North America ($-0.16 \text{ m yr}^{-1}$), Eurasia ($-0.24 \text{ m yr}^{-1}$), and Greenland ($0.19 \text{ m yr}^{-1}$). Albedo ($\alpha$), calculated in the ice-sheet model, is also based on Bintanja et al. (2002) and is given as follows:

$$\alpha(x, y, \text{mnth}) = \alpha_{\text{snow}} - (\alpha_{\text{snow}} - \alpha_{\text{b}}) \exp^{-15D(x, y, \text{mnth}-1)} - 0.015 M_{\text{prev}}(x, y). \quad (\text{D3})$$

Here, $\alpha_{\text{b}}$ represents the albedo without any snow, with values of 0.5 for bare ice, 0.2 for land, and 0.1 for water. The amount of melt (m) from the previous year is defined as $Melt_{\text{prev}}$. If snow is added on top, which increases the depth (m) of the firn layer ($D$), the albedo can increase up to $\alpha_{\text{snow}}$, which represents the albedo of fresh snow (0.85). Therefore, the albedo in the model varies between the background and snow albedo. The depth of the firn layer is calculated using the amount of snow added on top without melting.

Some of the melt and rainfall can refreeze in the model. Here, we use the approach by Huybrechts and de Wolde (1999), employing the total amount (m yr$^{-1}$) of liquid water ($l$), superimposed water ($s$), and precipitation ($P$).

$$s(x, y, \text{mnth}) = \max\{0, 0.012(T_0 - T(x, y, \text{mnth}))\} \quad (\text{D4})$$

$$l(x, y, \text{mnth}) = R(x, y, \text{mnth}) + M(x, y, \text{mnth}) \quad (\text{D5})$$

$$r(x, y, \text{mnth}) = \min\{\min\{S(x, y, \text{mnth}), L(x, y, \text{mnth})\}, P(x, y, \text{mnth})\} \quad (\text{D6})$$

$$P(x, y, \text{mnth}) \quad (\text{D6})$$

By combining snowfall, refreezing, and melt, the SMB can be calculated as follows:

$$\begin{aligned} \text{SMB}(x, y, \text{mnth}) = & S(x, y, \text{mnth}) + r(x, y, \text{mnth}) \\ & - M(x, y, \text{mnth}). \end{aligned} \quad (\text{D7})$$

**Code availability.** IMAU-ICE is described by Berends et al. (2022), and the most recent version of the source code can be found at https://github.com/IMAU-paleo/IMAU-ICE (last access: 31 July 2024). The version used in this study, along with the configuration files, is available on Zenodo (https://doi.org/10.5281/zenodo.11634769, Scherrenberg, 2024). Running the simulations requires additional files for CO$_2$ (see Bereiter et al., 2015), climate (PMIP3 database: https://esgf-node.ipsl.upmc.fr/search/cmip5-ipsl/, last access: 31 July 2024), insolation (Laskar et al., 2004), and initial topography (ETOPO: https://doi.org/10.7289/V5C8276M, Amante and Eakins, 2009; BedMachine: https://doi.org/10.5067/5XKQD5Y5V3VN, NSIDC, 2024). The LGM land and ice masks were obtained from Abe-Ouchi et al. (2015). Till friction angles were obtained from Gowan et al. (2019) and Laske and Masters (1997). For more information, contact the corresponding author.

**Data availability.** The results are available with 5 kyr (2D fields) and 100-year (scalar) output frequencies on Zenodo (https://doi.org/10.5281/zenodo.11634769, Scherrenberg, 2024). Additional 2D fields and higher output frequencies of up to 1 kyr can be requested by contacting the corresponding author.

**Supplement.** The supplement related to this article is available online at: https://doi.org/10.5194/cp-20-1761-2024-supplement.

**Author contributions.** MDWS conducted the simulations wrote the paper. The set-up for the experiments was created by RSWvdW, CJB, and MDWS. CJB provided model support. All authors contributed to the paper and the analysis of the results.

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

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
**Special issue statement.** This article is part of the special issue "Icy landscapes of the past". It is not associated with a conference.

**Acknowledgements.** NWO (Dutch Research Council) Exact and Natural Sciences provided support for the supercomputer facilities of the Dutch national supercomputer Snellius. We would also like to acknowledge the support from SurfSARA Computing and Networking Services.

**Financial support.** Meike D. W. Scherrenberg is supported by the Netherlands Earth System Science Centre (NESSC), which is funded by the Ministry of Education, Culture and Science (OCW; grant no. 024.002.001). Constantijn J. Berends is funded by PROTECT. This project has received funding from the European Union's Horizon 2020 research and innovation programme (grant no. 869304).

**Review statement.** This paper was edited by Lev Tarasov and reviewed by Niall Gandy and Fuyuki Saito.

## References

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, https://doi.org/10.1038/nature12374, 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, https://doi.org/10.5194/gmd-8-3621-2015, 2015.

Ahn, S., Khider, D., Lisiecki, L. E., and Lawrence, C. E.: A probabilistic Pliocene–Pleistocene stack of benthic $δ^{18}$O using a profile hidden Markov model, Dynam. Stat. Clim. Syst., 2, dxz002, https://doi.org/10.1093/climsys/dzx002, 2017.

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, https://doi.org/10.1007/s00382-019-04663-x, 2019.

Amante, C. and Eakins, B. W.: ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24, National Geophysical Data Center, NOAA [data set], https://doi.org/10.7289/V5C8276M, 2009.

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, https://doi.org/10.5194/cp-18-1883-2022, 2022.

Austermann, J., Wickert, A. D., Pico, T., Kingslake, J., Callaghan, K. L., and Creel, R. C.: Glacial isostatic adjustment shapes proglacial lakes over glacial cycles, Geophys. Res. Lett., 50, e2022GL101191, https://doi.org/10.1029/2022GL101191, 2022.

Batchelor, C. L., Margold, M., Krapp, M., Murton, D., 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, 3713, https://doi.org/10.1038/s41467-019-11601-2, 2016.

Beghin, P., Charbit, S., Dumas, C., Kageyama, M., and Ritz, C.: How might the North American ice sheet influence the northwestern Eurasian climate?, Clim. Past, 11, 1467–1490, https://doi.org/10.5194/cp-11-1467-2015, 2015.

Benn, D. I., Warren, C. R., and Mottram, R. H.: Calving processes and the dynamics of calving glaciers, Earth-Sci. Rev., 82, 143–179, https://doi.org/10.1016/J.EARSCIREV.2007.02.002, 2007.

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 CO$_2$ record from 800 to 600 kyr before present, Geophys. Res. Lett., 42, 542–549, https://doi.org/10.1002/2014GL061957, 2015.

Berends, C. J. and van de Wal, R. S. W.: A computationally efficient depression-filling algorithm for digital elevation models, applied to proglacial lake drainage, Geosci. Model Dev., 9, 4451–4460, https://doi.org/10.5194/gmd-9-4451-2016, 2016.

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, https://doi.org/10.5194/gmd-11-4657-2018, 2018.

Berends, C. J., de Boer, B., and van de Wal, R. S. W.: Reconstructing the evolution of ice sheets, sea level, and atmospheric CO$_2$ during the past 3.6 million years, Clim. Past, 17, 361–377, https://doi.org/10.5194/cp-17-361-2021, 2021a.

Berends, C. J., Köhler, P., Lourens, L. J., and van de Wal, R. S. W.: On the cause of the mid-Pleistocene transition, Rev. Geophys., 59, e2020RG000727, https://doi.org/10.1029/2020RG000727, 2021b.

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, https://doi.org/10.5194/gmd-15-5667-2022, 2022 (code available at: https://github.com/IMAU-paleo/IMAU-ICE, last access: 31 July 2024).

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.

Bonelli, S., Charbit, S., Kageyama, M., Woillez, M.-N., Ramestein, G., Dumas, C., and Quiquet, A.: Investigating the evolution of major Northern Hemisphere ice sheets during the last glacial-interglacial cycle, Clim. Past, 5, 329–345, https://doi.org/10.5194/cp-5-329-2009, 2009.

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.

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, https://pure.mpg.de/rest/items/item_2193290/component/file_2193291/content (last access: 2 February 2023), 2010.

Bueler, E. and Brown, J.: The shallow shelf approximation as a sliding law in a thermomechanically cou-
pled ice sheet model, J. Geophys. Res., 114, F03008,
https://doi.org/10.1029/2008JF001179, 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, https://doi.org/10.5194/gmd-8-1613-2015, 2015.

Carrivick, J. L. and Quincey, D. J.: Progressive increase in num-
ber and volume of ice-marginal lakes on the western margin of
the Greenland Ice Sheet, Global Planet. Change, 116, 156–163,
https://doi.org/10.1016/j.gloplacha.2014.02.009, 2014.

Carrivick, J. L., How, P., Lea, J. M., Sutherland, J. L., Grimes, M., Tweed, F. S., Cornford, S., Quincey, D. J., and Mallalieu, J.: Ice-Marginal Proglacial Lakes Across Greenland: Present Status and a Possible Future, Geophys. Res. Lett., 49, e2022GL099276, https://doi.org/10.1029/2022GL099276, 2022.

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, https://doi.org/10.5194/cp-3-15-2007, 2007.

Clark, P. and Pollard, D.: Origin of the Middle Pleistocene Transi-
tion by ice sheet erosion of regolith, Paleoceanography, 13, 1–9,
https://doi.org/10.1029/97PA02660, 1998.

Clark, P. U. and Tarasov, L.: Closing the sea level budget at the Last Glacial Maximum, P. Natl. Acad. Sci., 111, 15861–15862, https://doi.org/10.1073/pnas.1418970111, 2014.

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, https://doi.org/10.1007/s00382-003-0336-2, 2003.

Condron, A. and Winsor, P.: Meltwater routing and the Younger Dryas, P. Natl. Acad. Sci., 109, 19928–19933, https://doi.org/10.1073/pnas.1207381109, 2012.

Cutler, P. M., Mickelson, D. M., Colgan, P. M., MacAyeal, D. R., and Parizek, B. R.: Influence of the Great Lakes on the Dynamics of the Southern Laurentide Ice Sheet: Numerical Experiments, Geology, 29, 1039–1042, https://doi.org/10.1130/0091-7613(2001)029<1039:IOTGLO>2.0.CO;2, 2001.

Dalton, A. S., Margold, M., Stokes, C. R., Tarasov, L., Dyke, A. S., Adams, R. S., Allard, S., Arends, H. E., Atkinson, N., Attig, J. W., Barnett, P. J., Barnett, R. L., Batterson, M., Bernatchez, P., Borns, H. W., Breckenridge, A., Briner, J. P., Brouard, E., Campbell, J. E., Carlson, A. E., Clague, J. J., Curry, B. B., Daigneault, R. A., Dubé-Loubert, H., Easterbrook, D. J., Franzi, D. A., Friedrich, H. G., Funder, S., Gauthier, M. S., Gowan, A. S., Harris, K. L., Hétu, B., Hooyer, T. S., Jennings, C. E., Johnson, M. D., Kehew, A. E., Kelley, S. E., Kerr, D., King, E. L., Kjeldsen, K. K., Knaeble, A. R., Lajeunesse, P., Lakeman, T. R., Lamothe, M., Larson, P., Lavoie, M., Loope, H. M., Lowell, T. V., Lusardi, B. A., Manz, L., McMartin, I., Nixon, F. C., Occhietti, S., Parkhill, M. A., Piper, D. J., Pronk, A. G., Richard, P. J., Ridge, J. C., Ross, M., Roy, M., Seaman, A., Shaw, J., Stea, R. R., Teller, J. T., Thompson, W. B., Thorleifsson, L. H., Utting, D. J., Veillette, J. J., Ward, B. C., Weddle, T. K., and Wright, H. E.: An updated radiocarbon-based ice margin chronology for the last deglaciation of the North American Ice Sheet Complex, Quaternary Sci. Rev., 234, 106223, https://doi.org/10.1016/j.quascirev.2020.106223, 2020.

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, https://doi.org/10.1007/s00382-012-1562-2, 2013.

Denton, G. H., Anderson, R. F., Toggweiler, J., Edwards, R., Schaefer, J., and Putnam, A.: The last glacial termination, Science, 328, 1652–1656, https://doi.org/10.1126/science.1184119, 2010.

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., Lalellec, 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, https://doi.org/10.1007/s00382-012-1636-1, 2013.

Feldmann, J., Albrecht, T., Khroulev, C., Pattyn, F., and Levermann, A.: Resolution-dependent performance of grounding line motion in a shallow model compared with a full-Stokes model according to the MISMIP3d intercomparison, J. Glaciol., 60, 353–359, https://doi.org/10.3189/2014JoG13J093, 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, https://doi.org/10.5194/tc-14-3935-2020, 2020.

Flowers, G. E.: Modelling water flow under glaciers and ice sheets, P. Roy. Soc. A, 471, 20140907, https://doi.org/10.1098/rspa.2014.0907, 2015.

Forte, A., Moucha, R., Simmons, N. A., Grand, S. P., and Mitrovica, J. X.: Deep-mantle contributions to the surface dynamics of the North American continent, Tectonophysics, 481, 3–15, https://doi.org/10.1016/j.tecto.2009.06.010, 2010.

Ganopolski, A. and Calov, R.: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles, Clim. Past, 7, 1415–1425, https://doi.org/10.5194/cp-7-1415-2011, 2011.

Geirsdóttir, Á., Miller, G. H., Wattrus, N. J., Björnsson, H., and Thors, K.: Stabilization of glaciers terminating in closed water bodies: Evidence and broader implications, Geophys. Res. Lett., 35, L17502, https://doi.org/10.1029/2008GL034432, 2008.

Gomez, N., Weber, M. E., Clark, P. U., Mitrovica, J. X., and Han, H. K.: Antarctic ice dynamics amplified by Northern Hemisphere sea-level forcing, Nature, 587, 600–604, https://doi.org/10.1038/s41586-020-2916-2, 2020.
Gowan, E. J., Niu, L., Knorr, G., and Lohmann, G.: Geology datasets in North America, Greenland and surrounding areas for use with ice sheet models, Earth Syst. Sci. Data, 11, 375–391, https://doi.org/10.5194/essd-11-375-2019, 2019.

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, https://doi.org/10.1038/s41467-021-21469-w, 2021.

Grant, K. M., Rohling, E. J., Bronk Ramsey, C., Cheng, H., Edwards, R. L., Florindo, F., Heslop, D., Marra, F., Roberts, A. P., Tamisiea, M. E., and Williams, F.: Sea-level variability over five glacial cycles, Nat. Commun., 5, 5076, https://doi.org/10.1038/ncomms6076, 2014.

Hasenfratz, A. P., Jaccard, S. L., Martínez-García, A., Sigman, D. M., Hodell, D. A., Vance, D., Bernasconi, S. M., Kleiven, H. F., Haumann, F. A., and Haug, G. H.: The residence time of Southern Ocean surface waters and the 100 000-year ice age cycle, Science, 363, 1080–1084, https://doi.org/10.1126/science.aat7067, 2019.

Hearty, P. J., Kindler, P., Cheng, H., and Edwards, R. L.: A +20 m middle Pleistocene sea-level highstand (Bermuda and the Bahamas) due to partial collapse of Antarctic ice, Geology, 27, 375–378, https://doi.org/10.1130/0091-7613(1999)027<0375:AMMPSL>2.3.CO;2, 1999.

Hildes, D. H., Clarke, G. K., Flowers, G. E., and Marshall, S. J.: Subglacial erosion and englacial sediment transport modelled for North American ice sheets, Quaternary Sci. Rev., 23, 409–430, https://doi.org/10.1016/j.quascirev.2003.06.005, 2004.

Hinck, S., Gowan, E. J., Zhang, X., and Lohmann, G.: PISM-LakeCC: Implementing an adaptive proglacial lake boundary in an ice sheet model, The Cryosphere, 16, 941–965, https://doi.org/10.5194/tc-16-941-2022, 2022.

Hoffman, M. J. and Price, S.: Feedbacks between coupled subglacial hydrology and glacier dynamics, J. Geophys. Res.-Earth, 119, 414–436, https://doi.org/10.1002/2013JF002943, 2014.

Hughes, A. L. C., Gyllencreutz, R., Lohne, Ø. S., Mangerud, J., and Svendsen, J. I.: The last Eurasian ice sheets – a chronological database and time-slice reconstruction, DATED-1, Boreas, 45, 1–45, https://doi.org/10.1111/ber.12142, 2015.

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., Stockhouse, 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 Igm experiment, served by ESGF, WDCC at DKRZ, https://doi.org/10.1594/WDCC/CMIP5.MXEPlg, 2012.

Lambeck, K., Rouby, H., Purcell, A., Sun, Y., and Sambridge, M.: Sea level and global ice volumes from the Last Glacial Maximum

to the Holocene, P. Natl. Acad. Sci. USA, 111, 15296–15303, https://doi.org/10.1073/pnas.1411762111, 2014.

Laskar, J., Robutel, P., Joutel, F., 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, https://doi.org/10.1051/0004-6361:20041335, 2004.

Laske, G. and Masters, G.: A Global Digital Map of Sediment Thickness, EOS T. Am. Geophys. Union, https://igppweb.ucsd.edu/~gabi/sediment.html (last access: 30 July 2024), 1997.

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, https://doi.org/10.5194/tc-15-3229-2021, 2021.

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

Lepper, K., Buell, A. W., Fisher, T. G., and Lowell, T. V.: A chronology for glacial Lake Agassiz shorelines along Upham's namesake transect, Quatern. Res., 80, 88–98, https://doi.org/10.1016/j.yqres.2013.02.002, 2013.

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, https://doi.org/10.5194/cp-12-1225-2016, 2016.

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, https://doi.org/10.1175/JAS-D-15-0295.1, 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, https://doi.org/10.1038/s41561-022-00956-9, 2022.

Mangerud, J., Alexanderson, H., Birks, H. H., Paus, A., Perić, Z. M., and Svendsen, J. I.: Did the Eurasian ice sheets melt completely in early Marine Isotope Stage 3? New evidence from Norway and a synthesis for Eurasia, Quaternary Sci. Rev., 317, 108274, https://doi.org/10.1016/j.quascirev.2023.108136, 2023.

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, https://doi.org/10.5194/tc-5-727-2011, 2011.

McManus, J. F., Francois, R., Gherardi, J. M., Keigwin, L. D., and Brown-Leger, S.: Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes, Nature, 428, 834–837, https://doi.org/10.1038/nature02494, 2004.

Melanson, A., Bell, T., and Tarasov, L.: Numerical modelling of subglacial erosion and sediment transport and its application to the North American ice sheets over the Last Glacial cycle, Quaternary Sci. Rev., 68, 154–174, https://doi.org/10.1016/j.quascirev.2013.02.017, 2013.

Menviel, L.: The southern amplifier, Science, 363, 1040–1041, https://doi.org/10.1126/science.aaw7196, 2019.

Milankovitch, M.: Kanon der Erdbestrahlung und Seine Andwendung auf das Eiszeitenproblem, Royal Serbian Academy Special Publication 132, Royal Serbian Academy, Belgrade, Serbia, 1941.
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, https://doi.org/10.1017/jog.2019.42, 2019.

NSIDC: IceBridge BedMachine Greenland, Version 1, NSIDC [data set], https://doi.org/10.5067/5XKQD5Y5V3VN, 2024.

Obase, T., Abe-Ouchi, A., and Saito, F.: Abrupt climate changes in the last two deglaciations simulated with different Northern ice sheet discharge and insolation, Sci. Rep., 11, 22359, https://doi.org/10.1038/s41598-021-01651-2, 2021.

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, https://doi.org/10.1016/j.quascirev.2009.07.004, 2010.

Parrenin, F. and Paillard, D.: Amplitude and phase of glacial cycles from a conceptual model, Earth Planet. Sc. Lett., 214, 243–250, 2003.

Patton, H., Hubbard, A., Andreassen, K., Auriac, A., Whitehouse, P. L., Stroeven, A. P., Shackleton, C., Winsborrow, M., Heyman, J., and Hall, A. M.: Deglaciation of the Eurasian ice sheet complex, Quaternary Sci. Rev., 169, 148–172, https://doi.org/10.1016/j.quascirev.2017.05.019, 2017.

Pattyn, F.: The paradigm shift in Antarctic ice sheet modelling, Nat. Commun., 9, 2728, https://doi.org/10.1038/s41467-018-05003-z, 2018.

Pattyn, F., Schoof, C., Perichon, L., Hindmarsh, R. C. A., Bueler, E., de Fleurian, B., Durand, G., Gagliardini, O., Gladstone, R., Goldberg, D., Gudmundsson, G. H., Huybrechts, P., Lee, V., Nick, F. M., Payne, A. J., Pollard, D., Rybak, O., Saito, F., and Vieli, A.: Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP, The Cryosphere, 6, 573–588, https://doi.org/10.5194/tc-6-573-2012, 2012.

Pattyn, F., Perichon, L., Durand, G., Favier, L., Gagliardini, O., Hindmarsh, R. C. A., Zwinger, T., Albrecht, T., Cornford, S., Docquier, D., Fürst, J. J., Goldberg, D., Gudmundsson, G. H., Humbert, A., Hütten, M., Huybrechts, P., Jouvet, G., Kleiner, T., Larour, E., Martin, D., Morlighem, M., Payne, A. J., Pollard, D., Rückamp, M., Rybak, O., Seroussi, H., Thoma, M., and Wilkens, N.: Grounding-line migration in plan-view marine ice-sheet models: results of the ice2sea MISMIP3d intercomparison, J. Glaciol., 59, 410–422, https://doi.org/10.3189/2013JoG12J129, 2013.

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, https://doi.org/10.5194/cp-7-1089-2011, 2011.

Peltier, W. R.: The impulse response of a Maxwell Earth, Rev. Geophys., 12, 649–669, https://doi.org/10.1029/RG012i004p00649, 1974.

Pollard, D.: A retrospective look at coupled ice sheet–climate modeling, Climatic Change, 100, 173–194, https://doi.org/10.1007/s10584-010-9830-9, 2010.

Pritchard, H. D., Ligtenberg, S. R. M., Fricker, H. A., Vaughan, D. G., Van Den Broeke, M. R., and Padman, L.: Antarctic ice-sheet

loss driven by basal melting of ice shelves, Nature, 484, 502–505, https://doi.org/10.1038/nature10968, 2012.

Quiquet, A., Dumas, C., Paillard, D., Ramstein, G., Ritz, C., and Roche, D. M.: Deglacial Ice Sheet Instabilities Induced by Proglacial Lakes, Geophys. Res. Lett., 48, e2020GL092141, https://doi.org/10.1029/2020GL092141, 2021.

Raymo, M. E. and Mitrovica, J. X.: Collapse of polar ice sheets during the stage 11 interglacial, Nature, 483, 453–456, https://doi.org/10.1038/nature10891, 2012.

Rignot, E. and Jacobs, S.: Rapid Bottom Melting Widespread near Antarctic Ice Sheet Grounding Lines, Science, 296, 2020–2023, https://doi.org/10.1126/science.1070942, 2002.

Ritz, C., Edwards, T. L., Durand, G., Payne, A. J., Peyaud, V., and Hindmarsh, R. C. A.: Potential sea-level rise from Antarctic ice-sheet instability constrained by observations, Nature, 528, 115–118, https://doi.org/10.1038/nature16147, 2015.

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.: Scherrenberg et al. (2024) supplement (Climate of the past): Ice-sheet model code, and output of Northern Hemisphere ice-sheet evolution of the past 800 kyr, Zenodo [data set], https://doi.org/10.5281/zenodo.11634769, 2024.

Scherrenberg, M. D. W., Berends, C. J., Stap, L. B., and van de Wal, R. S. W.: Modelling feedbacks between the Northern Hemisphere ice sheets and climate during the last glacial cycle, Clim. Past, 19, 399–418, https://doi.org/10.5194/cp-19-399-2023, 2023.

Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res., 112, F03S28, https://doi.org/10.1029/2006JF000664, 2007.

Schoof, C.: Marine ice sheet stability, J. Fluid Mech., 698, 62–72, https://doi.org/10.1017/jfm.2012.43, 2012.

Shepherd, A., Ivins, E., Rignot, E., Smith, B., Van Den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., Nowicki, S., Payne, T., Scambos, T., Schlegel, N., Geruo, A., Agosta, C., Ahlstrøm, A., Babonis, G., Barletta, V., Blazquez, A., Bonin, J., Csatho, B., Cullather, R., Felikson, D., Fettweis, X., Forsberg, R., Gallee, H., Gardner, A., Gilbert, L., Groh, A., Gunter, B., Hanna, E., Harig, C., Helm, V., Horvath, A., Horwath, M., Khan, S., Kjeldsen, K. K., Konrad, H., Langen, P., Lecavalier, B., Loomis, B., Luthcke, S., McMillan, M., Melini, D., Mernild, S., Mohajerani, Y., Moore, P., Mouginot, J., Moyano, G., Muir, A., Nagler, T., Nield, G., Nilsson, J., Noel, B., Otosaka, I., Pattle, M. E., Peltier, W. R., Pie, N., Rietbroek, R., Rott, H., Sandberg-Sørensen, L., Sasgen, I., Save, H., Scheuchl, B., Schrama, E., Schröder, L., Seo, K. W., Simonsen, S., Slater, T., Spada, G., Sutterley, T., Talpe, M., Tarasov, L., Van De Berg, W. J., Van Der Wal, W., Van Wessem, M., Vishwakarma, B. D., Wiese, D., and Wouters, B.: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s41586-018-0179-y, 2018.

Siegenthaler, U., Stocker, T. F., Monnin, E., Luthi, D., Schwander, J., Stauffer, B., Raynaud, D., Barnola, J.-M., Fischer, H., and Masson-Delmotte, V.: Stable carbon cycle climate relationship during the Late Pleistocene, Science, 310, 1313–1317, 2005.

Sigman, D., Fripiat, F., Studer, A. S., Kemeny, P. C., Martínez-García, A., Hain, M. P., Ai, X., Wang, X., Ren, H., and Haug, G. H.: The Southern Ocean during the ice ages: A review
of the Antarctic surface isolation hypothesis, with comparison to the North Pacific, Quaternary Sci. Rev., 254, 106732, https://doi.org/10.1016/j.quascirev.2020.106732, 2021.

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, https://doi.org/10.1016/j.quascirev.2018.12.018, 2019.

Spratt, R. M. and Lisiecki, L. E.: A Late Pleistocene sea level stack, Clim. Past, 12, 1079–1092, https://doi.org/10.5194/cp-12-1079-2016, 2016.

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, https://doi.org/10.5194/cp-10-2135-2014, 2014.

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, https://doi.org/10.5194/gmd-6-819-2013, 2013.

Sugiyama, S., Minowa, M., Sakakibara, D., Skvarca, P., Sawagaki, T., Ohashi, Y., Naito, N., and Chikita, K.: Thermal structure of proglacial lakes in Patagonia, J. Geophys. Res.-Earth, 121, 2270–2286, 2016.

Sun, S., Pattyn, F., Simon, E. G., Albrecht, T., Cornford, S., Calov, R., Dumas, C., Gillet-Chaulet, F., Goelzer, H., Golledge, N. R., Greve, R., Hoffman, M. J., Humbert, A., Kazmierczak, E., Kleiner, T., Leguy, G. R., Lipscomb, W. H., Martin, D., Morlighem, M., Nowicki, S., Pollard, D., Price, S., Quiquet, A., Seroussi, H., Schlemm, T., Sutter, J., van de Wal, R. S. W., Winkelmann, R., and Zhang, T.: Antarctic ice sheet response to sudden and sustained ice-shelf collapse (ABUMIP), J. Glaciol., 66, 891–904, https://doi.org/10.1017/jog.2020.67, 2020.

Tabor, C. R. and Poulsen, C. J.: Simulating the mid-Pleistocene transition through regolith removal, Earth Planet. Sc. Lett., 434, 231–240, https://doi.org/10.1016/j.epsl.2015.11.034, 2016.

Tarasov, L. and Peltier, W. R.: Terminating the 100kyr Ice Age cycle, J. Geophys. Res., 102, 21665–21693, https://doi.org/10.1029/97JD01766, 1997.

Tarasov, L. and Peltier, W.: Arctic freshwater forcing of the Younger Dryas cold reversal, Nature, 435, 662–665, https://doi.org/10.1038/nature03617, 2005.

Tarasov, L. and Peltier, W. R.: A calibrated deglacial drainage chronology for the North American continent: Evidence of an Arctic trigger for the Younger Dryas, Quaternary Sci. Rev., 25, 659–688, 2006.

Teller, J. T., Leverington, D. W., and Mann, J. D.: Freshwater outbursts to the oceans from glacial Lake Agassiz and their role in climate change during the last deglaciation, Quaternary Sci. Rev., 21, 879–887, https://doi.org/10.1016/S0277-3791(01)00145-7, 2002.

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, https://doi.org/10.5194/cp-10-487-2014, 2014.

Upham, W.: The geology of central and western Minnesota. A preliminary report. By Warren Upham, assistant on the Geological and Natural History Survey of the state, under the direction of Prof. N. H. Winchell. [From the General Report of Progress for the Year 1879], St. Paul, The Pioneer Press Co., https://doi.org/10.1017/S0016756800117194, 1881.

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.

van Aalderen, V., Charbit, S., Dumas, C., and Quiquet, A.: Relative importance of the mechanisms triggering the Eurasian ice sheet deglaciation in the GRISLI2.0 ice sheet model, Clim. Past, 20, 187–209, https://doi.org/10.5194/cp-20-187-2024, 2024.

van Calcar, C. J., van de Wal, R. S. W., Blank, B., de Boer, B., and van der Wal, W.: Simulation of a fully coupled 3D glacial isostatic adjustment – ice sheet model for the Antarctic ice sheet over a glacial cycle, Geosci. Model Dev., 16, 5473–5492, https://doi.org/10.5194/gmd-16-5473-2023, 2023.

Velay-Vitow, J., Chandan, D., and Peltier, W. R.: Into the Holocene, anatomy of the Younger Dryas cold reversal and preboreal oscillation, Sci. Rep., 14, 3134, https://doi.org/10.1038/s41598-024-53591-2, 2024.

Warren, C. R. and Kirkbride, M. P.: Calving speed and climatic sensitivity of New Zealand lake-calving glaciers, Ann. Glaciol., 36, 173–178, https://doi.org/10.3189/172756403781816446, 2003.

Warren, C. R., Greene, D. R., and Glasser, N. F.: Glaciar Upsala, Patagonia: Rapid Calving Retreat in Fresh Water, Ann. Glaciol., 21, 311–316, https://doi.org/10.3189/S0260305500015998, 1995.

Weertman, J.: Stability of the Junction of an Ice Sheet and an Ice Shelf, J. Glaciol., 13, 3–11, https://doi.org/10.3189/S0022143000023327, 1974.

Willeit, M., Ganopolski, A., Calov, R., and Brovkin, V.: Mid-Pleistocene transition in glacial cycles explained by declining CO₂ and regolith removal, Sci. Adv., 5, eaav7337, https://doi.org/10.1126/sciadv.aav7337, 2019.