Influence of the choice of insolation forcing on the results of a conceptual glacial cycle model Inﬂuence 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 inﬂuence 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 inﬂuence of the use of different summer insolation forcings, having different contributions from precession and obliquity, on the model results. We show that some


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 glacialinterglacial 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 δ 18 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 0 , and keeping all the other model parameters constant. This is similar to the results of Paillard (1998), who obtained a frequency shift on the glacial cycles by using a linearly increasing deglaciation threshold. This is also coherent with the more recent results of Tzedakis et al. (2017), who demonstrated that the particular sequence of interglacials that happened over the Quaternary and the frequency shift from 41 to 100 kyr could be explained with a simple rule, taking into account a deglaciation threshold that increases over time, leading to "skipped" insolation peaks and longer cycles.

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;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 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;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 0 : 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.
The idea that the deglaciation threshold is linked to both insolation and ice volume is not new 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 2 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).

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 −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, W i is the mean insolation of day i in W m −2 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.

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 0 and I 0 . For all the simulations performed, we kept τ i = 9 kyr, τ g = 30 kyr, τ d = 12 kyr and I 0 = 0, as these parameters gave correct behaviour in previous studies (Parrenin and Paillard, 2003). No attempt was made to tune these parameters to obtain a behaviour closer to the data: on the contrary, focus was set more towards the influence of the deglaciation threshold parameter, V 0 , on the model results. To compare our model results to the data, we used the benthic δ 18 O stack "LR04" (Lisiecki and Raymo, 2005) as a proxy for ice volume, considering that most of the δ 18 O changes of benthic foraminifera represent changes in continental ice (Shackleton, 1967;Shackleton and Opdyke, 1973). Lower δ 18 O values correspond to lower ice volume. The model results and the LR04 curve were normalized to facilitate their comparison. In the following, "data" refers to the normalized δ 18 O stack curve LR04.
To estimate which model parameters lead to model results closer to the data, the definition of an objective criteria is needed. The choice of such a criteria is not straightforward, and the use of different criteria could have led to slightly different results. Our model is simple and does not aim at precisely reproducing the ice volume evolution, but rather at reproducing the main qualitative features, such as the shape and frequency of the oscillations. Therefore, we used a criteria based on the state of the system: glaciation or deglaciation. Similar results can be obtained using a simple correlation coefficient (see Sect. S2). The definition of the deglaciation state in the data is explained in Sect. 2.4. A critical point is that our model should be able to deglaciate at the right time, when a deglaciation is seen in the data. Conversely, the model should not produce deglaciation at periods where the data do not show deglaciation.
We defined a criteria for each of these two conditions, and assembled them in a global criteria. To determine if deglaciation is well reproduced by our model, we look at the state of the model (glaciation or deglaciation) at the time halfway between the start of the deglaciation and the end of the deglaciation. If the model state at that time is deglaciation, the deglaciation is considered as correctly reproduced. Otherwise, it is considered as a "missed" deglaciation. We simply defined the criteria c 1 as the fraction of deglaciations correctly reproduced. This criteria equals one when all the deglaciations take place at the right time, meaning that the model produces a deglaciation state every time deglaciation is seen in the data. To ensure that the model does not deglaciate too often, we looked at insolation maxima that are not associated with deglaciation in the data and ensured the corresponding model state was glaciation. We decided to look specifically at insolation maxima, as they are where the model is the most likely to deglaciate when it should not. A deglaciation seen in the model at a place where no deglaciation is seen in the data is considered as a "wrong transition". We defined the criteria c 2 as c 2 = 1 − w, with w the fraction of wrong transitions. The c 2 criteria equals one when the model does not deglaciate when it should not. To take into account these two conditions, the overall accuracy criteria c is defined as c = c 1 · c 2 . It is equal to one when all the deglaciations are correctly placed and when no additional deglaciation compared to the data take place.
In order to study the evolution of the optimal deglaciation threshold V 0 over the Quaternary, it was divided into five 500 kyr periods. The V 0 values that maximize the accuracy criteria for each time period and insolation forcing are called "optimal V 0 ". To determine the optimal V 0 threshold corresponding to each period and insolation forcing, several simulations were carried out and the parameter values maximizing the accuracy criteria c were selected. More precisely, for each insolation and period, 3500 simulations corresponding to different V 0 thresholds (from V 0 = 1.0 to V 0 = 8.0 with a step of 0.1) and different initial conditions (initial volume V init ranging from 0.0 to 5.0 with a step of 0.2, and the two possible initial states -glaciation or deglaciation) were performed. The numerical integration of the model equations was done with a fourth-order Runge-Kutta scheme.
For each insolation forcing, the best fit over the Quaternary is defined as the simulation over the whole Quaternary (0 to 2500 ka) with V 0 changing with time that is equal to the corresponding optimal V 0 at each time period. Addition- ally, simulations were performed to determine the optimal V 0 threshold obtained when the optimization procedure is carried out over the whole Quaternary. It is called V Q 0 in the following.

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

