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

HAL
open science

# Influence of the choice of insolation forcing on the results of a conceptual glacial cycle model

Gaëlle Leloup, Didier Paillard

► To cite this version:

Gaëlle Leloup, Didier Paillard. Influence of the choice of insolation forcing on the results of a conceptual glacial cycle model. Climate of the Past, 2022, 18 (3), pp.547-558. 10.5194/cp-18-547-2022. hal-03620611

HAL Id: hal-03620611

https://hal.science/hal-03620611

Submitted on 26 Mar 2022

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_2.png" /></figure>

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

# Influence of the choice of insolation forcing on the results of a conceptual glacial cycle model

Gaëlle Leloup<fnref n="1" />,<fnref n="2" /> and Didier Paillard<fnref n="2" />

**Correspondence:** Gaëlle Leloup (gaelle.leloup@lsce.ipsl.fr)

Received: 6 September 2021 – Discussion started: 13 September 2021

Revised: 11 February 2022 – Accepted: 18 February 2022 – Published: 25 March 2022

**Abstract.** Over the Quaternary, ice volume variations are “paced” by astronomy. However, the precise way in which the astronomical parameters influence the glacial-interglacial cycles is not clear. The origin of the 100 kyr cycles over the last 1 million years and of the switch from 40 to 100 kyr cycles over the Mid-Pleistocene Transition (MPT) remain largely unexplained. By representing the climate system as oscillating between two states, glaciation and deglaciation, switching once glaciation and deglaciation thresholds are crossed, the main features of the ice volume record can be reproduced (Parrenin and Paillard, 2012). However, previous studies have only focused on the use of a single summer insolation as input. Here, we use a simple conceptual model to test and discuss the influence of the use of different summer insolation forcings, having different contributions from precession and obliquity, on the model results. We show that some features are robust. Specifically, to be able to reproduce the frequency shift over the MPT, while having all other model parameters fixed, the deglaciation threshold needs to increase over time, independently of the summer insolation used as input. The quality of the model-data agreement however depends on the chosen type of summer insolation and time period considered.

## 1 Introduction

Palaeoclimate records over the Quaternary (last 2.6 Myr), such as ice cores (Jouzel et al., 2007) or marine cores (Lisiecki and Raymo, 2005) show a succession of oscillations. These oscillations are due to the build-up and retreat of

northern continental ice sheets, corresponding to respectively cold and warm periods known as the glacial-interglacial cycles. Over the last 1 million years, there is an alternation of long glaciations followed by quick deglaciations, leading to cycles of a period of 100 kyr. Yet, in the earliest part of the record, glacial-interglacial cycles are mostly dominated by a frequency of 41 kyr with lower amplitudes. Spectral analysis of the record (Hays et al., 1976) has revealed that astronomical frequencies are imprinted into the ice volume record, suggesting a strong link between insolation and the glacial-interglacial cycles.

The nature and physics of this link has been a central question since the discovery of previous warm and cold periods, and long before the obtention of continuous δ¹⁸O records. It is known that the variations of annual total energy are of too-low an amplitude to explain such changes (Croll, 1864). Therefore, focus has been set on seasonal variations. Croll (1864) assumed that glaciations were linked to colder winters. In contrast, the idea that the decisive element for glaciation was the presence of cold summers, due to reduced summer insolation, at latitudes of the Northern Hemisphere critical for ice sheet growth (65° N), was taken up by Milankovitch, who made it the key element of his ice age theory (Berger, 2021). For conceptual models, this raises the question of which insolation to use as input. When summer insolation is used, this questions the definition of summer: should it be defined as a specific single day, like the summer solstice; the astronomical summer between the two equinoxes; or a fixed number of days around the solstice? The choice leads to very different forcings with different contributions from obliquity and precession. For Earth system models and

<footnote n="1">¹Agence Nationale pour la gestion des déchets radioactifs (ANDRA), 1 Rue Jean Monnet, 92290 Châtenay-Malabry, France</footnote>

<footnote n="2">²Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ-Université Paris-Saclay, 91198 Gif-sur-Yvette, France</footnote>
climate models, insolation is computed at each time step for each grid area, and such choice of the input forcing is not needed. However, other modelling choices have to be made. For instance, several parameterizations are used to represent ice sheet surface melt (Robinson et al., 2010), like the positive degree day (PDD) method (Reeh, 1991), in which surface melt depends solely on air temperature, or the insolation temperature melt (ITM) method (Van den Berg et al., 2008), which takes into account the effect of both temperature and insolation. In both cases, the translation of insolation local and seasonal variations into ice sheet changes and ice age cycles remains an open modelling question.

The oblation of 100 kyr cycles is not possible with a linear theory like that of Milankovitch (Hays et al., 1976; Imbrie and Imbrie, 1980), and a form of non-linearity is needed. Indeed, there is no simple relation between insolation extrema and ice volume extrema. One of the largest deglaciations (termination V) occurred when insolation variations were minimal. Conversely, insolation variations were maximal at termination III, whereas the transition was rather small. In addition, the amplitude of summer insolation variations is maximal every 400 kyr, but this frequency is absent from the palaeoclimatic records. The 100 kyr cycles have been proposed to be linked to either eccentricity-driven variations of precession (Raymo, 1997; Lisiecki, 2010), obliquity (Huybers and Wunsch, 2005; Liu et al., 2008), or both (Huybers, 2011; Parrenin and Paillard, 2012), to internal oscillations phase locked to the astronomical forcing (Saltzman et al., 1984; Paillard, 1998; Gildor and Tziperman, 2000; Tziperman et al., 2006), to internal oscillations independent of the astronomical forcing (Saltzman and Sutera, 1987; Togweiler, 2008) or to period-doubling bifurcation (Verbitsky et al., 2018). Additionally, the Mid-Pleistocene Transition (MPT) and the switch from a 41 kyr dominated record to a 100 kyr one remain mostly unexplained.

Several conceptual models have been developed to try to solve these questions. Calder (1974) was the first to develop a simple conceptual model, linking insolation and ice volume variations. Imbrie and Imbrie (1980) also developed a conceptual model, where the rate of change was different in the case of a warming or cooling climate. Paillard (1998) developed the idea that the glacial-interglacial cycles can be seen as relaxation oscillations between multiple equilibria, like a glaciation and a deglaciation state. It was suggested that the criteria to trigger a deglaciation depends on both insolation and ice volume, whereas insolation alone seems able to trigger glaciation (Parrenin and Paillard, 2003). Here, we adapt and simplify the model of Parrenin and Paillard (2003) and extend it over the whole Quaternary.

One of the critical questions for conceptual models is to decide which insolation to use as input. Milankovitch's work utilized "caloric seasons" at 65° N, the half-year with the highest insolation. This was also the case in Calder's model, which used caloric seasons at 50° N as input. Imbrie and Imbrie (1980) used the July insolation. The use of summer

