Articles | Volume 20, issue 6
Research article
17 Jun 2024
Research article |  | 17 Jun 2024

New estimation of critical insolation–CO2 relationship for triggering glacial inception

Stefanie Talento, Matteo Willeit, and Andrey Ganopolski

It has been previously proposed that glacial inception represents a bifurcation transition between interglacial and glacial states and is governed by the nonlinear dynamics of the climate–cryosphere system. To trigger glacial inception, the orbital forcing (defined as the maximum of summer insolation at 65° N and determined by Earth’s orbital parameters) must be lower than a critical level, which depends on the atmospheric CO2 concentration. While paleoclimatic data do not provide a strong constraint on the dependence between CO2 and critical insolation, its accurate estimation is of fundamental importance for predicting future glaciations and the effect that anthropogenic CO2 emissions might have on them. In this study, we use the novel Earth system model of intermediate complexity CLIMBER-X with interactive ice sheets to produce a new estimation of the critical insolation–CO2 relationship for triggering glacial inception. We perform a series of experiments in which different combinations of orbital forcing and atmospheric CO2 concentration are maintained constant in time. We analyze for which combinations of orbital forcing and CO2 glacial inception occurs and trace the critical relationship between them, separating conditions under which glacial inception is possible from those where glacial inception is not materialized. We also provide a theoretical foundation for the proposed critical insolation–CO2 relation. We find that the use of the maximum summer insolation at 65° N as a single metric for orbital forcing is adequate for tracing the glacial inception bifurcation. Moreover, we find that the temporal and spatial patterns of ice sheet growth during glacial inception are not always the same but depend on the critical insolation and CO2 level. The experiments evidence the fact that during glacial inception, ice sheets grow mostly in North America, and only under low CO2 conditions are ice sheets also formed over Scandinavia. The latter is associated with a weak Atlantic Meridional Overturning Circulation (AMOC) for low CO2. We find that the strength of AMOC also affects the rate of ice sheet growth during glacial inception.

1 Introduction

Glacial cycles are the dominant mode of climatic variability over the last 2.7 million years. The timing of glaciations and deglaciations is primarily controlled by changes in Earth’s orbital parameters through the modulation of the amount of solar radiation received in high latitudes of the Northern Hemisphere (NH) during boreal summer (Milankovitch1941; Ganopolski2024).

Using an Earth system model of intermediate complexity (EMIC), Calov et al. (2005) proposed that glacial inception represents a bifurcation transition from interglacial to glacial states of the Earth system. When summer insolation at high latitudes of the NH falls below a certain threshold, the interglacial state becomes unstable and ice sheet growth begins. The process is amplified through nonlinear feedbacks, among which the snow and ice albedo and the elevation feedbacks are the dominant ones. Eventually, large ice sheets develop over the North American and Eurasian continents.

The threshold value for boreal summer insolation to trigger a glacial inception depends on the atmospheric concentration of greenhouse gases, among which CO2 is, by far, the most important (Köhler et al.2010). While paleoclimatic data do not constrain this dependence, its accurate estimation is of fundamental importance for predicting future glaciations and the effect that anthropogenic CO2 emissions might have on them (Archer and Ganopolski2005; Talento and Ganopolski2021).

The first attempt at tracing the relationship between insolation and CO2 for triggering glacial inception was performed in Archer and Ganopolski (2005) using the EMIC CLIMBER-2 (Petoukhov et al.2000; Ganopolski et al.2001). This first estimation was, however, rather crude as it was based only on a small set of relatively short experiments that did not allow for the detection of the bifurcation with sufficient accuracy. Ganopolski et al. (2016) updated the results also using CLIMBER-2 with a new methodology that allowed for increased accuracy. The authors proposed a functional relationship between summer maximum insolation at 65° N (smx65) and CO2 for triggering glacial inception. The relationship describes the critical smx65 for the onset of glaciation as linearly dependent on the logarithm of CO2: smx65cr=αln(CO2/CO20)+β, where CO2 is the atmospheric CO2 concentration, CO20 is a reference atmospheric CO2 value (the pre-industrial value), and α and β are empirical constants.

The shape of the dependency is in agreement with the facts that (i) the radiative forcing of CO2 follows a logarithmic structure and (ii) that in CLIMBER-2 the temperature response to radiative forcing of CO2 and orbital forcing is linear within the considered range.

The purpose of this work is to build upon the findings of Ganopolski et al. (2016) and generate a new estimation of the relationship between insolation and atmospheric CO2 critical for triggering glacial inception using the recently developed and more advanced CLIMBER-X model (Willeit et al.2022, 2023, 2024). While CLIMBER-X is also an EMIC that shares with CLIMBER-2 the principles of computational efficiency and usage of a statistical–dynamical atmospheric model, it represents a significant development compared to the CLIMBER-2 model used in previous studies (Ganopolski and Brovkin2017; Willeit et al.2019). In particular, CLIMBER-X has (i) a much higher horizontal resolution in the atmosphere and land models (5° × 5° versus  50° × 10° in a longitude–latitude grid), (ii) a 3D ocean model, (iii) more internally consistent components, and (iv) better treatment of individual processes (e.g., surface energy balance, precipitation, radiative transfer, sea ice dynamics, photosynthesis, vegetation dynamics). Furthermore, we investigate the temporal and spatial patterns of ice sheet growth during glacial inception and the effects of changes in the Atlantic Meridional Overturning Circulation (AMOC) strength on the rate of ice sheets growth. We also analyze the suitability of using the smx65 as a single metric for orbital forcing. A theoretical derivation of the critical insolation–CO2 relationship for triggering glacial inception is provided in Appendix A.

2 Model and experimental setup

2.1 Model description

We employ the new CLIMBER-X EMIC (Willeit et al.2022, 2023, 2024). The model was specifically designed to simulate the Earth system evolution on timescales from decades to hundreds of thousands of years. The climate component of CLIMBER-X consists of the 2.5-dimensional Semi-Empirical dynamical-Statistical Atmosphere Model (SESAM), the 3-dimensional frictional–geostrophic ocean model GOLDSTEIN, the dynamic–thermodynamic SImple Sea Ice Model (SISIM) and the land surface–vegetation model PALADYN. All the models from the climate core are discretized on a longitude–latitude grid with 5° × 5° horizontal resolution and use a daily time step. The ocean model has 23 unequally spaced vertical layers, with a 10 m top layer and layer thickness increasing with depth and reaching 500 m at the ocean bottom. CLIMBER-X incorporates the 3-dimensional thermomechanical ice sheet model SICOPOLIS (Greve et al.2017) via the high-resolution physically based surface energy and mass balance interface SEMIX (Willeit et al.2024) and a basal ice shelf melt module. In this study, SICOPOLIS is applied only to the NH with a 32 km horizontal resolution. The viscoelastic solid Earth model VILMA (Klemann et al.2008; Martinec et al.2018) is used to simulate the bedrock response to changes in loading by solving the sea level equation. A detailed description of the coupling between climate and ice sheet components is provided in Willeit et al. (2024). CLIMBER-X does not simulate internal decadal-scale variability by its design, but it can simulate the forced response at this timescale, such as response to volcanic eruptions and CO2 changes (Willeit et al.2022). The equilibrium climate sensitivity of CLIMBER-X is  3 °C, close to the best-guess estimate from the IPCC (IPCC 20212021).

