Articles | Volume 18, issue 3
Research article
25 Mar 2022
Research article |  | 25 Mar 2022

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

Gaëlle Leloup and Didier Paillard

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 Paillard2012). 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 Raymo2005) 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 δ18O records. It is known that the variations of annual total energy are of too-low an amplitude to explain such changes (Croll1864). 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 (Berger2021). 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 (Reeh1991), 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 Imbrie1980), 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 (Raymo1997; Lisiecki2010), obliquity (Huybers and Wunsch2005; Liu et al.2008), or both (Huybers2011; Parrenin and Paillard2012), to internal oscillations phase locked to the astronomical forcing (Saltzman et al.1984; Paillard1998; Gildor and Tziperman2000; Tziperman et al.2006), to internal oscillations independent of the astronomical forcing (Saltzman and Sutera1987; Toggweiler2008) 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 Paillard2003). 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 (Paillard2015). 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 Paillard2012). 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 V0, 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 (Paillard1998; Parrenin and Paillard2003, 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

(1) ( g ) d v d t = - I τ i + 1 τ g ( d ) d v d t = - I τ i - v τ d ,

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 (Raymo1997; Parrenin and Paillard2003, 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 Calov2011). 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.

(2) ( d ) to ( g ) : I < I 0 ( g ) to ( d ) : I + v > V 0

The idea that the deglaciation threshold is linked to both insolation and ice volume is not new (Parrenin and Paillard2003, 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 (Oerlemans1980; Pollard1982; 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 (MacAyeal1992; Marshall and Clark2002), ice sheet margin calving (Pollard1982), 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 Marshall1995; Ganopolski and Calov2011) 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 Parrenin2004; 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 (Paillard1998; Parrenin and Paillard2003). 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 Paillard2012). 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−2). 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−2. 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.


Here, Wi is the mean insolation of day i in W m−2 and βi equals one if Wi>τ 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−2 and ISI above 400 W m−2) 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: V0 and I0. For all the simulations performed, we kept τi=9 kyr, τg=30 kyr, τd=12 kyr and I0=0, as these parameters gave correct behaviour in previous studies (Parrenin and Paillard2003). 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, V0, on the model results. To compare our model results to the data, we used the benthic δ18O stack “LR04” (Lisiecki and Raymo2005) as a proxy for ice volume, considering that most of the δ18O changes of benthic foraminifera represent changes in continental ice (Shackleton1967; Shackleton and Opdyke1973). Lower δ18O 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 δ18O 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 c1 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 c2 as c2=1-w, with w the fraction of wrong transitions. The c2 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=c1c2. 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 V0 over the Quaternary, it was divided into five 500 kyr periods. The V0 values that maximize the accuracy criteria for each time period and insolation forcing are called “optimal V0”. To determine the optimal V0 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 V0 thresholds (from V0=1.0 to V0=8.0 with a step of 0.1) and different initial conditions (initial volume Vinit 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 V0 changing with time that is equal to the corresponding optimal V0 at each time period. Additionally, simulations were performed to determine the optimal V0 threshold obtained when the optimization procedure is carried out over the whole Quaternary. It is called V0Q in the following.

2.4 Definition of the deglaciation state in the data

To calculate our accuracy criteria c and therefore determine the optimal V0 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 δ18O 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 δ18O (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 V0 and corresponding accuracy

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

Figure 2Optimal deglaciation threshold V0 over the five 500 kyr periods for the four different summer insolation forcings at 65 N and optimal constant V0 threshold obtained when the optimization procedure is done over the whole Quaternary (V0Q). When several values of the deglaciation threshold V0 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 V0 lead to the same accuracy criteria, whereas in others there is only one V0 value maximizing the accuracy. When several V0 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 V0 thresholds. But the most striking feature is the increase of the optimal V0 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 V0 over the whole Quaternary V0Q 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 V0 threshold for each time period as well as to the fixed V0Q value maximizing the accuracy over the whole Quaternary is displayed in Fig. 3.

Figure 3Accuracy 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 V0Q).


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 V0Q value) is globally lower than the accuracy for each time period. This is due to the fact that the V0Q values obtained are lower than the optimal V0 values in the later part of the Quaternary and higher than the optimal V0 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 Raymo2005), 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 δ18O 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 Lisiecki2016), 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 δ18O as an ice volume proxy has already been largely debated (Shackleton1967; Chappell and Shackleton1986; Shackleton and Opdyke1973; 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 Raymo2005). We refer the reader to Raymo et al. (2018) for a review of possible biases in the interpretation of the LR04 benthic δ18O 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 V0 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 V0 threshold often imply a gradual cooling, due to lower pCO2 (Raymo1997; Paillard1998; Berger1978), but this idea is not supported by more recent data, which do not show a gradual decrease in the pCO2 trend, but rather only a decrease in minimal pCO2 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 Pollard1998; Clark et al.2006), or a sea ice switching mechanism (Gildor and Tziperman2000).

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 Gildor2003; Raymo and Nisancioglu2003). In contrast, others (Ashkenazy and Tziperman2004) 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 V0 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 4Best 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−2 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 V0. 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 V0 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 V0 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 V0 threshold over this period. Some differences are however noticeable between the different forcings. In particular the agreement for the ISI above 300 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 W m−2 forcing, the range of V0 values that allow us to correctly reproduce most of the terminations during the last 1 million years is reduced (only values of V0=3.9–4.0), whereas the results are more robust for the three other insolation forcings, with a broader range of working V0 values. The ISI above 300 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/τ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−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 W m−2 is used.