insolation gradually shifted towards the use of the summer solstice insolation, most probably as it is easier to compute thanks to the tables provided by Berger (1978) and allows one to obtain better fits for the more recent part of the records (Paillard, 2015). More recently, Huybers (2006) suggested that the integrated summer insolation (ISI) over a certain threshold could be better, as it would more closely follow PDDs, an important metric for ice sheet mass balance. Others have also used combinations of orbital parameters as a forcing (Imbrie et al., 2011; Parrenin and Paillard, 2012). However, most models only use one type of insolation forcing and do not consider the influence of other insolation forcing types on the model results.

In our work, we will consider several summer insolation forcings at 65° N (summer solstice insolation, caloric season and ISI over two different threshold values) in order to study their influence on the model results. These different summer insolation forcings have similar shape, but the respective contribution of obliquity and precession differ. The different insolation forcings have different performances in matching the data, depending on the time period considered. However, we show that some features are independent of the insolation forcing used as input. In particular, we are able to reproduce a switch from 41 kyr oscillations before the MPT to 100 kyr cycles afterwards, in agreement with the records for all insolation forcings, by varying a single parameter: the deglaciation threshold $V_0$, and keeping all the other model parameters constant. This is similar to the results of Paillard (1998), who obtained a frequency shift on the glacial cycles by using a linearly increasing deglaciation threshold. This is also coherent with the more recent results of Tzedakis et al. (2017), who demonstrated that the particular sequence of interglacials that happened over the Quaternary and the frequency shift from 41 to 100 kyr could be explained with a simple rule, taking into account a deglaciation threshold that increases over time, leading to "skipped" insolation peaks and longer cycles.

## 2 Conceptual model and methods

### 2.1 Conceptual model

The model used in our study is an adapted and simplified version of the conceptual model of Parrenin and Paillard (2003). For the glacial-interglacial cycles, it is not a new idea that the climate system can be represented by relaxation oscillations between multiple equilibria, like a glaciation and a deglaciation state (Paillard, 1998; Parrenin and Paillard, 2003, 2012). The aim of conceptual models is not to explicitly model and represent physical processes but rather to help us understand critical aspects of the climate system. Here, we do not intend to explicitly represent the numerous physical processes involved in ice sheet volume evolution, affecting surface mass balance, ice discharge to the ocean and bottom melt of grounded ice. Instead, we represent the climate
system with two distinct states of evolution: the “glaciation state” (g) and “deglaciation state” (d). We make the assumption that the evolution of the ice sheet volume in these two states can be simply described by two terms. The first one, common to the glaciation and deglaciation states, is a linear relation to the summer insolation: when the insolation is above average, the ice sheet melts, whereas when the insolation is low enough, the ice sheet grows. The second term, specific to the system state, represents an evolution trend linked to the system state: a slow glaciation in (g) state and a rapid deglaciation in (d) state.

In our model, the evolution of the ice volume in these two states is described by