Even though CLIMBER-X has a relatively coarse spatial resolution and a simplified treatment of some processes, the number of physical processes included in the model and its performance for the present-day climate are comparable with state-of-the-art (CMIP6-type) climate models (Willeit et al.2022, 2023). CLIMBER-X is reasonably successful in reproducing the present-day climatology of variables relevant for the formation of ice sheets. However, near-surface summer temperatures (a key variable for determining ablation of ice sheets) are up to 5 °C too warm over the eastern portion of North America when compared to observations (Willeit et al.2024). Preliminary experiments show that such biases preclude the realistic development of North American ice sheets. Therefore, similarly to Ganopolski et al. (2010), we implemented a 2 m temperature bias correction in the surface mass balance calculation (see Willeit et al.2024). In addition, a uniform offset of 0.5 °C is applied globally in the surface mass balance scheme. As discussed in Willeit et al. (2024), the precipitation biases in CLIMBER-X are less problematic for the simulation of the surface mass balance of NH ice sheets.

2.1.1 Experimental setup

We perform a series of equilibrium experiments in which different orbital parameters (obliquity, eccentricity and precession) and the atmospheric CO2 concentration are maintained constant for the whole duration of the simulation. All NH ice sheets are fully interactive and influence climate through albedo, elevation, sea level (which affects the land–sea mask) and freshwater flux into the ocean (Willeit et al.2024). Since the ocean model in CLIMBER-X is based on the rigid-lid approximation we always enforce the net annual global surface freshwater flux being zero. In the case of interactive ice sheets, we still enforce the global annual surface freshwater flux being zero but then adjust the ocean volume each year to match the change in global ice volume (Willeit et al.2022). The Antarctic ice sheet is not interactive but fixed at its present-day state. The experiments are performed using an acceleration technique with asynchronous coupling between the climate and ice sheet model components. An acceleration factor of 10 is used; i.e., the climate model component is updated every 10 years of the ice sheet model. If we were to calculate global ocean volume change by integrating surface freshwater fluxes, then a 10-fold acceleration would require us to aggregate these fluxes over 10 years, but since in our model the ocean volume change is not driven by surface freshwater fluxes this is not necessary and the freshwater fluxes are computed as in the non-accelerated case. Willeit et al. (2024) demonstrated that acceleration factors up to  10 are acceptable and provide realistic results of the last glacial inception. Each model simulation is run for 100 000 years of the ice sheet model (corresponding to 10 000 years of the climate model) starting from pre-industrial equilibrium conditions with the observed present-day Greenland ice sheet with a uniform ice temperature of 10 °C.

We select 19 combinations of orbital parameters corresponding to the 19 local minima of the smx65 time series of the last 400 kyr, corresponding to values in the range between 430 and 496 W m−2 (Fig. 1). Therefore, selected combinations of orbital parameters have similar precession but very different eccentricity and obliquity values. For each of these 19 sets of orbital parameters, we produce a series of equilibrium experiments designed to detect the critical insolation–CO2 relationship for triggering inception. The procedure is, for each set of orbital parameters, as follows.

  1. We start from a simulation using a (fairly high) atmospheric CO2 concentration of 500 ppm, for which glacial inception is not expected to occur under any orbital configuration.

  2. If the NH ice volume (v) is <X m s.l.e. (meters in sea level equivalent) during the whole simulation, we reduce CO2 by 5 ppm and repeat the simulation.

If, on the contrary, v is Xms.l.e. at any moment in the simulation, we stop the iteration and define the insolation–CO2 combination used in this case as critical for glacial inception.

In this study, we select a value of X=15 m s.l.e. That value corresponds to an ice volume roughly twice the value of the present-day Greenland ice sheet volume. As in CLIMBER-X ice volume over Greenland does not increase substantially under glacial inception conditions, and 15 ms.l.e. is thus associated with considerable ice sheet growth outside Greenland.

To validate the model suitability for tracing glacial inception bifurcation transitions, similarly to Ganopolski et al. (2016), we test the model performance against two paleoclimatic constraints. The first constraint is that, for a present-day orbital configuration and atmospheric CO2 concentration of 280 ppm (named the Control simulation, Table 1), the model should not produce any substantial ice growth outside of Greenland. The second constraint is that, on the contrary, for the same 280 ppm level of atmospheric CO2 concentration but an orbital configuration corresponding to the end of MIS11 (398 kyr BP), significant ice growth, indicative of glacial inception, should occur. Under the current model version these two constraints are met: after 100 kyr of simulation the NH ice volume totals, including the Greenland ice sheet, 8.3 ms.l.e. for the present-day case and 21 ms.l.e. for the MIS11 configuration.

Additional experiments with prescribed present-day ice sheets or applying freshwater hosing and extraction in the North Atlantic Ocean to alter AMOC strength are also conducted for selected orbital forcing–CO2 configurations (Table 1).

Figure 1Summer maximum insolation at 65° N (smx65) in the period 400 kyr BP–present, computed using Laskar et al. (2004) data for Earth orbital parameters. Marked in red are all 19 local minima of the time series, whose orbital configurations are selected for the CLIMBER-X experiments.


Table 1Summary of selected experiment configurations.

Download Print Version | Download XLSX

3 Results

3.1 Critical smx65–CO2 relationship for triggering glacial inception

The usage of smx65 as a proxy for orbital forcing is common practice and considered a good metric for the analysis of glacial inception processes (Ganopolski et al.2016; Paillard1998; Leloup and Paillard2022), as 65° N latitude coincides with the presence of land potentially capable of supporting the initial growth of ice sheets.

Considering smx65 to be an orbital forcing indicator, Fig. 2 shows the critical smx65–CO2 combinations for triggering glacial inception in the 19 individual sets of experiments. Accounting for the fact that the radiative effect of CO2 is logarithmic and that the temperature response to insolation is linear, we fit the data (via least squares) with curves of the shape

(1) smx 65 cr = α ln CO 2 280 + β .

