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

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 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 obtention 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; Toggweiler, 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₀, 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

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}\left(g\right)\frac{\mathrm{d}v}{\mathrm{d}t}=-\frac{I}{{\mathit{\tau }}_{\mathrm{i}}}+\frac{\mathrm{1}}{{\mathit{\tau }}_g}\\ \left(d\right)\frac{\mathrm{d}v}{\mathrm{d}t}=-\frac{I}{{\mathit{\tau }}_{\mathrm{i}}}-\frac{v}{{\mathit{\tau }}_d}\end{array},\right.\end{array}$$

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 V₀: 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.

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}\left(d\right)\mathrm{to}\left(g\right):I<I_{\mathrm{0}}\\ \left(g\right)\mathrm{to}\left(d\right):I+v>V_{\mathrm{0}}\end{array}\right.\end{array}$$

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 CO₂ 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.

$$\mathrm{ISI}\left(\mathit{\tau }\right)=\sum _i{\mathit{\beta }}_i(W_i\cdot \mathrm{86}\mathrm{400})$$

Here, W_i is the mean insolation of day i in W m⁻² and β_i equals one if W_i>τ 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.

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.

2.3 Optimal model parameters

The conceptual model relies on a small number of model parameters: τ_i, τ_g, τ_d and the two thresholds for the switch from a glaciation state to a deglaciation state and the inverse: V₀ and I₀. For all the simulations performed, we kept τ_i=9 kyr, τ_g=30 kyr, τ_d=12 kyr and I₀=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₀, on the model results. To compare our model results to the data, we used the benthic δ¹⁸O stack “LR04” (Lisiecki and Raymo, 2005) as a proxy for ice volume, considering that most of the δ¹⁸O changes of benthic foraminifera represent changes in continental ice (Shackleton, 1967; Shackleton and Opdyke, 1973). Lower δ¹⁸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 δ¹⁸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₁ 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₂ as $c_{\mathrm{2}}=\mathrm{1}-w$, with w the fraction of wrong transitions. The c₂ 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_{\mathrm{1}}\cdot c_{\mathrm{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₀ over the Quaternary, it was divided into five 500 kyr periods. The V₀ values that maximize the accuracy criteria for each time period and insolation forcing are called “optimal V₀”. To determine the optimal V₀ 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₀ thresholds (from V₀=1.0 to V₀=8.0 with a step of 0.1) and different initial conditions (initial volume V_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₀ changing with time that is equal to the corresponding optimal V₀ at each time period. Additionally, simulations were performed to determine the optimal V₀ threshold obtained when the optimization procedure is carried out over the whole Quaternary. It is called $V_{\mathrm{0}}^{\mathrm{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₀ 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 δ¹⁸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 δ¹⁸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 corresponding 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₀ and corresponding accuracy

For each insolation, we computed the deglaciation threshold V₀ maximizing the accuracy for each of the five 500 kyr periods and the fixed $V_{\mathrm{0}}^{\mathrm{Q}}$ value maximizing the accuracy over the whole Quaternary. The results are displayed in Fig. 2.

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

In some cases, several values of the deglaciation threshold V₀ lead to the same accuracy criteria, whereas in others there is only one V₀ value maximizing the accuracy. When several V₀ 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₀ thresholds. But the most striking feature is the increase of the optimal V₀ 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₀ over the whole Quaternary $V_{\mathrm{0}}^{\mathrm{Q}}$ is between 3.4 and 4 for each insolation type. It is a value in the 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₀ threshold for each time period as well as to the fixed $V_{\mathrm{0}}^{\mathrm{Q}}$ value maximizing the accuracy over the whole Quaternary is displayed in Fig. 3.

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_{\mathrm{0}}^{\mathrm{Q}}$).

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_{\mathrm{0}}^{\mathrm{Q}}$ value) is globally lower than the accuracy for each time period. This is due to the fact that the $V_{\mathrm{0}}^{\mathrm{Q}}$ values obtained are lower than the optimal V₀ values in the later part of the Quaternary and higher than the optimal V₀ 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 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 δ¹⁸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 δ¹⁸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 δ¹⁸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₀ 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₀ threshold often imply a gradual cooling, due to lower pCO₂ (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₂ trend, but rather only a decrease in minimal pCO₂ 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₀ 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.

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

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₀. 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₀ 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₀ 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₀ threshold over this period. Some differences are however noticeable between the different forcings. In particular the agreement for the ISI above 300 W m⁻² 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 W m⁻² forcing, the range of V₀ values that allow us to correctly reproduce most of the terminations during the last 1 million years is reduced (only values of V₀=3.9–4.0), whereas the results are more robust for the three other insolation forcings, with a broader range of working V₀ values. The ISI above 300 W m⁻² 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 $\mathrm{d}V/\mathrm{d}t=-I/{\mathit{\tau }}_{\mathrm{i}}$ 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 dependance 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 W m⁻² 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 W m⁻² is used.

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₀ is calibrated over the [100–1000] ka period and maintained for the whole simulation (V₀=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₀ over the last cycle (raised to V₀=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₀ 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.

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.

3.3 Reflections about the future

To model future natural evolutions of the climate system, possible evolutions of the V₀ 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₀ 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 glaciation, which is determined by the I₀ parameter. In our model, a broad range of I₀ parameter values lead to satisfying results over the last 1 million years. These different values projected 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 dynamics 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 simulation 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.