Optimal deglaciation threshold V 0 and corresponding accuracy
For each insolation, we computed the deglaciation threshold V 0 maximizing the accuracy for each of the five 500 kyr periods and the fixed V Q 0 value maximizing the accuracy over the whole Quaternary. The results are displayed in Fig. 2.
In some cases, several values of the deglaciation threshold V 0 lead to the same accuracy criteria, whereas in others there is only one V 0 value maximizing the accuracy. When several V 0 values are equivalent, the mean value was plotted and the other possible values are represented with error bars.
Over the same time period, different insolation values lead to slightly different optimal V 0 thresholds. But the most striking feature is the increase of the optimal V 0 threshold over time, and more specifically for the last 1 million years, after the MPT. This feature is valid regardless of the insolation type. The optimal V 0 over the whole Quaternary V Q 0 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 0 threshold for each time period as well as to the fixed V Q 0 value maximizing the accuracy over the whole Quaternary is displayed in Fig. 3.
It is first noticeable that the accuracy is higher over the last 1 million-year period, regardless the input summer insolation used. In the last 1 million years, the summer solstice insolation as input produces the best results. However, this is no longer the case for older time periods: the solstice insolation gives the worst results at the start of the Quaternary. The accuracy obtained for the whole Quaternary period (fixed V Q 0 value) is globally lower than the accuracy for each time period. This is due to the fact that the V Q 0 values obtained are lower than the optimal V 0 values in the later part of the Quaternary and higher than the optimal V 0 values in the earliest part of the Quaternary, leading to a poorer representation of both of these periods.
In the earlier part of the Quaternary (periods earlier than 1.5 Ma), the results are less robust. This is due to increased uncertainties in the LR04 record, and the associated definition of interglacial periods, which affects our accuracy criteria. In their classification, Tzedakis et al. (2017) stated that the definition of the valley depth needed to separate several interglacials was quite straightforward in the earlier part of 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 δ 18 O LR04 curve includes at the same time an ice volume and deep water temperature component. Ice volume and sea level reconstructions do exist (Bintanja et al., 2005;Spratt and Lisiecki, 2016), but are however limited to the more recent part of the Quaternary and do not allow for the investigation of the pre-MPT period. The use of δ 18 O as an ice volume proxy has already been largely debated (Shackleton, 1967;Chappell and Shackleton, 1986;Shackleton and Opdyke, 1973;Clark et al., 2006), and recent studies (Elderfield et al., 2012) have shown that the temperature component may be as large as 50 %. Furthermore, the stack was tuned to insolation (Lisiecki and Raymo, 2005). We refer the reader to Raymo et al. (2018) for a review of possible biases in the interpretation of the LR04 benthic δ 18 O stack as an ice volume and sea level reconstruction. All these reasons encourage us to remain at a qualitative level to fit the data.
The increase in the optimal deglaciation threshold V 0 over the MPT is however a robust feature that does not depend on the input summer insolation forcing used. This is coherent with the work of Paillard (1998) and Tzedakis et al. (2017), who obtained a frequency shift from 41 to 100 kyr cycles with an increasing deglaciation threshold. The classical hypothesis and possible physical meanings of the rise of the V 0 threshold often imply a gradual cooling, due to lower pCO 2 (Raymo, 1997;Paillard, 1998;Berger, 1978), but this idea is not supported by more recent data, which do not show a gradual decrease in the pCO 2 trend, but rather only a decrease in minimal pCO 2 values (Hönisch et al., 2009;Yan et al., 2019). Other hypothesis imply a gradual erosion of a thick regolith layer that exposed the ice sheet to a higher friction bedrock, allowing for a thicker ice sheet to develop (Clark and Pollard, 1998;Clark et al., 2006), or a sea ice switching mechanism (Gildor and Tziperman, 2000).
The overall lowest accuracy in the older part of the record suggests that our non-linear threshold model is less adapted for this time period. Indeed, the ice volume might respond more linearly to the insolation forcing before the MPT, as some studies suggest (Tziperman and Gildor, 2003;Raymo and Nisancioglu, 2003). In contrast, others (Ashkenazy and Tziperman, 2004) have shown that the 41 kyr pre-MPT oscillations are in fact significantly asymmetric and therefore suggested that oscillations both before and after the MPT could be explained as the self-sustained variability of the climate system, phase locked to the astronomical forcing. In our model, the deglaciation threshold V 0 changes over time, leading to an amplitude change of the cycles. However, it cannot be excluded that the mechanisms behind the pre-and post-MPT glacial cycles are structurally different, and that they cannot be explained with the same physics and the same equations in our model. Moreover, this kind of conceptual model has mainly been constructed in order to explain the non-linear features of the 100 kyr cycles, and it is therefore not surprising that its agreement with the data on the pre-MPT period is less satisfying.