Figure 2CLIMBER-X estimation of critical summer maximum insolation at 65° N (smx65) and CO2 for triggering glacial inception at the threshold of 15 ms.l.e. (red dots). The red line corresponds to the least-square fit following a logarithmic shape. The CLIMBER-2 estimation from Ganopolski et al. (2016) is depicted in blue. Big red circles indicate the experiments Hsmx65_LCO2, Ismx65_ICO2 and Lsmx65_HCO2.


The least-squares fit of the data points produces values of α=75 W m−2 and β=465W m−2. The fit is good, with a coefficient of determination R2=0.95. The largest error in the fit is produced for an intermediate insolation case (462 W m−2, corresponding to the minimum in summer solstice insolation at 209 ka), in which the difference between the least-squares fit and the CLIMBER-X result for the critical insolation is 8 W m−2. This deviation from the best fit is explained by a very high obliquity at that time (see Sect. 3.3). The new CLIMBER-X estimation for α and β is very similar to the one generated with the CLIMBER-2 model and presented in Ganopolski et al. (2016). For that study the criterion for defining glacial inception was comparable (growth of the NH ice sheets of more than 13.5 ms.l.e. was required), and the values for α and β were 77 and 466 W m−2, respectively. The robustness of the estimated parameters of the critical insolation–CO2 relationship is also supported by the theoretical analysis presented in Appendix A, where it is shown that for climate conditions similar to pre-industrial, the value of α is about 80 W m−2. Moreover, the value derived from the results presented by Abe-Ouchi et al. (2013) is 83 W m−2 (note that this number was not reported in Abe-Ouchi et al.2013, but can be calculated from their Fig. 2). As was shown in Ganopolski et al. (2016), the value of β is reasonably well constrained by paleoclimate data since the critical insolation curve must pass between rather close insolation values corresponding to the end of MIS11 and present insolation. At the same time, paleo-data provide no constraint on the value of α. This is why, in Talento and Ganopolski (2021), we used a very conservative approach by accepting as “valid” any α values from the range 150 to 0 W m−2; i.e., we assumed relative uncertainties of up to 100 %. The results of the present study strongly indicate that the uncertainty is much smaller, likely not higher than 20 %. Such a reduction of the uncertainty range would also significantly reduce the uncertainties in the projections of the timing of future glaciations for different anthropogenic CO2 emissions.

In order to better understand glacial inception at the critical insolation–CO2 threshold and to evaluate similarities and differences for diverse insolation–CO2 pairings, we select three experiments and analyze the temporal and spatial evolution of ice sheet growth. The three experiments (Hsmx65_LCO2, Lsmx65_HCO2 and Ismx65_ICO2) correspond to either the highest, lowest or an intermediate value of the 19 smx65 time series local minima, respectively. In each experiment, the atmospheric CO2 corresponds to the critical CO2 for triggering glacial inception (Table 1). While, by definition, all three experiments Hsmx65_LCO2, Lsmx65_HCO2 and Ismx65_ICO2 reach at least 15 ms.l.e. NH ice volume at some point during the simulation, they do it with different temporal and spatial dynamics (Figs. 3 and 4).

Figure 3NH ice volume (m s.l.e.) temporal evolution at the critical combination of smx65 and CO2 for triggering glacial inception (threshold 15 ms.l.e. level, dashed line) for three selected experiments (see Table 1).


For the experiment Hsmx65_LCO2 (smx65 =496W m−2 and CO2= 190 ppm) the NH ice volume time series crosses the 15 ms.l.e. threshold 67 kyr into the simulation and reaches a quasi-equilibrium value of 16.6 ms.l.e. only after 90 kyr (Fig. 3). At the moment of crossing the glacial inception threshold (15 ms.l.e.), the ice sheets have significantly developed outside of Greenland over Iceland, Svalbard, Novaya Zemlya and the Scandinavian Peninsula as well as in North America, north of Hudson Bay, from 90 to 70° W (Fig. 4a). The Scandinavian ice sheet development is aided by a weak AMOC (Fig. 5) and subsequent weakened meridional heat transport from the tropics towards that region. The AMOC reaches quasi-equilibrium conditions in the second half of the simulation, with a maximum magnitude of  16 Sv. At the threshold crossing time, the AMOC is weak and shallow (Fig. 6a) and deep convection in the North Atlantic Ocean only occurs in the Norwegian Sea (Fig. 6d).

Figure 4Ice thickness (a–c) for Hsmx65_LCO2, Ismx65_ICO2 and Lsmx65_HCO2 at the time when each experiment reaches the 15 ms.l.e. threshold.

Figure 5Maximum Atlantic overturning (Sv) temporal evolution at the critical combination of smx65 and CO2 for triggering glacial inception (threshold 15 ms.l.e. level) for three selected experiments (see Table 1).


In the experiment Ismx65_ICO2 (smx65 =464.8W m−2 and CO2= 285 ppm) the NH ice volume time series shows a faster increasing rate at the start of the experiment than the Hsmx65_LCO2 experiment but stagnates at 14.5 ms.l.e. between 40 and 60 kyr of simulation. At 60 kyr of simulation, the time series shows an abrupt increase, crossing the 15 ms.l.e. threshold at 63 kyr. After that, the time series stagnates again at 17.5 ms.l.e. before showing another abrupt increase towards the end of the simulated period. In the 100 kyr of simulation, quasi-stationary conditions are not reached (Fig. 3). In this experiment, the AMOC is unstable and oscillates between  25 and  16 Sv (Fig. 5). It should be noted that the timescale of these simulated oscillations is distorted by the climate acceleration factor. The AMOC weakening episodes precede the abrupt NH ice volume growth and are thus likely triggers. At the moment when the NH ice volume is 15 ms.l.e., ice sheets cover the north of Hudson Bay from 100 to 70° W (covering almost the totality of Baffin Island), and there is full ice coverage of Iceland, Svalbard and Novaya Zemlya and partial ice coverage of Scandinavia (Fig. 4c). The moderate ice growth over Scandinavia is related to the AMOC weakening episodes (Fig. 5) and the consequential decreased heat transport from the tropics, as such ice growth is inhibited in simulations with a stable AMOC (see Sect. 3.2 with extraction simulation Ismx65_ICO2 0.01 Sv). At the time of the 15 ms.l.e. crossing, the Atlantic overturning circulation penetrates to  3000 m of depth and deep-water formation does not occur in the Labrador Sea and south of Greenland (Fig. 6b, e).