$$ \left\{ \begin{array}{l} (g) \frac{dv}{dt} = -\frac{I}{\tau_i} + \frac{1}{\tau_g} \\ (d) \frac{dv}{dt} = -\frac{I}{\tau_i} - \frac{v}{\tau_d} \end{array} \right. , \qquad (1) $$

where v represents the normalized ice volume. τi, τd and τg are time constants. I is the normalized summer insolation forcing at 65° N. The implicit assumption is made that the global ice volume changes are mainly driven by the Northern Hemisphere ice sheet waning and waxing, as we focus on the effect of insolation changes at high northern latitudes. This is not a limiting assumption as data suggest much lower sea level variations due to the Antarctic ice sheet than due to those of the Northern Hemisphere over glacial-interglacial timescales. For example, the contribution of the Antarctic ice sheet to the ~120 m Last Glacial Maximum sea level drop is estimated to be between 10 and 35 m sea level equivalent (Lambeck et al., 2014). Other effects like thermal expansion, small glaciers and ice caps are estimated to be around 3 to 4 m sea level equivalent. Furthermore, it has been suggested (Jakob et al., 2020) that the growth of larger Northern Hemisphere ice sheets since the start of the Quaternary and the associated sea level drop has a stabilizing effect on the East Antarctic ice sheet, as it limits its exposure to warm ocean waters.

A critical point is to define the criteria for the switch between the glaciation and deglaciation states. To enter the deglaciation state, both ice volume and insolation seem to play a role (Raymo, 1997; Parrenin and Paillard, 2003, 2012), as terminations occur after considerable build-up of ice sheet over the last 1 million years. To represent the role of both ice volume and insolation in the triggering of deglaciations, the condition to switch from (g) to (d) state uses a linear combination of ice volume and insolation. The deglaciation is triggered when the combination crosses a defined threshold V0: the deglaciation threshold. As in the work of Parrenin and Paillard (2003), this allows transitions to occur with moderate insolation when the ice volume is large enough and reciprocally. Conversely, glacial inception seems to depend on orbital forcing alone (Khodri et al., 2001; Ganopolski and Calov, 2011). Therefore, the condition to switch from the deglaciation state to the glaciation state is based on insolation only: it is possible to enter glaciation when the insolation

becomes low enough.

$$ \left\{ \begin{array}{l} (d) \text{ to } (g) : I < I_0 \\ (g) \text{ to } (d) : I + v > V_0 \end{array} \right. \qquad (2) $$

The idea that the deglaciation threshold is linked to both insolation and ice volume is not new (Parrenin and Paillard, 2003, 2012), and is similar to that developed by Tzedakis et al. (2017), where the threshold for a complete deglaciation decreases with time as the system accumulates instability, with ice sheets becoming more sensitive to insolation increase. As the ice sheet grows and extends to lower latitudes, the insolation needed to reach a negative mass balance decreases. This idea has been confirmed with modelling studies (Abe-Ouchi et al., 2013). Several reasons can explain this increase of instability over time: delayed bedrock rebound and exposure to higher temperature as the ice sheet sinks (Oerlemans, 1980; Pollard, 1982; Abe-Ouchi et al., 2013), increase in basal sliding as the ice sheet grows and the base becomes warmer as it is more isolated from the cold surface temperature (MacAyeal, 1992; Marshall and Clark, 2002), ice sheet margin calving (Pollard, 1982), decrease of the ice sheet albedo due to either snow ageing (Gallée et al., 1992) or increase of dust deposition as the ice sheet expands (Peltier and Marshall, 1995; Ganopolski and Calov, 2011) and increase of atmospheric CO2 due to the release of deep ocean carbon when the Antarctic ice sheet extends over the continental shelves (Paillard and Parrenin, 2004; Bouttes et al., 2012).

## 2.2 Summer insolation

The conceptual model defined previously uses summer insolation as input. It is therefore important to consider which insolation should be used. Insolation is usually taken at 65° N, a typical latitude for Northern Hemisphere ice sheets. Several articles have used the summer solstice daily insolation at 65° N as input insolation forcing (Paillard, 1998; Parrenin and Paillard, 2003). Others (Tzedakis et al., 2017) have used the caloric season, as defined by Milankovitch (1941), defined as the half-year with the highest insolation. It is also possible to use as input a linear combination of orbital parameters (Imbrie et al., 2011; Parrenin and Paillard, 2012). Huybers (2006) also defined the integrated summer insolation (ISI) above a threshold.

Here, for the first time, we want to examine the effect of this choice on the dynamics of a conceptual model. We therefore use four different summer insolation forcings, and compare the model results obtained with each of them. We use the summer solstice insolation, the caloric season and ISI above two different thresholds (300 and 400 W m⁻²). Experiments were also conducted for the summer solstice insolation at 50° N instead of 65° N, but these are not presented here as they do not change the conclusions obtained with the forcings at 65° N. We refer the reader to Sect. S1 in the Supplement for information on these experiments.
These four insolation forcings differ, and have different contributions from precession and obliquity, as can be seen in Fig. 1. For example, the summer solstice insolation has a low contribution from obliquity (41 kyr) compared to the caloric season, or ISI above 300 W m⁻². The summer solstice insolation and the caloric season were computed using the AnalySeries software (Paillard et al., 1996). The ISI for a threshold τ (noted as ISI above τ, or ISI(τ) in this article) was defined by Huybers (2006) as the sum of the insolation on days exceeding this threshold.

$$ \text{ISI}(\tau) = \sum_{i} \beta_{i} (W_{i} \cdot 86400) $$

Here, $W_i$ is the mean insolation of day i in W m⁻² and $\beta_i$ equals one if $W_i > \tau$ and zero otherwise.

To compute ISI(τ), one first needs to compute the daily insolation on day i, and then sum over the year for the days exceeding the threshold. Here, we developed a Python code based on the MATLAB code provided by Huybers (2006). Unlike the insolation files provided by Huybers (2006), which used the solution from Berger and Loutre (1991), we use the orbital parameter value of the Laskar et al. (2004) solution for the calculation. This results in slightly different results for deeper time periods as the calculation of orbital parameters also differ with these two estimations.

## 2.3 Optimal model parameters

The conceptual model relies on a small number of model parameters: $\tau_i$, $\tau_g$, $\tau_d$ and the two thresholds for the switch from a glaciation state to a deglaciation state and the inverse: $V_0$ and $I_0$. For all the simulations performed, we kept $\tau_i = 9$ kyr, $\tau_g = 30$ kyr, $\tau_d = 12$ kyr and $I_0 = 0$, as these parameters gave correct behaviour in previous studies (Parrenin and Paillard, 2003). No attempt was made to tune these parameters to obtain a behaviour closer to the data: on the contrary, focus was set more towards the influence of the deglaciation threshold parameter, $V_0$, on the model results. To compare our model results to the data, we used the benthic $\delta^{18}$O stack “LR04” (Lisiecki and Raymo, 2005) as a proxy for ice volume, considering that most of the $\delta^{18}$O changes of benthic foraminifera represent changes in continental ice (Shackleton, 1967; Shackleton and Opdyke, 1973). Lower $\delta^{18}$O values correspond to lower ice volume. The model results and the LR04 curve were normalized to facilitate their comparison. In the following, “data” refers to the normalized $\delta^{18}$O stack curve LR04.

To estimate which model parameters lead to model results closer to the data, the definition of an objective criteria is needed. The choice of such a criteria is not straightforward, and the use of different criteria could have led to slightly different results. Our model is simple and does not aim at precisely reproducing the ice volume evolution, but rather at reproducing the main qualitative features, such as the shape

and frequency of the oscillations. Therefore, we used a criteria based on the state of the system: glaciation or deglaciation. Similar results can be obtained using a simple correlation coefficient (see Sect. S2). The definition of the deglaciation state in the data is explained in Sect. 2.4. A critical point is that our model should be able to deglaciate at the right time, when a deglaciation is seen in the data. Conversely, the model should not produce deglaciation at periods where the data do not show deglaciation.

We defined a criteria for each of these two conditions, and assembled them in a global criteria. To determine if deglaciation is well reproduced by our model, we look at the state of the model (glaciation or deglaciation) at the time halfway between the start of the deglaciation and the end of the deglaciation. If the model state at that time is deglaciation, the deglaciation is considered as correctly reproduced. Otherwise, it is considered as a “missed” deglaciation. We simply defined the criteria $c_1$ as the fraction of deglaciations correctly reproduced. This criteria equals one when all the deglaciations take place at the right time, meaning that the model produces a deglaciation state every time deglaciation is seen in the data. To ensure that the model does not deglaciate too often, we looked at insolation maxima that are not associated with deglaciation in the data and ensured the corresponding model state was glaciation. We decided to look specifically at insolation maxima, as they are where the model is the most likely to deglaciate when it should not. A deglaciation seen in the model at a place where no deglaciation is seen in the data is considered as a “wrong transition”. We defined the criteria $c_2$ as $c_2 = 1 - w$, with w the fraction of wrong transitions. The $c_2$ criteria equals one when the model does not deglaciate when it should not. To take into account these two conditions, the overall accuracy criteria c is defined as $c = c_1 \cdot c_2$. It is equal to one when all the deglaciations are correctly placed and when no additional deglaciation compared to the data take place.

In order to study the evolution of the optimal deglaciation threshold $V_0$ over the Quaternary, it was divided into five 500 kyr periods. The $V_0$ values that maximize the accuracy criteria for each time period and insolation forcing are called “optimal $V_0$”. To determine the optimal $V_0$ threshold corresponding to each period and insolation forcing, several simulations were carried out and the parameter values maximizing the accuracy criteria c were selected. More precisely, for each insolation and period, 3500 simulations corresponding to different $V_0$ thresholds (from $V_0 = 1.0$ to $V_0 = 8.0$ with a step of 0.1) and different initial conditions (initial volume $V_{\text{init}}$ ranging from 0.0 to 5.0 with a step of 0.2, and the two possible initial states – glaciation or deglaciation) were performed. The numerical integration of the model equations was done with a fourth-order Runge-Kutta scheme.

For each insolation forcing, the best fit over the Quaternary is defined as the simulation over the whole Quaternary (0 to 2500 ka) with $V_0$ changing with time that is equal to the corresponding optimal $V_0$ at each time period. Addition-
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) The four different summer insolation forcings at 65° N (summer solstice, caloric season, ISI above 300 W m⁻² and ISI above 400 W m⁻²) displayed over the Quaternary. The insolation forcings are centred and normalized by their standard deviation. (b) Corresponding spectral analysis over the Quaternary, normalized by standard deviation.</figcaption></figure>

ally, simulations were performed to determine the optimal $V_0$ threshold obtained when the optimization procedure is carried out over the whole Quaternary. It is called $V_0^Q$ in the following.

## 2.4 Definition of the deglaciation state in the data

To calculate our accuracy criteria *c* and therefore determine the optimal $V_0$ threshold over a given period, a definition of the deglaciation in the data is needed. We based our definition on the interglacial classification developed by Tzedakis et al. (2017). Tzedakis et al. (2017) differentiates interglacials, continued interglacials and interstadials based on a detrended version of the LR04 stack curve. A period is considered as an interglacial if its isotopic $\delta^{18}$O is below a threshold (3.68 ‰). Two interglacials are considered to be separated if they are separated by a local maximum above a threshold (3.92 ‰). This definition differs from the usual characterization of terminations, sometimes leading to several interglacials in the same isotopic stage. The definition of Tzedakis et al. (2017) is for interglacials, and as our focus is not on interglacial periods but rather on deglaciations, we adapted it. We defined deglaciations as periods of decreasing $\delta^{18}$O (and thus, ice volume, in our assumptions) preceding the interglacials. The start of the deglaciation is taken as the last local maxima above the threshold of 3.92 ‰, and the end of the deglaciation

is taken as the first local minima below the interglacial threshold of 3.68 ‰. The deglaciation periods in the data correspond to the time between the deglaciation starts and deglaciation ends are displayed with blue shading in Figs. 4 and 5.

# 3 Results and discussion

## 3.1 Optimal deglaciation threshold $V_0$ and corresponding accuracy

For each insolation, we computed the deglaciation threshold $V_0$ maximizing the accuracy for each of the five 500 kyr periods and the fixed $V_0^Q$ value maximizing the accuracy over the whole Quaternary. The results are displayed in Fig. 2.

In some cases, several values of the deglaciation threshold $V_0$ lead to the same accuracy criteria, whereas in others there is only one $V_0$ value maximizing the accuracy. When several $V_0$ values are equivalent, the mean value was plotted and the other possible values are represented with error bars.

Over the same time period, different insolation values lead to slightly different optimal $V_0$ thresholds. But the most striking feature is the increase of the optimal $V_0$ threshold over time, and more specifically for the last 1 million years, after the MPT. This feature is valid regardless of the insolation type. The optimal $V_0$ over the whole Quaternary $V_0^Q$ is between 3.4 and 4 for each insolation type. It is a value in the
<figure><img src="image_5.png" /></figure>

**Figure 2.** Optimal deglaciation threshold $V_0$ over the five 500 kyr periods for the four different summer insolation forcings at 65° N and optimal constant $V_0$ threshold obtained when the optimization procedure is done over the whole Quaternary ($V_0^Q$). When several values of the deglaciation threshold $V_0$ maximize the accuracy criteria, the mean value is plotted and the other possible values are represented with error bars.

middle of the highest values that best fit the latter part of the record and the lowest values that best fit the earliest part of the record.

For each insolation, the accuracy corresponding to the optimal $V_0$ threshold for each time period as well as to the fixed $V_0^Q$ value maximizing the accuracy over the whole Quaternary is displayed in Fig. 3.

It is first noticeable that the accuracy is higher over the last 1 million-year period, regardless the input summer insolation used. In the last 1 million years, the summer solstice insolation as input produces the best results. However, this is no longer the case for older time periods: the solstice insolation gives the worst results at the start of the Quaternary. The accuracy obtained for the whole Quaternary period (fixed $V_0^Q$ value) is globally lower than the accuracy for each time period. This is due to the fact that the $V_0^Q$ values obtained are lower than the optimal $V_0$ values in the later part of the Quaternary and higher than the optimal $V_0$ values in the earliest part of the Quaternary, leading to a poorer representation of both of these periods.

In the earlier part of the Quaternary (periods earlier than 1.5 Ma), the results are less robust. This is due to increased uncertainties in the LR04 record, and the associated definition of interglacial periods, which affects our accuracy criteria. In their classification, Tzedakis et al. (2017) stated that the definition of the valley depth needed to separate several interglacials was quite straightforward in the earlier part of

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

**Figure 3.** Accuracy over the five 500 kyr periods for the four different summer insolation forcings at 65° N, and accuracy criteria over the Quaternary when the optimization procedure is done over the whole Quaternary (constant $V_0^Q$).

the record, whereas for earlier time periods (ages older than 1.5 Ma), their method led to several “borderline cases”. Additionally, slightly different choices for the interglacial and interstadial thresholds would have led to a different population of interglacial. For us, this would lead to a different definition of deglaciation, and thus a different accuracy criteria. In addition, the resolution of the LR04 curve decreases with increasing age (Lisiecki and Raymo, 2005), and for periods with lower resolution or more uncertain age matching, the amplitude of the peaks might be reduced (Tzedakis et al., 2017).

Moreover, the $\delta^{18}$O LR04 curve includes at the same time an ice volume and deep water temperature component. Ice volume and sea level reconstructions do exist (Bintanja et al., 2005; Spratt and Lisiecki, 2016), but are however limited to the more recent part of the Quaternary and do not allow for the investigation of the pre-MPT period. The use of $\delta^{18}$O as an ice volume proxy has already been largely debated (Shackleton, 1967; Chappell and Shackleton, 1986; Shackleton and Opdyke, 1973; Clark et al., 2006), and recent studies (Elderfield et al., 2012) have shown that the temperature component may be as large as 50 %. Furthermore, the stack was tuned to insolation (Lisiecki and Raymo, 2005). We refer the reader to Raymo et al. (2018) for a review of possible biases in the interpretation of the LR04 benthic $\delta^{18}$O stack as an ice volume and sea level reconstruction. All these reasons encourage us to remain at a qualitative level to fit the data.

The increase in the optimal deglaciation threshold $V_0$ over the MPT is however a robust feature that does not depend on the input summer insolation forcing used. This is coherent with the work of Paillard (1998) and Tzedakis et al. (2017),
who obtained a frequency shift from 41 to 100 kyr cycles with an increasing deglaciation threshold. The classical hypothesis and possible physical meanings of the rise of the $V_0$ threshold often imply a gradual cooling, due to lower $pCO_2$ (Raymo, 1997; Paillard, 1998; Berger, 1978), but this idea is not supported by more recent data, which do not show a gradual decrease in the $pCO_2$ trend, but rather only a decrease in minimal $pCO_2$ values (Hönisch et al., 2009; Yan et al., 2019). Other hypothesis imply a gradual erosion of a thick regolith layer that exposed the ice sheet to a higher friction bedrock, allowing for a thicker ice sheet to develop (Clark and Pollard, 1998; Clark et al., 2006), or a sea ice switching mechanism (Gildor and Tziperman, 2000).

The overall lowest accuracy in the older part of the record suggests that our non-linear threshold model is less adapted for this time period. Indeed, the ice volume might respond more linearly to the insolation forcing before the MPT, as some studies suggest (Tziperman and Gildor, 2003; Raymo and Nisancioglu, 2003). In contrast, others (Ashkenazy and Tziperman, 2004) have shown that the 41 kyr pre-MPT oscillations are in fact significantly asymmetric and therefore suggested that oscillations both before and after the MPT could be explained as the self-sustained variability of the climate system, phase locked to the astronomical forcing. In our model, the deglaciation threshold $V_0$ changes over time, leading to an amplitude change of the cycles. However, it cannot be excluded that the mechanisms behind the pre- and post-MPT glacial cycles are structurally different, and that they cannot be explained with the same physics and the same equations in our model. Moreover, this kind of conceptual model has mainly been constructed in order to explain the non-linear features of the 100 kyr cycles, and it is therefore not surprising that its agreement with the data on the pre-MPT period is less satisfying.

## 3.2 Best fit over the Quaternary

Our conceptual model is able to reproduce qualitatively well the data (LR04 normalized curve) over the whole Quaternary. The model's best fit over the Quaternary for each insolation forcing, as defined in Sect. 2.3, is displayed in Fig. 4. It is able to reproduce the frequency shift from a dominant 41 kyr period before 1 Ma to longer cycles afterwards, as observed in the data, and thus by varying only one parameter during the whole simulation length: the deglaciation threshold $V_0$. Like the previous model of Parrenin and Paillard (2003), from which our model is adapted, ours quickly looses its sensitivity to the initial conditions (after no more than 200 kyr). For every input forcing, longer cycles are obtained in the last part of the Quaternary (last Myr). Figure 4 displays the results over the whole Quaternary with the $V_0$ threshold being set to its optimal value for each 500 kyr period, while Fig. 5 displays the results over the last 1 million years with the $V_0$ threshold being set to its optimal value for the [0–1000] ka period.

For the last 1 million years, it is possible to reproduce all terminations with the right timing, apart from the last deglaciation, for all insolation forcings, by using a single value of the $V_0$ threshold over this period. Some differences are however noticeable between the different forcings. In particular the agreement for the ISI above $300 \text{ W m}^{-2}$ forcing is not as good as for the other forcings: termination V (around 420 ka) is triggered later compared to in the data, while termination III (around 240 ka) is triggered too early. For the ISI above $300 \text{ W m}^{-2}$ forcing, the range of $V_0$ values that allow us to correctly reproduce most of the terminations during the last 1 million years is reduced (only values of $V_0 = 3.9–4.0$), whereas the results are more robust for the three other insolation forcings, with a broader range of working $V_0$ values. The ISI above $300 \text{ W m}^{-2}$ forcing has a low precession component, which explains why it is less successful in reproducing the data over the last 1 million years. Experiments with our model setup have shown that a summer forcing with no precession component could not successfully reproduce the data over the post-MPT period as accurately as the four forcings presented here that contain both precession and obliquity (see Sect. S3).

Despite the accurate timing of terminations, the spectral analysis of the model results over the last 1 million years differs from the spectral analysis of the data. For all forcings except the summer solstice insolation, obliquity continues to dominate after the MPT. The spectral analysis shows secondary and third peaks of lower frequency, but does not show a sharp 100 kyr cyclicity as in the LR04 record. Compared to the data, all the model outputs over the post-MPT period have a more pronounced obliquity and precession component and a less pronounced 100 kyr component. This feature is most probably due to the model formulation, and more specifically the direct dependence of ice volume evolution to insolation via the $dV/dt = -I/\tau_1$ term. This is one of the limits of our conceptual model. While the criteria on the switch to deglaciations allow us to reproduce the deglaciations at the right timing, the direct dependence of ice volume change to the insolation forcing is definitely too simplistic and probably produces an overestimated dependency of the ice evolution to the astronomical forcing for the latter part of the record.

In the first part of the Quaternary (2.6 to 1 Ma), the spectral analysis of the data is dominated by a 41 kyr (obliquity) peak. It is also the case for the model results, for each type of insolation. However, the model outputs also show a precession component (19 to 23 kyr), especially for the summer solstice and the ISI above $400 \text{ W m}^{-2}$ forcings, which does not exist on the data. Although our model reproduces the features of the record qualitatively well, there are some noticeable model-data mismatches that occur for all insolation types around 1100 and 2030 ka. Additionally, the amplitude of the oscillations produced is slightly too high in the older part of the Quaternary, especially when the ISI above $300 \text{ W m}^{-2}$ is used.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Best model fit over the whole Quaternary and corresponding spectral analysis. (b) presents the best fit of the model for the different summer insolation used as input, compared to the data. The results for the summer solstice insolation, caloric season, ISI above 300 and ISI above 400 W m⁻² are displayed. The data (normalized LR04 curve) are in black. The blue shading presents deglaciation periods in the data and the yellow shading deglaciation periods in the model. This results in a green shading when deglaciations are seen in the model and data at the same time. (a) presents the spectral analysis of the best fit solution over the last 1 million years. (c) presents the spectral analysis over the older part of the Quaternary (before 1 Ma).</figcaption></figure>

On the oldest part of the Quaternary, the caloric season forcing and the ISI above 300 W m⁻² seem to perform better in reproducing the frequency of the oscillations than the summer solstice insolation and ISI above 400 W m⁻². They correspond to the forcing with the strongest obliquity (41 kyr) component, which might explain why they more successfully represent this part of the record, which is dominated by obliquity.

Over the last 1 million years, the highest accuracy is obtained with the summer solstice insolation as input forcing ($c = 0.92$ for the summer solstice insolation, $c = 0.82$ for the caloric season and ISI above 400 W m⁻² and $c = 0.87$ for the ISI above 300 W m⁻²). This is mainly due to the fact that for the three other forcings (caloric season, ISI above 300 and ISI above 400 W m⁻²), the last deglaciation occurs one insolation peak too early, around 50 ka. However, if one computes the accuracy over the last 1 million years excluding the last 100 kyr (period from 100 ka to 1000 ka), the accuracy is similar for all four forcings (0.92 for summer solstice, caloric season and ISI above 400 W m⁻² and 0.89 for ISI above 300 W m⁻²). Figure 5 displays two different alternatives. For the first one (full line), the optimal $V_0$ is calibrated over the [100–1000] ka period and maintained for the whole simulation ($V_0 = 5.0$ for the summer solstice insolation, 4.5 for the caloric season, 4.0 for the ISI above 300 W m⁻² and 4.6 for the ISI above 400 W m⁻²). Except with the summer solstice insolation as input, the three other insolation forcings

fail to accurately reproduce the last deglaciation, as the last deglaciation occurs one insolation peak too early. The second alternative (dashed line) is to raise the deglaciation threshold $V_0$ over the last cycle (raised to $V_0 = 5.5$ from 100 ka onward). In this case, the model does accurately reproduce the last deglaciation for all insolation forcings. This suggests that we might need to further raise the $V_0$ threshold in order to model a theoretical “natural” future (without anthropogenic influence) with longer cycles. A cycle is not enough to conclude about a trend, but this should be envisaged for future natural scenarios. This is coherent with the idea of Paillard (1998), who used a linearly increasing deglaciation threshold over the Quaternary.

## 3.3 Reflections about the future

To model future natural evolutions of the climate system, possible evolutions of the $V_0$ threshold should be considered. However, we do not exclude the fact that variations of other parameters, which were kept constant in this study, could vary in the future. For instance, different $I_0$ thresholds have to be considered. The fate of the next thousands of years in a natural scenario cannot be ascertained with our conceptual model. Indeed, this depends on the onset (or not) of a deglaciation, which is determined by the $I_0$ parameter. In our model, a broad range of $I_0$ parameter values lead to satisfying results over the last 1 million years. These different values pro-
<figure><img src="image_8.png" /></figure>

**Figure 5.** Normalized model results over the last 1 million years, with the different summer insolation forcings: insolation at the summer solstice, caloric season, ISI above 300 W m⁻² and ISI above 400 W m⁻². The full coloured lines are the best fit computed over the 100–1000 ka period for each input insolation. The dashed lines represent the same solution, but with an increased V₀ threshold for the last deglaciation. The data (normalized LR04 stack curve) are in black. The blue shading represents deglaciation periods in the data. The yellow shading represents deglaciation periods in the model outputs (case of the increased V₀ threshold). This results in a green shading when deglaciations are seen in the model and data at the same time.

jected into the future lead either to, in one case, a glaciation start at present, or, in the other case, a continued deglaciation, switching to a glaciation state only around 50 kyr after present at the next insolation minima. This is due to the current particular astronomical configuration, with very low eccentricity, which leads to high summer insolation minima, where the threshold for glaciation might not be reached. This is coherent with previous studies (Berger and Loutre, 2002; Paillard, 2001), which suggested that without anthropogenic forcing, the present interglacial might have lasted 50 kyr.

However, this exercise is purely academic, as we are not taking into account the role of anthropogenic CO₂, which would affect the glaciation and deglaciation thresholds (Archer and Ganopolski, 2005; Talento and Ganopolski, 2021). Furthermore, our conceptual model cannot be extended outside the Quaternary, as the ice volume variations considered are exclusively those of the Northern Hemisphere, and our model is, by construction, unable to represent projected future Antarctic ice sheet mass loss.

## 4 Conclusions

We have used a conceptual model with very few tunable parameters that represents the climatic system with multiple equilibria and relaxation oscillation. Only one parameter was varied, the deglaciation threshold parameter V₀. We used different summer insolation as input for our conceptual model: the summer solstice insolation, the caloric season and ISI

over two different thresholds. With all these forcings, which have different contributions from obliquity and precession, we are able to reproduce the features of the ice volume over the whole Quaternary. More specifically, we are able to represent the MPT and the switch from a 41 kyr dominated record to larger cycles by raising the deglaciation threshold and keeping the other model parameters constant. This rise in the deglaciation threshold is valid regardless of the type of summer insolation forcing used as input. However, the data agreement is less satisfying before the MPT. This suggests the possibility that climate mechanisms might be structurally different before and after the MPT, with a more linear behaviour in pre-MPT conditions. This highlights that models are designed to answer rather specific questions, and a model built specifically to explain 100 kyr cycles might be less efficient in a more linear setting. More generally, this kind of glacial-interglacial conceptual model is designed to explain the main features of the Quaternary time period characterized by the waning and waxing of Northern Hemisphere ice sheets under the influence of changing astronomical parameters. In our case, this raises the question of which physical phenomena are responsible for making deglaciations “harder” to start in the latest part of the Quaternary compared to the earliest part. This kind of model is however unlikely to be directly applicable in a more general context, like during the Pliocene and earlier periods, or in the context of future climates under the long-term persistence of anthropogenic CO₂ (Archer and Ganopolski, 2005; Talento and Ganopolski, 2021). In order
to tackle such questions, it would be critical to gain a deeper
understanding of the natural evolution of the other forcings
involved in the climate system and most notably of the dy-
namics of the carbon cycle. Conceptual models are likely to
pave the way in this direction (Paillard, 2017). Indeed, just
as in the case of the Quaternary, a full mechanistic simula-
tion of the many processes at work is currently out of reach
and modelling work can only be very exploratory. Here, we
have shown that some robust features are required to explain
Quaternary ice age cycles. Similar conceptual modelling on
a wider temporal scope over the Cenozoic could help better
explain the connections between astronomical forcing, the
carbon cycle, ice sheets and the climate. This would help us
imagine what the Anthropocene might be like.

**Code and data availability.** The model code, insolation input files, spectral analysis and code needed to reproduce the figures are available for download: https://doi.org/10.5281/zenodo.6045532 (Leloup, 2022).

**Supplement.** The supplement related to this article is available online at: https://doi.org/10.5194/cp-18-547-2022-supplement.

**Author contributions.** GL and DP designed the study. GL performed the simulations, and also wrote the manuscript under the supervision of DP.

**Competing interests.** The contact author has declared that neither they nor their co-author has any competing interests.

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

**Special issue statement.** This article is part of the special issue "A century of Milankovic's theory of climate changes: achievements and challenges (NPG/CP inter-journal SI)". It is not associated with a conference.

**Acknowledgements.** We acknowledge the use of the LSCE storage and computing facilities. We also thank the reviewers for their helpful comments and suggestions.

**Financial support.** This research has been supported by ANDRA (contract no. 20080970).

**Review statement.** This paper was edited by Marie-France Loutre and reviewed by Mikhail Verbitsky, Andrey Ganopolski, and one anonymous referee.

References

Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M., 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.

Archer, D. and Ganopolski, A.: A movable trigger: Fossil fuel CO₂ and the onset of the next glaciation, Geochem. Geophy. Geosy., 6, Q05003, https://doi.org/10.1029/2004GC000891, 2005.

Ashkenazy, Y. and Tziperman, E.: Are the 41 kyr glacial oscillations a linear response to Milankovitch forcing?, Quaternary Sci. Rev., 23, 1879–1890, https://doi.org/10.1016/j.quascirev.2004.04.008, 2004.

Berger, A.: Long-term variations of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362–2367, 1978.

Berger, A.: Milankovitch, the father of paleoclimate modeling, Clim. Past, 17, 1727–1733, https://doi.org/10.5194/cp-17-1727-2021, 2021.

Berger, A. and Loutre, M.-F.: Insolation values for the climate of the last 10 million years, Quaternary Sci. Rev., 10, 297–317, https://doi.org/10.1016/0277-3791(91)90033-Q, 1991.

Berger, A. and Loutre, M.-F.: Climate: An Exceptionally Long Interglacial Ahead?, Science, 297, 1287–1288, https://doi.org/10.1126/science.1076120, 2002.

Bintanja, R., van de Wal, R., and Oerlemans, J.: Modelled atmospheric temperatures and global sea levels over the past million years, Nature, 437, 125–128, https://doi.org/10.1038/nature03975, 2005.

Bouttes, N., Paillard, D., Roche, D. M., Waelbroeck, C., Kageyama, M., Lourantou, A., Michel, E., and Bopp, L.: Impact of oceanic processes on the carbon cycle during the last termination, Clim. Past, 8, 149–170, https://doi.org/10.5194/cp-8-149-2012, 2012.

Calder, N.: Arithmetic of ice ages, Nature, 252, 216–218, https://doi.org/10.1038/252216a0, 1974.

Chappell, J. and Shackleton, N. J.: Oxygen isotopes and sea level, Nature, 324, 137–140, 1986.

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

Clark, P., Archer, D., Pollard, D., Blum, J. D., Rial, J., Brovkin, V., Mix, A., Pisias, N., and Roy, M.: The middle Pleistocene transition: characteristics, mechanisms, and implications for long-term changes in atmospheric pCO₂, Quaternary Sci. Rev., 25, 3150–3184, https://doi.org/10.1016/j.quascirev.2006.07.008, 2006.

Croll, J.: XIII. On the physical cause of the change of climate during geological epochs, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 28, 121–137, https://doi.org/10.1080/14786446408643733, 1864.

Elderfield, H., Ferretti, P., Greaves, M., Crowhurst, S., Mc-cave, I. N., Hodell, D., and Piotrowski, A. M.: Evolu-tion of Ocean Temperature and Ice Volume Through the Mid-Pleistocene Climate Transition, Science, 337, 704–709, https://doi.org/10.1126/science.1221294, 2012.

Gallée, H., Van Yperselb, J. P., Fichefet, T., Marsiat, I., Tricot, C., and Berger, A.: Simulation of the last glacial cycle by a cou-
pled, sectorially averaged climate-ice sheet model: 2. Response to insolation and CO₂ variations, J. Geophys. Res.-Atmos., 97, 15713–15740, https://doi.org/10.1029/92JD01256, 1992.

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.

Gildor, H. and Tziperman, E.: Sea ice as the glacial cycles’ Climate switch: role of seasonal and orbital forcing, Paleooceanography, 15, 605–615, https://doi.org/10.1029/1999PA000461, 2000.

Hays, J., Imbrie, J., and Shackleton, N. J.: Variations in the Earth’s Orbit: Pacemaker of the Ice Ages, Science, 194, 1121–1132, https://doi.org/10.1126/science.194.4270.1121, 1976.

Hönisch, B., Hemming, N. G., Archer, D., Siddall, M., and McManus, J. F.: Atmospheric Carbon Dioxide Concentration Across the Mid-Pleistocene Transition, Science, 324, 1551–1554, https://doi.org/10.1126/science.1171477, 2009.

Huybers, P.: Early Pleistocene glacial cycles and the integrated summer insolation forcing, Science, 313, 508–511, https://doi.org/10.1126/science.1125249, 2006.

Huybers, P.: Combined obliquity and precession pacing of late Pleistocene deglaciations, Nature, 480, 229–232, https://doi.org/10.1038/nature10626, 2011.

Huybers, P. and Wunsch, C.: Obliquity pacing of the late Pleistocene glacial terminations, Nature, 434, 491–494, https://doi.org/10.1038/nature03401, 2005.

Imbrie, J. and Imbrie, J. Z.: Modeling the Climatic Response to Orbital Variations, Science, 207, 943–953, https://doi.org/10.1126/science.207.4434.943, 1980.

Imbrie, J. Z., Imbrie-Moore, A., and Lisiecki, L. E.: A phase-space model for Pleistocene ice volume, Earth Planet. Sci. Lett., 307, 94–102, 2011.

Jakob, K. A., Wilson, P. A., Pross, J., Ezard, T. H. G., Fiebig, J., Repschläger, J., and Friedrich, O.: A new sea-level record for the Neogene/Quaternary boundary reveals transition to a more stable East Antarctic Ice Sheet, P. Natl. Acad. Sci. USA, 117, 30980–30987, https://doi.org/10.1073/pnas.2004209117, 2020.

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J., Chappel-laz, J., Fischer, H., Gallet, J.-C., Johnsen, S., Leuenberger, M., Louergue, L., Lüthi, D., Oerter, H., Parrenin, F., Raisbeck, G., and Wolff, E.: Orbital and Millennial Antarctic Climate Variability over the Past 800,000 Years, Science, 317, 793–796, https://doi.org/10.1126/science.1141038, 2007.

Khodri, M., Leclainche, Y., Ramstein, G., Braconnot, P., Marti, O., and Cortijo, E.: Simulating the amplification of orbital forcing by ocean feedbacks in the last glaciation, Nature, 410, 570–574, 2001.

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., and Levrard, B.: A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285, 2004.

Leloup, G.: Data and code for results and figures in Climate of the Past “Influence of the choice of insolation forcing on the results of a conceptual glacial cycle model”, Zenodo [code and data set], https://doi.org/10.5281/zenodo.6045532, 2022

Lisiecki, L. E.: Links between eccentricity forcing and the 100,000-year glacial cycle, Nat. Geosci., 3, 349–352, https://doi.org/10.1038/ngeo828, 2010.

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic δ¹⁸O records, Paleooceanography, 20, PA1003, 2005.

Liu, Z., Cleaveland, L., and Herbert, T.: Early onset and origin of 100-kyr cycles in Pleistocene tropical SST records, Earth Planet. Sci. Lett., 265, 703–715, https://doi.org/10.1016/j.epsl.2007.11.016, 2008.

MacAyeal, D. R.: Irregular oscillations of the West Antarctic ice sheet, Nature, 359, 29–32, https://doi.org/10.1038/359029a0, 1992.

Marshall, S. J. and Clark, P. U.: Basal temperature evolution of North American ice sheets and implications for the 100-kyr cycle, Geophys. Res. Lett., 29, 67-1-67-4, https://doi.org/10.1029/2002GL015192, 2002.

Milankovitch, M.: Kanon der Erdbestrahlung und seine Anwendung auf das Eiszeitenproblem, Royal Serbian Academy Special Publication, 133, 1–633, 1941.

Oerlemans, J.: Model experiments on the 100,000-yr glacial cycle, Nature, 287, 430–432, https://doi.org/10.1038/287430a0, 1980.

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

Paillard, D.: Glacial cycles: toward a new paradigm, Rev. Geophys., 39, 325–346, https://doi.org/10.1029/2000RG000091, 2001.

Paillard, D.: Quaternary glaciations: from observations to theories, Quaternary Sci. Rev., 107, 11–24, https://doi.org/10.1016/j.quascirev.2014.10.002, 2015.

Paillard, D.: The Plio-Pleistocene climatic evolution as a consequence of orbital forcing on the carbon cycle, Clim. Past, 13, 1259–1267, https://doi.org/10.5194/cp-13-1259-2017, 2017.

Paillard, D. and Parrenin, F.: The Antarctic ice sheet and the triggering of deglaciations, Earth Planet. Sci. Lett, 227, 263–271, https://doi.org/10.1016/j.epsl.2004.08.023, 2004.

Paillard, D., Labeyrie, L., and Yiou, P.: Macintosh Program performs time-series analysis, Eos Transactions, AGU, 77, 379–379, https://doi.org/10.1029/96EO00259, 1996.

Parrenin, F. and Paillard, D.: Amplitude and phase of glacial cycles from a conceptual model, Earth Planet. Sci. Lett., 214, 243–250, https://doi.org/10.1016/S0012-821X(03)00363-7, 2003.

Parrenin, F. and Paillard, D.: Terminations VI and VIII (~530 and ~720 kyr BP) tell us the importance of obliquity and precession in the triggering of deglaciations, Clim. Past, 8, 2031–2037, https://doi.org/10.5194/cp-8-2031-2012, 2012.

Peltier, W. R. and Marshall, S.: Coupled energy-balance/ice-sheet model simulations of the glacial cycle: A possible connection between terminations and terrigenous dust, J. Geophys. Res.-Atmos., 100, 14269–14289, https://doi.org/10.1029/95JD00015, 1995.

Pollard, D.: A simple ice sheet model yields realistic 100 kyr glacial cycles, Nature, 296, 334–338, https://doi.org/10.1038/296334a0, 1982.

Raymo, M.: The timing of major climate terminations, Paleooceanography, 12, 577–585, https://doi.org/10.1029/97PA01169, 1997.
Raymo, M. E. and Nisancioglu, K. H.: The 41 kyr world: Milankovitch's other unsolved mystery, Paleooceanography, 18, 1011, https://doi.org/10.1029/2002PA000791, 2003.

Raymo, M. E., Kozdon, R., David, E., Lisiecki, L., and Ford, H. L.: The accuracy of mid-Pliocene $\delta^{18}$O-based ice volume and sea level reconstructions, Earth-Sci. Rev., 177, 291–302, https://doi.org/10.1016/j.earscirev.2017.11.022, 2018.

Reeh, N.: Parameterization of Melt Rate and Surface Temperature on the Greenland Ice Sheet, Polarforschung, 59, 113–128, 1991.

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, https://doi.org/10.5194/tc-4-129-2010, 2010.

Saltzman, B. and Sutera, A.: The Mid-Quaternary Climatic Transition as the Free Response of a Three-Variable Dynamical Model, J. Atmos. Sci., 44, 236–241, https://doi.org/10.1175/1520-0469(1987)044<0236:TMQCTA>2.0.CO;2, 1987.

Saltzman, B., Hansen, A., and Maasch, K.: The Late Quaternary Glaciations as the Response of a Three-Component Feedback System to Earth-Orbital Forcing, J. Atmos. Sci., 41, 3380–3389, https://doi.org/10.1175/1520-0469(1984)041<3380:TLQGAT>2.0.CO;2, 1984.

Shackleton, N.: Oxygen Isotope Analyses and Pleistocene Temperatures Re-assessed, Nature, 215, 15–17, https://doi.org/10.1038/215015a0, 1967.

Shackleton, N. J. and Opdyke, N. D.: Oxygen isotope and palaeomagnetic stratigraphy of Equatorial Pacific core V28-238: Oxygen isotope temperatures and ice volumes on a 105 year and 106 year scale, Quaternary Res., 3, 39–55, https://doi.org/10.1016/0033-5894(73)90052-5, 1973.

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.

Talento, S. and Ganopolski, A.: Reduced-complexity model for the impact of anthropogenic $CO_2$ emissions on future glacial cycles, Earth Syst. Dynam., 12, 1275–1293, https://doi.org/10.5194/esd-12-1275-2021, 2021.

Toggweiler, J. R.: Origin of the 100,000-year timescale in Antarctic temperatures and atmospheric $CO_2$, Paleooceanography, 23, PA2211, https://doi.org/10.1029/2006PA001405, 2008.

Tzedakis, P., Crucifix, M., Mitsui, T., and Wolff, E. W.: A simple rule to determine which insolation cycles lead to interglacials, Nature, 542, 427–432, https://doi.org/10.1038/nature21364, 2017.

Tziperman, E. and Gildor, H.: On the mid-Pleistocene transition to 100-kyr glacial cycles and the asymmetry between glaciation and deglaciation times, Paleooceanography, 18, 1–1–1–8, https://doi.org/10.1029/2001pa000627, 2003.

Tziperman, E., Raymo, M. E., Huybers, P., and Wunsch, C.: Consequences of pacing the Pleistocene 100 kyr ice ages by nonlinear phase locking to Milankovitch forcing, Paleooceanography, 21, PA4206, https://doi.org/10.1029/2005PA001241, 2006.

Van den Berg, J., van de Wal, R., and Oerlemans, H.: A mass balance model for the Eurasian Ice Sheet for the last 120,000 years, Global Planet. Change, 61, 194–208, https://doi.org/10.1016/j.gloplacha.2007.08.015, 2008.

Verbitsky, M. Y., Crucifix, M., and Volobuev, D. M.: A theory of Pleistocene glacial rhythmicity, Earth Syst. Dynam., 9, 1025–1043, https://doi.org/10.5194/esd-9-1025-2018, 2018.

Yan, Y., Bender, M. L., Brook, E. J., Clifford, H. M., Kemeny, P. C., Kurbatov, A. V., Mackay, S., Mayewski, P. A., Ng, J., Severinghaus, J. P., and Higgins, J. A.: Two-million-year-old snapshots of atmospheric gases from Antarctic ice, Nature, 574, 663–666, https://doi.org/10.1038/s41586-019-1692-3, 2019.