Best fit over the Quaternary
Our conceptual model is able to reproduce qualitatively well the data (LR04 normalized curve) over the whole Quaternary. The model's best fit over the Quaternary for each insolation forcing, as defined in Sect. 2.3, is displayed in Fig. 4. It is able to reproduce the frequency shift from a dominant 41 kyr period before 1 Ma to longer cycles afterwards, as observed in the data, and thus by varying only one parameter during the whole simulation length: the deglaciation threshold V 0 . Like the previous model of Parrenin and Paillard (2003), from which our model is adapted, ours quickly looses its sensitivity to the initial conditions (after no more than 200 kyr). For every input forcing, longer cycles are obtained in the last part of the Quaternary (last Myr). Figure 4 displays the results over the whole Quaternary with the V 0 threshold being set to its optimal value for each 500 kyr period, while Fig. 5 displays the results over the last 1 million years with the V 0 threshold being set to its optimal value for the [0-1000] ka period.
For the last 1 million years, it is possible to reproduce all terminations with the right timing, apart from the last deglaciation, for all insolation forcings, by using a single value of the V 0 threshold over this period. Some differences are however noticeable between the different forcings. In particular the agreement for the ISI above 300 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 V 0 values that allow us to correctly reproduce most of the terminations during the last 1 million years is reduced (only values of V 0 = 3.9-4.0), whereas the results are more robust for the three other insolation forcings, with a broader range of working V 0 values. The ISI above 300 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 V 0 is calibrated over the [100-1000] ka period and maintained for the whole simulation (V 0 = 5.0 for the summer solstice insolation, 4.5 for the caloric season, 4.0 for the ISI above 300 W m −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 V 0 over the last cycle (raised to V 0 = 5.5 from 100 ka onward). In this case, the model does accurately reproduce the last deglaciation for all insolation forcings. This suggests that we might need to further raise the V 0 threshold in order to model a theoretical "natural" future (without anthropogenic influence) with longer cycles. A cycle is not enough to conclude about a trend, but this should be envisaged for future natural scenarios. This is coherent with the idea of Paillard (1998), who used a linearly increasing deglaciation threshold over the Quaternary.

Reflections about the future
To model future natural evolutions of the climate system, possible evolutions of the V 0 threshold should be considered. However, we do not exclude the fact that variations of other parameters, which were kept constant in this study, could vary in the future. For instance, different I 0 thresholds have to be considered. The fate of the next thousands of years in a natural scenario cannot be ascertained with our conceptual model. Indeed, this depends on the onset (or not) of a glaciation, which is determined by the I 0 parameter. In our model, a broad range of I 0 parameter values lead to satisfying results over the last 1 million years. These different values pro-Figure 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 −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 V 0 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 0 threshold). This results in a green shading when deglaciations are seen in the model and data at the same time.
jected into the future lead either to, in one case, a glaciation start at present, or, in the other case, a continued deglaciation, switching to a glaciation state only around 50 kyr after present at the next insolation minima. This is due to the current particular astronomical configuration, with very low eccentricity, which leads to high summer insolation minima, where the threshold for glaciation might not be reached. This is coherent with previous studies (Berger and Loutre, 2002;Paillard, 2001), which suggested that without anthropogenic forcing, the present interglacial might have lasted 50 kyr.
However, this exercise is purely academic, as we are not taking into account the role of anthropogenic CO 2 , 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.

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 0 . 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 2 (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.
Code and data availability. The model code, insolation input files, spectral analysis and code needed to reproduce the figures are available for download: https://doi.org/10.5281/zenodo.6045532 (Leloup, 2022).
Author contributions. GL and DP designed the study. GL performed the simulations, and also wrote the manuscript under the supervision of DP.
Competing interests. The contact author has declared that neither they nor their co-author has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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