The Lsmx65_HCO2 experiment (smx65 =430.8W m−2 and CO2= 460 ppm) shows the fastest NH ice volume time series increase at the start of the simulation, crossing the 15 ms.l.e. threshold after only 18 kyr of simulation. After that, the ice volume slightly decreases, stabilizing at 14.5 ms.l.e. (Fig. 3). At the moment of crossing of the threshold, the ice sheet coverage over North America is larger than in the other two experiments, while no ice sheet develops over Scandinavia (Fig. 4c). For this case, the AMOC reaches equilibrium with a maximum strength of 27 Sv in the second half of the simulation (Fig. 5). Deep convection is widespread in the North Atlantic, occurring over the Labrador Sea, south of Greenland and the Norwegian Sea (Fig. 6f). The Atlantic overturning circulation reaches depths of 4000 m (Fig. 6c).

From these experiments it is clear that the timescales involved with glacial development in the vicinity of the glacial inception bifurcation point are extremely long and that even 100 kyr of simulation might not be enough to reach stationary conditions. It is also clear that, in CLIMBER-X, glacial inception can occur for different AMOC configurations. The spatial patterns of glaciation are different and related to the AMOC intensity, with a weak AMOC transporting less heat from the tropics towards Scandinavia and facilitating ice growth there. We see that in quasi-equilibrium conditions higher atmospheric CO2 concentrations produce a stronger AMOC in CLIMBER-X. This was also observed in CLIMBER-2 (Bouttes et al.2012; Gottschalk et al.2019) and general circulation models, starting from Stouffer and Manabe (2003), who found AMOC strengthening under doubling and quadrupling of CO2 and significant AMOC weakening under CO2 halving in multimillennial simulations. Recently, Bonan et al. (2022) described the AMOC response to instantaneous CO2 quadrupling in several state-of-the-art climate models. Mots of the runs were not long enough, and the AMOC continues to evolve through the runs, but at least in one model (CESM1), the AMOC is appreciably stronger under CO2 quadrupling. For a lower CO2 concentration, there are two systematic studies (Oka et al.2012; Galbraith and de Lavergne2019) where it has been shown that global cooling causes a weakening of the AMOC.

Figure 6Annual mean Atlantic Ocean meridional overturning circulation streamfunction (a–c) and maximum (March) mixed layer depth (d–f) for the experiments Hsmx65_LCO2, Ismx65_ICO2 and Lsmx65_HCO2 at the time when each experiment reaches the 15 ms.l.e. threshold.

To get a better understanding on how different spatial patterns of ice growth arise with different combination of insolation and CO2, we analyze the climatic conditions leading to the ice development. Figure 7 displays the annual mean surface mass balance (SMB) anomalies with respect to Control at the moment when the simulations reach the 15 ms.l.e. of NH ice volume. For the three experiments, positive SMB anomalies with respect to Control are found over Greenland and over Baffin Island, the anomalies being larger in magnitude for the Lsmx65_HCO2 case. For the Hsmx65_LCO2 experiment, the largest positive anomalies occur over the Scandinavian Peninsula.

Figure 7Annual surface mass balance anomaly with respect to Control for the experiments Hsmx65_LCO2, Ismx65_ICO2 and Lsmx65_HCO2 (a–c), respectively, at the time when each experiment reaches the 15 ms.l.e. threshold. Continental and ice sheet margins are show in black.

Given the fact that surface temperature and precipitation fields might be directly impacted by the presence of an ice sheet, we generate auxiliary experiments in which the conditions are the same as in Hsmx65_LCO2, Ismx65_ICO2 and Lsmx65_HCO2 but without interactive ice sheets (i.e., the ice sheets remain constant and equal to the initial present-day conditions; Table 1) and denote them with an additional FixIce label.

Boreal summer (June–August; JJA) 2 m temperature and annual precipitation anomalies with respect to Control in simulations without interactive ice sheets for the NH between 120° W and 60° E are displayed in Fig. 8. The largest temperature differences between experiments are seen over the northern North Atlantic and are related to the differences in AMOC strength and related northward oceanic heat transport. The AMOC is the weakest in Hsmx65_LCO2_FixIce, and as a result, Scandinavia is colder in this experiment by 2–4 °C than in the Control experiment. This explains why the ice sheets developed in Scandinavia in Hsmx65_LCO2 but not in Lsmx65_HCO2. Over northern North America, temperatures are rather similar in all experiments, which is explained by the fact that summer temperatures are the main factor controlling the mass balance of ice sheets. Still, temperatures in Hsmx65_LCO2 are lower by ca. 2 °C than in Lsmx65_HCO2_FixIce. The reason why glacial inception requires lower temperatures in experiments with higher insolation is straightforward. Unlike CO2, insolation affects the surface mass balance of ice sheets not only through temperature but also directly through absorbed shortwave radiation (e.g., Robinson et al.2010). This is why, to compensate for the effect of higher insolation on surface ice sheet melt in the experiment Hsmx65_LCO2, summer temperatures should be lower than in the experiments where glacial inception is simulated under lower summer insolation.

Figure 8Summer (JJA) temperature anomalies (a–c) and annual precipitation anomalies (d–f) with respect to Control_FixIce for experiments Hsmx65_LCO2_FixIce, Ismx65_ICO2_FixIce and Lsmx65_CO2_FixIce. Control continental margins (present-day conditions) are shown in black.

Annual precipitation anomalies relative to Control shown in Fig. 8d–f also mostly reflect changes in AMOC, although globally there is more precipitation under a higher CO2 level. Over Scandinavia, precipitation changes counteract temperature changes, but the effect of temperature is dominant. This is why, in spite of significantly lower annual precipitation, the ice sheet in Scandinavia grows in the experiment Hsmx65_LCO2 but not in Lsmx65_HCO2. Precipitation differences in northern North America are small.

In general, a comparison of these three experiments shows that summer climate conditions corresponding to the initiation of glaciation are similar in North America despite very different insolation and CO2 in these experiments. This explains rather similar patterns of initial ice sheet growth in North America. Climate conditions are more different in Scandinavia, and this is why glaciation develops here simultaneously with glaciation in North America only under sufficiently low CO2.

3.2 Role of AMOC change in glacial inception

As shown above, the AMOC state is different for different combinations of insolation–CO2 close to the glacial inception limit (Fig. 6a–c). The question therefore arises of what role the AMOC plays for glacial inception. To explore the relative importance of the AMOC state compared to the summer insolation and the atmospheric CO2 concentration in shaping the conditions for glacial inception, we have performed additional simulations in order to separate the contribution of the different factors. For that we have taken the extreme cases, Hsmx65_LCO2 and Lsmx65_HCO2, and run three sets of model simulations with prescribed present-day ice sheets:

  • present-day orbital configuration with CO2 concentration corresponding to HCO2 and LCO2,

  • pre-industrial CO2 concentration with orbital configuration corresponding to Hsmx65 and Lsmx65,

  • pre-industrial CO2 concentration and orbital configuration with and without 0.1 Sv freshwater hosing in the North Atlantic between 50–70° N to get two states with weak and strong AMOC ( 10 Sv difference roughly corresponding to the difference between Hsmx65_LCO2 and Lsmx65_HCO2 experiments) under the same background climate conditions.