On the oldest part of the Quaternary, the caloric season forcing and the ISI above 300 W m−2 seem to perform better in reproducing the frequency of the oscillations than the summer solstice insolation and ISI above 400 W m−2. 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−2 and c=0.87 for the ISI above 300 W m−2). This is mainly due to the fact that for the three other forcings (caloric season, ISI above 300 and ISI above 400 W m−2), 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−2 and 0.89 for ISI above 300 W m−2). Figure 5 displays two different alternatives. For the first one (full line), the optimal V0 is calibrated over the [100–1000] ka period and maintained for the whole simulation (V0=5.0 for the summer solstice insolation, 4.5 for the caloric season, 4.0 for the ISI above 300 W m−2 and 4.6 for the ISI above 400 W m−2). 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 V0 over the last cycle (raised to V0=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 V0 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 5Normalized 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−2 and ISI above 400 W m−2. 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 V0 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 V0 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 V0 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 I0 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 I0 parameter. In our model, a broad range of I0 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 Loutre2002; Paillard2001), 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 CO2, which would affect the glaciation and deglaciation thresholds (Archer and Ganopolski2005; Talento and Ganopolski2021). 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 V0. 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 CO2 (Archer and Ganopolski2005; Talento and Ganopolski2021). 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 (Paillard2017). 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.

Code and data availability

The model code, insolation input files, spectral analysis and code needed to reproduce the figures are available for download: (Leloup2022).


The supplement related to this article is available online at:

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.


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.


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.


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,, 2013. a, b

Archer, D. and Ganopolski, A.: A movable trigger: Fossil fuel CO2 and the onset of the next glaciation, Geochem. Geophy. Geosy., 6, Q05003,, 2005. a, b

Ashkenazy, Y. and Tziperman, E.: Are the 41 kyr glacial oscillations a linear response to Milankovitch forcing?, Quaternary Sci. Rev., 23, 1879–1890,, 2004. a

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

Berger, A.: Milankovitch, the father of paleoclimate modeling, Clim. Past, 17, 1727–1733,, 2021. a

Berger, A. and Loutre, M.-F.: Insolation values for the climate of the last 10 million years, Quaternary Sci. Rev., 10, 297–317,, 1991. a

Berger, A. and Loutre, M.-F.: Climate: An Exceptionally Long Interglacial Ahead?, Science, 297, 1287–1288,, 2002. a

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

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

Calder, N.: Arithmetic of ice ages, Nature, 252, 216–218,, 1974. a

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

Clark, P. and Pollard, D.: Origin of the Middle Pleistocene Transition by ice sheet erosion of regolith, Paleoceanography, 13, 1–9,, 1998. a

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 pCO2, Quaternary Sci. Rev., 25, 3150–3184,, 2006. a, b

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,, 1864. a, b

Elderfield, H., Ferretti, P., Greaves, M., Crowhurst, S., Mccave, I. N., Hodell, D., and Piotrowski, A. M.: Evolution of Ocean Temperature and Ice Volume Through the Mid-Pleistocene Climate Transition, Science, 337, 704–709,, 2012. a

Gallée, H., Van Yperselb, J. P., Fichefet, T., Marsiat, I., Tricot, C., and Berger, A.: Simulation of the last glacial cycle by a coupled, sectorially averaged climate-ice sheet model: 2. Response to insolation and CO2 variations, J. Geophys. Res.-Atmos., 97, 15713–15740,, 1992. a

Ganopolski, A. and Calov, R.: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles, Clim. Past, 7, 1415–1425,, 2011. a, b

Gildor, H. and Tziperman, E.: Sea ice as the glacial cycles’ Climate switch: role of seasonal and orbital forcing, Paleoceanography, 15, 605–615,, 2000. a, b

Hays, J., Imbrie, J., and Shackleton, N. J.: Variations in the Earth's Orbit: Pacemaker of the Ice Ages, Science, 194, 1121–1132,, 1976. a, b

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