These three sets of two simulations each allow separating the effect of atmospheric CO2, summer insolation and AMOC state on the summer temperature in the NH. The results clearly show that summer temperatures at the locations where the ice sheets start to grow during glacial inceptions (i.e., Scandinavia and northeastern Canada) are much more sensitive to CO2 and orbital forcing than to changes in the AMOC (Fig. 9). The AMOC state has a generally larger impact on summer temperatures over Scandinavia (Fig. 9c), at least partly explaining why the inception patterns show the largest differences there (Fig. 4). It has to be noted that changes in CO2 and insolation also affect the AMOC. In particular, the AMOC is stronger by 9 Sv in the experiment with high CO2 compared to low CO2 and by 3 Sv in experiments with high insolation compared to low insolation. However, as seen from Fig. 9c, these changes can contribute only a little to the direct effect of CO2 (Fig. 9a) and insolation (Fig. 9b). Thus, at the potential locations of glacial inception, AMOC changes provide a positive, but not very strong, feedback to both primary forcings.

Figure 9Summer (JJA) temperature differences due to (a) CO2 (HCO2–LCO2), (b) insolation (Hsmx65–Lsmx65) and (c) AMOC state (strong–weak).

To further investigate the role of AMOC strength change at the glacial inception bifurcation limits, we produce a series of water hosing and extraction experiments (Table 1). The freshwater forcing is applied in the North Atlantic Ocean between 50 and 70° N and compensated for over the tropical (30° S–30° N) Pacific Ocean to prevent a global salinity drift. The magnitude of the freshwater flux is 0.01 Sv, either positive (added freshwater, leading to a weakening of AMOC) or negative (extracted freshwater, leading to a strengthening of AMOC). The hosing and extraction forcing is applied for the whole duration of the experiments. Note that 0.01 Sv is a very small value compared to a typical stability threshold for the AMOC in the CLIMBER-X model (Willeit et al.2022).

The freshwater extraction produces less NH ice sheet volume in all three cases (Fig. 10). In particular, in the experiments Hsmx65_LCO2 and Ismx65_ICO2 the extraction in enough to prevent the NH ice volume from crossing the 15 ms.l.e. threshold used for the definition of glacial inception. In opposition, the freshwater hosing flux produces in all three cases faster and larger ice growth in the NH. The largest change in terms of NH ice volume is produced for Ismx65_ICO2, reaching a magnitude of  150 ms.l.e. In this case, the freshwater hosing flux produces the largest AMOC weakening of all three experiments and even an almost complete AMOC shutdown between 50 and 80 kyr of the simulation (Fig. 11).

Figure 10NH ice volume (m s.l.e.) temporal evolution in water hosing and extraction experiments (Table 1).


Figure 11Maximum Atlantic overturning (Sv) temporal evolution in water hosing and extraction experiments (Table 1).


This set of hosing experiments shows that the AMOC (in particular its stability) plays a role in triggering glacial inception and the rate of ice growth. In general, the weaker the AMOC the larger and faster the development of the NH ice sheets for a given insolation–CO2 pairing.

3.3 Suitability of smx65 as a proxy for orbital forcing

The usage of smx65 as a proxy for orbital forcing in tracing the glacial inception bifurcation is justified by the relevance of 65° N latitude for the development of ice sheets (Ganopolski2024). However, as with any simple metric for such a complex phenomenon as the influence of orbital forcing on Earth climate, it cannot be perfect. This is why it is not surprising that although most of the points in Fig. 2 are located very close to the proposed logarithmic curves, several deviate more significantly, for example the point corresponding to insolation at 209 ka, as it was discussed above.

To evaluate the possibility of improving the proxy for orbital forcing, we separate the smx65 time series into its obliquity and precessional components (as in Jackson and Broccoli2003). The decomposition is done via a least-squares fit to the data of the last 400 kyr, expressing the smx65 time series in the form

(2) smx 65 = a 0 + a 1 ecc sin ω + a 2 ( obl - obl ) ,

with ecc, ω and obl denoting eccentricity, longitude of perihelion and obliquity, respectively; coefficients a0, a1 and a2 are derived from the least-squares fit and the overbar indicates time averaging. The precessional (smx65p) and obliquity (smx65o) components are defined as


We analyze different linear combinations of smx65p and smx65o as alternative proxies for orbital forcing. We wish to find γ, which maximizes the goodness of fit (R2 coefficient of determination) between the logarithm of critical CO2 and the generalized proxy for summer insolation in the form smx65p+γsmx65o.

Any value of γ between 1 and 2 yields an extremely good fit (R2 higher than 0.95), and thus the quantities smx65p+γsmx65o for those values of γ are good proxies for orbital forcing (Fig. 12). Values of γ lower than 1, in which the contribution of obliquity is considered relatively less important than that of precession, are not optimal proxies. Similarly, values of γ higher than 2 (i.e., the obliquity contribution is weighted more than the precessional one) are also not adequate. The best fit is found for γ≈1.3. Although this fit is only marginally better than for γ=1, corresponding to smx65, and thus of little importance for the main results presented in this paper, it still can help to explain some deviations of individual results from the logarithmic curve. Indeed, as was discussed in Ganopolski (2024), the smx65 metric contains the smallest contribution of obliquity compared to other proposed proxies for orbital forcing. The results of the analysis presented above suggest that the relative contribution of obliquity is underestimated in smx65 by ca. 30 %. Since the amplitude of obliquity contribution to smx65 is only about 15 W m−2, the generalized proxy for summer insolation for γ=1.3 would not deviate from smx65 by more than 5 W m−2, which is very small compared to smx65 values. However, it has an implication for the “outlier” in Fig. 2 corresponding to insolation at 209 ka, when obliquity was close to its maximum value during the entire past 800 kyr. The correction upwards by 5 W m−2 would bring the point corresponding to 209 ka insolation much closer to the logarithmic curve.

Figure 12Goodness of fit (R2 coefficient of determination) derived from the least-squares fit between ln (CO2,critical) and smx65p+γsmx65o as a function of γ, considering CO2 values critical for glacial inception at the 15 ms.l.e. level.


4 Summary and conclusions

We used the recently developed EMIC CLIMBER-X to produce a new estimation of the critical insolation–CO2 relationship for triggering glacial inception. The use of CLIMBER-X for tracing the glacial inception bifurcation constitutes a further development from previous attempts produced with the coarser-resolution and simpler CLIMBER-2 model (Ganopolski et al.2016; Archer and Ganopolski2005).

The study is based on sets of equilibrium experiments in which the orbital parameters and atmospheric CO2 concentration are kept constant in time. The sets are based on orbital configurations corresponding to all the local minima of summer maximum insolation at 65° N (smx65) of the last 400 kyr. Therefore, our experiments reflect a diverse arrangement of orbital parameters in configurations potentially favorable for glacial inception. The estimation of the glacial inception bifurcation is done within a 5 ppm resolution in the CO2 phase space.

The new estimation of the critical insolation–CO2 relationship for triggering glacial inception follows a logarithmic dependence, explained by the logarithmic shape of the radiative forcing of CO2. When considering that glacial inception occurs if the NH ice volume develops at least 15 ms.l.e., the new CLIMBER-X estimation for the critical summer maximum insolation at 65° N and CO2 relationship is

(5) smx 65 cr = - 75 ln CO 2 280 + 465 .

This new estimation is close to the one produced with CLIMBER-2 in Ganopolski et al. (2016), who used a similar threshold for NH ice volume. Moreover, we showed that the summer maximum insolation at 65° N is a skillful single metric for tracing the glacial inception bifurcation.

We also showed that the temporal and spatial patterns of glacial inception depend on the combination of insolation and critical CO2 concentration. In all simulations of glacial inception, ice sheets developed north of Hudson Bay in North America, as well as in Iceland and Svalbard. Considerable ice growth over the Scandinavian Peninsula is only observed in experiments in which the Atlantic Meridional Overturning Circulation (AMOC) significantly weakens, which occurs in either low-CO2 cases (stable weak AMOC) or intermediate-CO2 situations (unstable AMOC with substantial weakening episodes). In general, the lower the CO2, the weaker the AMOC and, thus, the weaker the heat transport from the tropics towards European high latitudes, facilitating the Scandinavian ice sheet development. The time needed for reaching the glacial inception bifurcation threshold varies between 18 000 and almost 100 000 years, highlighting the long timescales of the climate–cryosphere system at the bifurcation limit.

Appendix A: Theoretical derivation of the critical insolation–CO2 relationship

To derive the critical insolation–CO2 relationship for triggering glacial inception, we consider the surface energy balance of snow cover at the location where the nucleation of the northern continental ice sheet begins during glacial inceptions. Important assumptions are that

  • this location remains the same for the different combinations of insolation and CO2 and that

  • glacial inception begins from the formation of extensive perennial snowfields rather than from the spreading of glaciers from high mountains.

Both these assumptions are consistent with the results of our model simulations. Surface daily mass balance of the snowpack during the melt season is described by the following equation:

(A1) S abs + R - R + H = M L ,

where Sabs is the absorbed shortwave radiation, R is the downward longwave radiation flux near the Earth's surface, R is the longwave radiation emitted by the snow surface, H is the sensible heat flux (positive downward), M is the snowmelt rate in kgm-2s-1 and L is the latent heat of fusion. Surface latent heat flux can be neglected during the melt season, as it is an order of magnitude smaller than shortwave and longwave radiation fluxes (e.g., Ettema et al.2010). Since we consider only the period of snowmelt, surface temperature during this period is 0 °C. In this case, the terms on the left-hand side of Eq. (A1) can be expressed as follows:


where T is the surface air temperature (in °C), S is the insolation at the top of the atmosphere, αs is the surface albedo, τ is the atmospheric transmissivity, F(0) is the downward longwave radiation for a surface air temperature of 0 °C in W m−2, T0=273.15K, σ=5.67×10-8Wm-2K-4 is the Stefan–Boltzmann constant, ϵ 1 is the snow emissivity for longwave radiation, and α and β are empirical parameters. In the nucleation center of glacial inception (northern Canada and Baffin Island), the melt season is only about 2 months and approximately coincides with the period of positive air temperatures. This allows us to make further simplifications when integrating Eq. (A1) over the melt season after substituting Eqs. (A2), (A3), (A4) and (A5):

(A6) S ( 1 - α s ) τ + F ( 0 ) + β T - σ T 0 4 + γ T = M L ,

where T is the surface air temperature, S is daily insolation and M is snowmelt, all averaged over the melt season. Assuming that the duration of the melt season at the location of glacial inception is about 2 months (which is similar to present climate conditions in the area where glacial inception is simulated in CLIMBER-X), the seasonal surface air temperature evolution can be approximated by a sinusoidal function with zero air temperatures at the beginning and end of the melt season and the average insolation and temperature can be expressed through their maximum annual values as


where kS 0.9 and kT 0.7 (see Fig. A1). Maximum summer air temperature is parameterized as

(A9) T max = T max * + δ ln CO 2 CO 2 0 + μ ( smx 65 - smx 65 * ) ,

where Tmax* is the maximum summer temperature at the location of glacial inception for pre-industrial conditions, namely for CO2= CO20 and smx65=smx65*, where CO2 is the atmospheric CO2 concentration, CO20= 280 ppm and smx65*=479W m−2 is the present-day maximum insolation at 65° N. The parameter μ represents the sensitivity of summer temperature at the location of glacial inception to summer insolation, i.e.,

(A10) μ = Δ T max Δ smx 65 .

Figure A1Typical seasonal variations of surface air temperature (a) and insolation (b) at the location of glacial inception.


The parameter δ in Eq. (A9) represents the sensitivity of maximum summer temperature to CO2 concentration. It can be expressed through other known parameters as δ=parfcs, where pa is the “polar amplification” of maximum summer temperature equal to ΔTmax/ΔTg, where ΔTg is the global air temperature anomaly caused by CO2 concentration change; rf is the parameter determining the radiative forcing of CO2,

(A11) Δ R CO 2 = r f ln CO 2 CO 2 0 ;

and cs is the specific climate sensitivity (e.g., Rohling et al.2012), which relates the global temperature change to radiative forcing:

(A12) Δ T g = c s Δ R .

Based on IPCC 2021 (2021) and climate model simulations pa=1.5, rf=5.7W m−2 and cs=0.75°C W−1 m2, which is equivalent to an equilibrium climate sensitivity of 3 °C. Thus δ=6.4°C. To derive a critical relationship between insolation and CO2 for triggering glacial inception, we make use of the fact that the glacial inception occurs when the total snowmelt minus refreezing during the melt season does not exceed annual snowfall P (in kg m−2), which for the considered region is less than 200 kg m−2:

(A13) M Δ t M = P 1 - f r ,

where ΔtM is the duration of the melt season and fr is an average refreezing fraction of snowmelt. Thus, Eq. (A6) can be rewritten as