Huybers, P.: Early Pleistocene glacial cycles and the integrated summer insolation forcing, Science, 313, 508–511,, 2006. a, b, c, d, e

Huybers, P.: Combined obliquity and precession pacing of late Pleistocene deglaciations, Nature, 480, 229–232,, 2011. a

Huybers, P. and Wunsch, C.: Obliquity pacing of the late Pleistocene glacial terminations, Nature, 434, 491–494,, 2005. a

Imbrie, J. and Imbrie, J. Z.: Modeling the Climatic Response to Orbital Variations, Science, 207, 943–953,, 1980. a, b, c

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

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

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J., Chappellaz, J., Fischer, H., Gallet, J.-C., Johnsen, S., Leuenberger, M., Loulergue, 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,, 2007. a

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

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

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

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],, 2022 a

Lisiecki, L. E.: Links between eccentricity forcing and the 100,000-year glacial cycle, Nat. Geosci., 3, 349–352,, 2010. a

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records, Paleoceanography, 20, PA1003, 2005. a, b, c, d

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

MacAyeal, D. R.: Irregular oscillations of the West Antarctic ice sheet, Nature, 359, 29–32,, 1992. a

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

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

Oerlemans, J.: Model experiments on the 100,000-yr glacial cycle, Nature, 287, 430–432,, 1980. a

Paillard, D.: The timing of Pleistocene glaciations from a simple multiple-state climate model, Nature, 391, 378–381,, 1998. a, b, c, d, e, f, g, h

Paillard, D.: Glacial cycles: toward a new paradigm, Rev. Geophys., 39, 325–346,, 2001. a

Paillard, D.: Quaternary glaciations: from observations to theories, Quaternary Sci. Rev., 107, 11–24,, 2015. a

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

Paillard, D. and Parrenin, F.: The Antarctic ice sheet and the triggering of deglaciations, Earth Planet. Sci. Lett, 227, 263–271,, 2004. a

Paillard, D., Labeyrie, L., and Yiou, P.: Macintosh Program performs time-series analysis, Eos Transactions, AGU, 77, 379–379,, 1996. a

Parrenin, F. and Paillard, D.: Amplitude and phase of glacial cycles from a conceptual model, Earth Planet. Sci. Lett., 214, 243–250,, 2003. a, b, c, d, e, f, g, h, i, j

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,, 2012. a, b, c, d, e, f, g

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

Pollard, D.: A simple ice sheet model yields realistic 100 kyr glacial cycles, Nature, 296, 334–338,, 1982. a, b

Raymo, M.: The timing of major climate terminations, Paleoceanography, 12, 577–585,, 1997. a, b, c

Raymo, M. E. and Nisancioglu, K. H.: The 41 kyr world: Milankovitch's other unsolved mystery, Paleoceanography, 18, 1011,, 2003. a

Raymo, M. E., Kozdon, R., David, E., Lisiecki, L., and Ford, H. L.: The accuracy of mid-Pliocene δ18O-based ice volume and sea level reconstructions, Earth-Sci. Rev., 177, 291–302,, 2018. a

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

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

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,<0236:TMQCTA>2.0.CO;2, 1987. a

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,<3380:TLQGAT>2.0.CO;2, 1984. a

Shackleton, N.: Oxygen Isotope Analyses and Pleistocene Temperatures Re-assessed, Nature, 215, 15–17,, 1967. a, b

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,, 1973. a, b

Spratt, R. M. and Lisiecki, L. E.: A Late Pleistocene sea level stack, Clim. Past, 12, 1079–1092,, 2016.  a

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

Toggweiler, J. R.: Origin of the 100,000-year timescale in Antarctic temperatures and atmospheric CO2, Paleoceanography, 23, PA2211,, 2008. a

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,, 2017. a, b, c, d, e, f, g, h, i

Tziperman, E. and Gildor, H.: On the mid-Pleistocene transition to 100-kyr glacial cycles and the asymmetry between glaciation and deglaciation times, Paleoceanography, 18, 1-1–1-8,, 2003. a

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, Paleoceanography, 21, PA4206,, 2006. a

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

Verbitsky, M. Y., Crucifix, M., and Volobuev, D. M.: A theory of Pleistocene glacial rhythmicity, Earth Syst. Dynam., 9, 1025–1043,, 2018. a

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

Short summary
Over the last 2.6 Myr, the Quaternary period has been marked by the alternation of extended and reduced Northern Hemisphere ice sheets, known as glacial-interglacial cycles. With a simple model, we are able to reproduce the main features of the ice volume evolution, like the switch of periodicity, from 41 kyr cycles to 100 kyr cycles, observed in the data after 1 Ma. The quality of the model-data agreement depending on the input insolation and period considered is discussed.