(A14) k S ( 1 - α s ) τ smx 65 + k T ( β + γ ) δ ln CO 2 CO 2 0 + μ smx 65 = P L ( 1 - f r ) Δ t M + σ T 0 4 - F ( 0 ) - k T ( β + γ ) T max * - μ smx 65 * ,

which can be rewritten in the form

(A15) smx 65 cr = A - B ln CO 2 CO 2 0 .

This is identical to the relationship proposed in Ganopolski et al. (2016). Here

(A16) A = P L ( 1 - f r ) - 1 Δ t M - 1 + σ T 0 4 - F ( 0 ) + k T ( β + γ ) ( μ smx 65 * - T max * ) k S ( 1 - α s ) τ + k T μ ( β + γ )


(A17) B = k T δ ( β + γ ) ( k S ( 1 - α s ) τ + k T μ ( β + γ ) .

We will start from determining B, which can be calculated using the following values of model parameters: αs=0.7, which is a typical albedo of melting snow; τ=0.5 is taken from Robinson et al. (2010), β=4Wm-2°C-1 (derived from the regression of model data for the considered region), γ=15 W m−2 °C−1 (middle of the range from Braithwaite2009) and μ=0.08 W−1 m2 °C (from climate model simulations for different orbital parameters). However, since kS(1-αs)τ=0.13kT(β+γ)=0.9, Eq. (A17) for B can be significantly simplified, and only two parameters are required:

(A18) B = δ μ 80 W m - 2 .

Equation (A16) for A contains a number of parameters which are not accurately known. However, knowledge of all these parameters is not required. It is only essential that the value of A for different glacial inceptions should remain approximately constant. This requirement is satisfied since only the total snowfall depends on the values of insolation and CO2 in this formula. However, the contribution of the latent heat term PL(1−frtM to the value of A is only about 10 %–20 %, and, according to climate model simulations, changes in annual precipitation for a wide range of CO2 concentrations are of the same magnitude. Thus, the total variations of A across the considered range of climate conditions should not exceed a few percent. This allows us to determine A using two paleoclimate constraints, namely that glacial inceptions occurred at the end of MIS19 and 11, while no glacial inception occurred in the recent past (see also Fig. 3b in Ganopolski et al.2016). These empirical constraints suggest that under the present-day maximum summer insolation of 479 W m−2, glacial inception can only occur if the CO2 concentration is around 240 ppm. This gives the following value:

(A19) A = 479 + 80 ln 240 280 467 W m - 2 .

Note that the choice of critical CO2 for present-day insolation of 240 ± 20 ppm results in an uncertainty of A of only ±6 W m−2. Equation (A15) now can be rewritten as

(A20) smx 65 cr [ W m - 2 ] = - 80 ln CO 2 [ ppm ] 280 + 467 ,

which is numerically very similar to Ganopolski et al. (2016) and the results of the present study.

Equation (A18) for B has a clear physical sense: the slope of the critical relationship between insolation and the logarithm of CO2 concentration is equal to the ratio between local temperature sensitivities to CO2 (δ) and summer insolation (μ). This value has dimension W m−2. If one wants to compare the effect of insolation to the effect of radiative forcing of CO2 in the same units of W m−2, then it would be B/rf=15; i.e., 1 W m−2 of CO2 forcing is approximately equivalent to 15 W m−2 of the maximum summer insolation in determining the glacial inception. This number is consistent with the results of earlier studies (Calov and Ganopolski2005; Abe-Ouchi et al.2013).

Code and data availability

Code and data used in this study can be obtained from (Talento, 2023).

Author contributions

ST performed the model simulations and made the figures, with contributions from MW. All authors contributed to the analysis of results and preparation of the manuscript.

Competing interests

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


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims 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.


The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research, and the Land Brandenburg for supporting this project by providing resources on the high-performance computer system at the Potsdam Institute for Climate Impact Research.

Financial support

This research has been supported by the Nationale Genossenschaft für die Lagerung radioaktiver Abfälle (NAGRA), the Swiss National Cooperative for the Disposal of Radioactive Waste, and the German climate modeling project PalMod supported by the German Federal Ministry of Education and Research (BMBF) as a Research for Sustainability initiative (FONA) (grant nos. 01LP1920B, 01LP1917D, 01LP2305B).

The publication of this article was funded by the Open Access Fund of the Leibniz Association.

Review statement

This paper was edited by Qiuzhen Yin and reviewed by two anonymous referees.


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

Archer, D. and Ganopolski, A.: A movable trigger: Fossil fuel CO2 and the onset of the next glaciation, Geochem. Geophys. Geosy., 6, 1–7,, 2005. a, b, c

Bonan, D. B., Thompson, A. F., Newsom, E. R., Sun, S., and Rugenstein, M.: Transient and Equilibrium Responses of the Atlantic Overturning Circulation to Warming in Coupled Climate Models: The Role of Temperature and Salinity, J. Climate, 35, 5173–5193,, 2022. a

Bouttes, N., Roche, D. M., and Paillard, D.: Systematic study of the impact of fresh water fluxes on the glacial carbon cycle, Clim. Past, 8, 589–607,, 2012. a

Braithwaite, R. J.: Calculation of sensible-heat flux over a melting ice surface using simple climate data and daily measurements of ablation, Ann. Glaciol., 50, 9–15,, 2009. a

Calov, R. and Ganopolski, A.: Multistability and hysteresis in the climate-cryosphere system under orbital forcing, Geophys. Res. Lett., 32, L21717,, 2005. a

Calov, R., Ganopolski, A., Claussen, M., Petoukhov, V., and Greve, R.: Transient simulation of the last glacial inception. Part I: glacial inception as a bifurcation in the climate system, Clim. Dynam., 24, 545–561,, 2005. a

Ettema, J., van den Broeke, M. R., van Meijgaard, E., and van de Berg, W. J.: Climate of the Greenland ice sheet using a high-resolution climate model – Part 2: Near-surface climate and energy balance, The Cryosphere, 4, 529–544,, 2010. a

Galbraith, E. and de Lavergne, C.: Response of a comprehensive climate model to a broad range of external forcings: relevance for deep ocean ventilation and the development of late Cenozoic ice ages, Clim. Dynam., 52, 653–679,, 2019. a

Ganopolski, A.: Toward generalized Milankovitch theory (GMT), Clim. Past, 20, 151–185,, 2024. a, b, c

Ganopolski, A. and Brovkin, V.: Simulation of climate, ice sheets and CO2 evolution during the last four glacial cycles with an Earth system model of intermediate complexity, Clim. Past, 13, 1695–1716,, 2017. a

Ganopolski, A., Petoukhov, V., Rahmstorf, S., Brovkin, V., Claussen, M., Eliseev, A., and Kubatzki, C.: CLIMBER-2: a climate system model of intermediate complexity. Part II: model sensitivity, Clim. Dynam., 17, 735–751,, 2001. a

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

Ganopolski, A., Winkelmann, R., and Schellnhuber, H. J.: Critical insolation–CO2 relation for diagnosing past and future glacial inception, Nature, 529, 200–203,, 2016. a, b, c, d, e, f, g, h, i, j, k, l

Gottschalk, J., Battaglia, G., Fischer, H., Frölicher, T. L., Jaccard, S. L., Jeltsch-Thömmes, A., Joos, F., Köhler, P., Meissner, K. J., Menviel, L., Nehrbass-Ahles, C., Schmitt, J., Schmittner, A., Skinner, L. C., and Stocker, T. F.: Mechanisms of millennial-scale atmospheric CO2 change in numerical model simulations, Quaternary Sci. Rev., 220, 30–74,, 2019. a

Greve, R., Calov, R., and Herzfeld, U. C.: Projecting the response of the Greenland ice sheet to future climate change with the ice sheet model SICOPOLIS, Low Temperature Science, 75, 117–129,, 2017. a

IPCC 2021: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, in press, 2021. a, b

Jackson, C. S. and Broccoli, A. J.: Orbital forcing of Arctic climate: Mechanisms of climate response and implications for continental glaciation, Clim. Dynam., 21, 539–557,, 2003. a

Klemann, V., Martinec, Z., and Ivins, E. R.: Glacial isostasy and plate motion, J. Geodyn., 46, 95–103,, 2008. a

Köhler, P., Bintanja, R., Fischer, H., Joos, F., Knutti, R., Lohmann, G., and Masson-Delmotte, V.: What caused Earth's temperature variations during the last 800,000 years? Data-based evidence on radiative forcing and constraints on climate sensitivity, Quaternary Sci. Rev., 29, 129–145,, 2010. a

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,, 2004. a

Leloup, G. and Paillard, D.: Influence of the choice of insolation forcing on the results of a conceptual glacial cycle model, Clim. Past, 18, 547–558,, 2022. a

Martinec, Z., Klemann, V., van der Wal, W., Riva, R. E., Spada, G., Sun, Y., Melini, D., Kachuck, S. B., Barletta, V., Simon, K., A, G., and James, T. S.: A benchmark study of numerical implementations of the sea level equation in GIA modelling, Geophys. J. Int., 215, 389–414,, 2018. a

Milankovitch, M.: Kanon der Erdbestrahlung und Seine Andwendung auf das Eiszeitenproblem, vol. 33, Spec. Publ. 132, R. Serbian Acad., Belgrade, 633 pp., 1941. a

Oka, A., Hasumi, H., and Abe-Ouchi, A.: The thermal threshold of the Atlantic meridional overturning circulation and its control by wind stress forcing during glacial climate, Geophys. Res. Lett., 39, 1–6,, 2012. a

Paillard, D.: The timing of Pleistocene glaciations from a simple multiple-state climate model, Nature, 391, 378–381,, 1998. a

Petoukhov, V., Ganopolski, a., Brovkin, V., Claussen, M., Eliseev, A., Kubatzki, C., and Rahmstorf, S.: CLIMBER-2: a climate system model of intermediate complexity. Part I: model description and performance for present climate, Clim. Dynam., 16, 1–17,, 2000. a

Robinson, A., Calov, R., and Ganopolski, A.: An efficient regional energy-moisture balance model for simulation of the Greenland Ice Sheet response to climate change, The Cryosphere, 4, 129–144,, 2010. a, b

Rohling, E. J., Sluijs, A., Dijkstra, H. A., Köhler, P., Van De Wal, R. S., Von Der Heydt, A. S., Beerling, D. J., Berger, A., Bijl, P. K., Crucifix, M., Deconto, R., Drijfhout, S. S., Fedorov, A., Foster, G. L., Ganopolski, A., Hansen, J., Hönisch, B., Hooghiemstra, H., Huber, M., Huybers, P., Knutti, R., Lea, D. W., Lourens, L. J., Lunt, D., Masson-Demotte, V., Medina-Elizalde, M., Otto-Bliesner, B., Pagani, M., Pälike, H., Renssen, H., Royer, D. L., Siddall, M., Valdes, P., Zachos, J. C., and Zeebe, R. E.: Making sense of palaeoclimate sensitivity, Nature, 491, 683–691,, 2012. a

Stouffer, R. J. and Manabe, S.: Equilibrium response of thermohaline circulation to large changes in atmospheric CO2 concentration, Clim. Dynam., 20, 759–773,, 2003. a

Talento, S.: Code and data: New estimation of critical insolation – CO2 relationship for triggering glacial inception, OSF [code and data],, 2023. 

Talento, S. and Ganopolski, A.: Reduced-complexity model for the impact of anthropogenic CO2 emissions on future glacial cycles, Earth Syst. Dynam., 12, 1275–1293,, 2021. a, b

Willeit, M., Ganopolski, A., Calov, R., and Brovkin, V.: Mid-Pleistocene transition in glacial cycles explained by declining CO2 and regolith removal, Science Advances, 5, 1–9,, 2019. a

Willeit, M., Ganopolski, A., Robinson, A., and Edwards, N. R.: The Earth system model CLIMBER-X v1.0 – Part 1: Climate model description and validation​​​​​​​​​​​​​​, Geosci. Model Dev., 15, 5905–5948,, 2022. a, b, c, d, e, f

Willeit, M., Ilyina, T., Liu, B., Heinze, C., Perrette, M., Heinemann, M., Dalmonech, D., Brovkin, V., Munhoven, G., Börker, J., Hartmann, J., Romero-Mujalli, G., and Ganopolski, A.: The Earth system model CLIMBER-X v1.0 – Part 2: The global carbon cycle, Geosci. Model Dev., 16, 3501–3534,, 2023. a, b, c

Willeit, M., Calov, R., Talento, S., Greve, R., Bernales, J., Klemann, V., Bagge, M., and Ganopolski, A.: Glacial inception through rapid ice area increase driven by albedo and vegetation feedbacks, Clim. Past, 20, 597–623,, 2024. a, b, c, d, e, f, g, h, i

Short summary
To trigger glacial inception, the summer maximum insolation at high latitudes in the Northern Hemisphere must be lower than a critical value. This value is not constant but depends on the atmospheric CO2 concentration. Paleoclimatic data do not give enough information to derive the relationship between the critical threshold and CO2. However, knowledge of such a relation is important for predicting future glaciations and the impact anthropogenic CO2 emissions might have on them.