Insolation evolution and ice volume legacies determine interglacial and glacial intensity

Interglacials and glacials represent low and high ice volume end-members of ice age cycles. While progress has been made in our understanding of how and when transitions between these states occur, their relative intensity has been lacking an explanatory framework. With a simple quantitative model, we show that over the last 800 000 years interglacial intensity can be described as a function of the strength of the previous glacial and the summer insolation at high latitudes in both hemispheres during the deglaciation. Since the precession components in the boreal and austral insolations counteract each other, the amplitude increase in obliquity cycles after 430 000 years ago is imprinted in interglacial intensities, contributing to the manifestation of the so-called Mid-Brunhes Event. Glacial intensity is also linked to the strength of the previous interglacial, the time elapsed from it, and the evolution of boreal summer insolation. Our results suggest that the memory of previous climate states and the time course of the insolation are crucial for understanding interglacial and glacial intensities.

1 Introduction

The most prominent climate signal in the last 800 000 years (800 kyr) is the alternating pattern of glacials (with large ice sheets at high northern latitudes and an extended Antarctic Ice Sheet – AIS) and interglacials (with little Northern Hemisphere – NH – ice outside Greenland, Past Interglacials Working Group of PAGES, 2016) (Fig. 1d–g). It is generally accepted that this pattern is initiated and paced by changes in insolation caused by astronomical changes in orbit and axial tilt (Milankovitch, 1941) (Fig. 1a–c), modulated by strong feedbacks, including those of the carbon cycle (Fig. 1d) and ice sheet albedo. Both conceptual and sophisticated models have successfully reproduced many of the features of the observed record (e.g. Huybers, 2011; Parrenin and Paillard, 2012; Verbitsky et al., 2018; Berger et al., 1999; Abe-Ouchi et al., 2013; Willeit et al., 2019). However, it remains challenging to state simply which features of astronomical forcing determine the timing of glacial terminations or the amplitude of glacial cycles. As a result, the holy grail of Quaternary palaeoclimate – to take only the external Milankovitch forcing and predict accurately the sequence and strength of glacial cycles – remains elusive. Recently a simple rule (Tzedakis et al., 2017) was successful in predicting the timing of the occurrence of interglacials but did not address their intensities.

Figure 1 Changes in climate conditions over the last 800 kyr. (a) Eccentricity (red) and precession parameter (black) (Laskar et al., 2004). (b) Obliquity (Laskar et al., 2004). (c) Caloric summer half-year insolation at 65^∘ N (orange) and at 65^∘ S (grey) based on the orbital solution of Laskar et al. (2004). (d) Compilation of atmospheric CO₂ records from Antarctic ice cores (Bereiter et al., 2015, and references therein; Nehrbass-Ahles et al., 2020; Shin et al., 2020; Bauska et al., 2021). (e) EPICA Dome C sea salt Na flux, a proxy for sea ice extent (Wolff et al., 2006). (f) Global average surface temperature as a temperature deviation from the present (average over 0–5 kyr BP) (Snyder, 2016). (g) Stack of benthic δ¹⁸O records, LR04 (Lisiecki and Raymo, 2005). Ice core data are on the AICC2012 timescale (Bazin et al., 2013).

Here, interglacial or glacial intensity (strength) refers to the magnitude of a measurable property that can be extracted from proxy climate records. Various metrics integrating global climate effects have been employed, including surface temperature (Snyder, 2016; Past Interglacials Working Group, 2016), the oxygen isotope ratio (δ¹⁸O) in benthic Foraminifera in a stack of globally distributed records (Lisiecki and Raymo, 2005), sea-level reconstructions based on corals (Dutton et al., 2015), hydraulic control models of semi-isolated basins (Grant et al., 2014), and δ¹⁸O of seawater (Elderfield et al., 2012). In terms of interglacial strength, the metrics converge in their broad trends over the last 800 kyr (Fig. 1): before the Mid-Brunhes Event (MBE; ∼430 kyr before present – BP), interglacials (Marine Isotope Stages – MISs – 19c, 17, 15e, 15a, and 13a) were weaker (cooler, higher benthic δ¹⁸O, atmospheric CO₂ lower than pre-industrial concentrations) (Berger and Wefer, 2003; Tzedakis et al., 2009). The strongest interglacials (MIS11c, 9e, 5e, and 1) occurred after the MBE, although MIS7e and MIS7c are closer in intensity to pre-MBE interglacials. With respect to temperature and sea level, MIS5e and MIS11c are the most prominent interglacials, followed by MIS9e and MIS1 (Elderfield et al., 2012; Grant et al., 2014; Dutton et al., 2015; Snyder, 2016; Past Interglacials Working Group of PAGES, 2016). The change in interglacial intensity at the MBE has been attributed to different factors, including an increase in the amplitude of obliquity cycles (Fig. 1b) (EPICA Community Members, 2004; Mitsui and Boers, 2022), increased instability of the West Antarctic Ice Sheet (Holden et al., 2011), changes in southern westerlies and Southern Ocean ventilation (Yin, 2013), and a reduction in the volume of interglacial Antarctic Bottom Water, contributing to a greater release of CO₂ from the deep ocean (Barth et al., 2018). Differences in the intensities of individual interglacials have been discussed in terms of the contribution of insolation and greenhouse gases (Yin and Berger, 2012; Obase et al., 2021) and also in relation to the ice sheet size of the preceding glacial (Raymo, 1997; Paillard, 1998; Berger and Wefer, 2003; Lang and Wolff, 2011). However, a simple explanatory framework providing a systematic understanding of these differences and trends has been lacking.

With respect to glacial strength, sea-level reconstructions (Elderfield et al., 2012; Grant et al., 2014; Shakun et al., 2015; Spratt and Lisiecki, 2016) and glacial geologic evidence (Hughes and Gibbard, 2018; Batchelor et al., 2019; Hughes et al., 2020) indicate that the largest increases in global ice volume occurred during MIS12, 16, 2, and 6, the smallest during MIS15b, 7d, 14, and 8, with MIS10, 18, and 20 somewhere in between. A decrease in boreal summer insolation is the primary trigger for glacial inception (Ganopolski et al., 2016), leading to rapid glacial advances in continental interiors and mountains (Hughes and Gibbard, 2018; Hughes et al., 2020). Within a glacial period, “excess” ice volume accumulates during intervals of low eccentricity–precession forcing, leading to weak summer insolation maxima and minimum ice ablation (Raymo, 1997; Paillard, 1998). These observations underline the importance of the evolution of boreal insolation and time for ice sheet buildup (Parrenin and Paillard, 2003) for differences in glacial strength.

In this work, we focus on the interglacial and glacial intensities over the last 800 kyr, where there is a broad consensus on which δ¹⁸O peaks correspond to interglacials (Past Interglacials Working Group of PAGES, 2016). The extension of our models beyond 800 kyr BP will be the subject of a future study. As a stepping stone to the more difficult target of predicting the amplitude of a long succession of glacial cycles from the forcing, here we address the question of predictability from one glacial trough to the next interglacial peak and vice versa over the last 800 kyr. If we know the astronomical forcing and some aspect of the preceding climate history, can we predict the strength of each interglacial and of each glacial? In anthropomorphic terms, could a bystander at a glacial maximum predict the strength of the interglacial that occurs about 20 kyr later, and could a bystander at the end of an interglacial predict the strength of the glacial maximum that might occur as long as 100 kyr later? We consider 11 interglacials over the last 800 kyr following the Past Interglacials Working Group of PAGES (2016) and also take into account the interglacial MIS21 for predicting the glacial intensity in MIS20. We use the LR04 benthic δ¹⁸O stack as a reference record (Lisiecki and Raymo, 2005) on the premise that the convolved ice volume and deep-water temperature signal is a robust integrated metric of intensities (Past Interglacials Working Group of PAGES, 2016). Thus, the whole purpose of the present paper is to predict the amplitude of δ¹⁸O from the insolation, the only external driver of the climate system. Atmospheric CO₂ or other climate feedbacks are considered agents in between insolation and δ¹⁸O changes, and their potential role is discussed in the final section.

2 Data and methods

Given that uncertainties of the order of ±20 m in sea-level reconstructions from deconvolved records of deep-water temperatures and δ¹⁸O of seawater (e.g. Elderfield et al., 2012) exceed differences in interglacial highstands (Dutton et al., 2015), we use the LR04 benthic δ¹⁸O stack record (Lisiecki and Raymo, 2005) (mean standard error of 0.06 ‰) as an integrated metric of interglacial intensities. Although benthic δ¹⁸O records contain local deep-water temperature and hydrographic effects, averaging several globally distributed records removes some of the regional variability. Interglacial intensities are measured by local minima in the LR04 stack at 2, 123, 217, 239, 329, 405, 491, 575, 610, 696, 780, and 858 kyr BP respectively (Fig. 2a as well as Fig. 4a). Glacial intensities are measured by local maxima in the same record at 18, 140, 223, 252, 341, 433, 536, 585, 630, 718, and 794 kyr BP (Fig. 2a as well as Fig. 4a). According to Elderfield et al. (2012), during glacials deep-water temperatures rapidly decline and stabilize as they reach the ocean-water freezing temperature, resembling a square wave function; the saw-tooth character of the benthic δ¹⁸O records, therefore, reflects the slow ice sheet buildup and rapid decay. By extension, differences in δ¹⁸O_max between glacial maxima in the LR04 record largely reflect differences in ice volume. We assume that the orbital tuning in the LR04 δ¹⁸O stack record is essentially correct, at least on orbital timescales. Thus, we take it for granted which insolation peak induces which interglacial. Under this assumption, we explore the relationships between the amplitude (not the timing) of δ¹⁸O peaks and the insolation forcing (see Tzedakis et al., 2017, for the timing, which explains how one or two obliquity cycles are skipped without having terminations).

Figure 2 Modelling interglacial intensities. (a) LR04 δ¹⁸O. The red circles indicate the minima of δ¹⁸O (δ¹⁸O_min) at each interglacial and the blue triangles the maxima (δ¹⁸O_max) at glacials. See below for the grey strips and the dashed lines. (b) Caloric summer half-year insolation at 65^∘ N (F_N, black) and 65^∘ S (F_S, green). The average of the two (magenta) is also shown. The blue dashed lines show timings t_s at which the caloric summer half-year insolation at 65^∘ N exceeds the average 5845 GJ m⁻² (black horizontal line) and the red dashed lines show timings t_e at which the insolation falls back below the average. Each termination starts roughly around t_s, and it is completed around t_e. Exceptionally, Termination III starts after the local insolation minimum at 254 kyr BP (orange dotted line), responding to the second rise in the insolation. The grey strips show the termination intervals [t_s, t_e] based on the insolation curve. (c) Integral of the caloric summer insolation anomaly between t_s and t_e at 65^∘ N, I_N (black cross), the integral at 65^∘ S for the same period, I_S (green diamond), and the average $I_{\text{AV}}=\frac{\mathrm{1}}{\mathrm{2}}(I_{\mathrm{N}}+I_{\mathrm{S}})$ (magenta square). (d) δ¹⁸O_max. (e) Predictions by linear regression models with explanatory variables in panels (c) and (d): Model 1 with I_N (black cross); Model 2 with both I_N and I_S with their own coefficients (blue diamond with cross); Model 3 with I_AV (magenta squares).

Figure 3 Scatter plot in the space: (a) (δ¹⁸O_max, I_N, δ¹⁸O_min), (b) (δ¹⁸O_max, I_AV, δ¹⁸O_min), and (c) (δ¹⁸O_max, β₂I_N+β₃I_S, δ¹⁸O_min). (a) We can separate the four strong interglacials (MIS11c, 9e, 5e, and 1) from the weak interglacials in the plane of explanatory variables (δ¹⁸O_max, I_N). In this case, however, the weak interglacial MIS15e is close to the strong interglacials in the ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}-I_{\mathrm{N}}$ plane. Replacing I_N with a mixture of I_N and I_S (b, c) leads to a clear separation of the four strongest interglacials (MIS11c, 9e, 5e, and 1) from the rest.

The caloric summer half-year insolation represents the amount of insolation integrated over the caloric summer half of the year, defined such that any day of the summer half receives more insolation than any day of the winter half (Milankovitch, 1941; Berger, 1978). Near 65^∘ N, the variance of this measure has almost equal contributions from climatic precession and the obliquity. Since ice ablation depends both on the strength of irradiation and the time period over which it is high (length of the summer season) (Milankovitch, 1941), we use this insolation metric in our models of interglacial and glacial intensities (Sect. 3). The caloric summer half-year insolation at 65^∘ N and at 65^∘ S is calculated at every 1 kyr by using R package palinsol (Crucifix, 2016) based on the orbital solution of Laskar et al. (2004). Both insolations have almost the same average $\stackrel{\mathrm{‾}}{F}=\mathrm{5.845}$ GJ m⁻² over the last 1×10⁶ years (Myr). The solar constant used is 1368 W m⁻².

In Sect. 3 we consider multiple linear regression models for interglacial intensities as well as for glacial intensities. We compare several candidate models using the Bayesian information criterion (BIC). The BIC is a common criterion for statistical model selection (Raftery, 1995) and is derived from the Bayes factor: the ratio of the marginal likelihoods of two competing models. Among different models that explain a dataset, the model with the lowest BIC is preferred. Any additional parameters in a model generally result in better or equally good fits to the data but can cause overfitting. Hence the BIC penalizes the increase in the number of parameters. When comparing two models (say model j against model k), the difference in BIC ($\mathrm{\Delta }\text{BIC}={\text{BIC}}_{\mathrm{k}}-{\text{BIC}}_{\mathrm{j}}$) is interpreted as the strength of the statistical evidence for model j against model k. As a rule of thumb (Raftery, 1995), the evidence is considered weak (or not worth more than a bare mention) if $\mathrm{0}<\mathrm{\Delta }\text{BIC}<\mathrm{2}$, positive if $\mathrm{2}<\mathrm{\Delta }\text{BIC}<\mathrm{6}$, strong if $\mathrm{6}<\mathrm{\Delta }\text{BIC}<\mathrm{10}$, and very strong if ΔBIC>10.

Table 1 Coefficients and statistics of the regression models for interglacial intensity δ¹⁸O_min (corresponding to Fig. 2). Model 1 (${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\mathrm{N}}$), Model 2 (${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\mathrm{N}}+{\mathit{\beta }}_{\mathrm{3}}{I}_{\mathrm{S}}$), and Model 3 (${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\text{AV}}$). The overall F test provides a p value of less than 0.01 in each model, which rejects the null hypothesis that none of the variables in the model is significant. The asterisks indicate the significance of each coefficient: ^∗ for $p\in (\mathrm{0.01},\mathrm{0.05}]$, ${}^{\ast \ast }$ for $p\in (\mathrm{0.001},\mathrm{0.01}]$, and ${}^{\ast \ast \ast }$ for $p\in [\mathrm{0},\mathrm{0.001}]$.

3 Results

3.1 Interglacial intensities

Figure 2a shows the LR04 benthic δ¹⁸O record over the last 800 kyr, where interglacial peaks (δ¹⁸O_min) are indicated by red circles and glacial peaks (δ¹⁸O_max) by blue triangles. In Fig. 2b, the caloric summer half-year insolation at 65^∘ N (F_N(t), black) and that at 65^∘ S (F_S(t), green) are shown respectively. The average of the two insolations ($\frac{\mathrm{1}}{\mathrm{2}}(F_{\mathrm{N}}+F_{\mathrm{S}}$), magenta) yields a variation correlated with obliquity (R=0.998) because the climatic precessions across the two hemispheres cancel each other (Milankovitch, 1941; Berger, 1978). Comparing Fig. 2a and b, we observe that each termination starts near the time when the boreal summer insolation F_N(t) exceeds its average $\stackrel{\mathrm{‾}}{F}=\mathrm{5.845}$ GJ m⁻². Exceptionally, the start of Termination III leading to MIS7e is delayed relative to the crossing point between F_N(t) and $\stackrel{\mathrm{‾}}{F}$ by about 11 kyr, possibly because the rise in the insolation F_N(t) stops (and slightly decreases) around a local minimum at 254 kyr BP. Termination III starts after the local minimum at 254 kyr BP, responding to the second rise in the insolation. We denote those upward crossing points between F_N(t) and $\stackrel{\mathrm{‾}}{F}$ associated with terminations by t_s (including an exception t_s=254 kyr BP for Termination III; the case without this exception is discussed below). We further observe that each termination is completed near the time (t_e) when the insolation F_N(t) falls below the average $\stackrel{\mathrm{‾}}{F}$. While most of the deglaciation intervals contain one insolation maximum, the interval leading to MIS17 contains two insolation maxima and the interval leading to MIS13a three insolation maxima.

We seek a simple regression relation describing interglacial intensity, represented by δ¹⁸O_min. We postulate that the magnitude of each termination is related to the boreal summer insolation anomaly integrated during the insolation-based termination period [t_s, t_e] at 65^∘ N and the austral summer insolation (at 65^∘ S) anomaly integrated during the same period:

$$I_{\mathrm{N}}=\sum _{t=t_{\mathrm{s}}}^{t_{\mathrm{e}}}\left(F_{\mathrm{N}}\right(t)-\stackrel{\mathrm{‾}}{F}),I_{\mathrm{S}}=\sum _{t=t_{\mathrm{s}}}^{t_{\mathrm{e}}}\left(F_{\mathrm{S}}\right(t)-\stackrel{\mathrm{‾}}{F}).$$

These quantities shown in Fig. 2c are used as explanatory variables for the interglacial intensities. In addition to these, we consider the average, $I_{\text{AV}}=\frac{\mathrm{1}}{\mathrm{2}}(I_{\mathrm{N}}+I_{\mathrm{S}})$, as an alternative explanatory variable. We then include the previous glacial intensity δ¹⁸O_max as another explanatory variable (Fig. 2d), based on the observation that strong interglacials often follow strong glacials (Lang and Wolff, 2011). Figure 3 shows that these variables are able to separate the strongest interglacials from the rest. In the present framework, the timing and intensity of glacial maxima (t_s and δ¹⁸O_max) are taken from the palaeodata, while the end of complete deglaciations, t_e, is known from a predictive model (Tzedakis et al., 2017) that indicates which insolation cycles lead to complete deglaciations.

We perform a linear regression of δ¹⁸O_min with equation ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\mathrm{N}}+{\mathit{\beta }}_{\mathrm{3}}{I}_{\mathrm{S}}$. We refer to this case in which β₂ and β₃ are chosen independently as Model 2. Model 1 does not use the term in I_S (i.e. β₃≡0). We also consider Model 3 defined by the regression equation ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\text{AV}}$, which is a special case of Model 2, in which the coefficients of I_N and I_S are equal. In all three models, the actual and predicted δ¹⁸O_min are strongly correlated (Fig. 2e and Table 1). Model 2, which has an extra parameter, gives the highest correlation of 0.94 (more variability explained), while Model 3 using the average insolation integral I_AV provides a similar level of correlation, 0.92, to Model 2. In all three models, all the explanatory variables are significant (Table 1). According to the BIC (Sect. 2.3), Model 2 and Model 3 have strong evidence against Model 1 (ΔBIC>8), indicating that insolations in both hemispheres are crucial for the differences in interglacial intensities. On the other hand, the evidence of Model 2 against Model 3 is weak (ΔBIC<2); that is, they are almost equally supported by the data.

Figure 4 Modelling glacial intensities. (a) LR04 δ¹⁸O. The red circles indicate δ¹⁸O_min at each interglacial and the blue triangles the maxima δ¹⁸O_min at glacials. The time intervals between them are shaded. Note that the data are plotted inversely to Fig. 2a, with glacial maxima above interglacial minima. (b) Caloric summer insolation at 65^∘ N. The grey shading is the same as in panel (a). (c) δ¹⁸O_min for each interglacial. (d) Total time L during which the caloric summer insolation is below a threshold 5.7 GJ m⁻² between the interglacial peak and the glacial peak. (e) Time span T between the interglacial peak and the glacial peak. (f) $\mathrm{1}-e^{-T/\mathrm{25}}$. (g) Prediction of δ¹⁸O_max from the linear regression relation with explanatory variables in panels (d) and (f) (R=0.90).

Predictions of all the regression models are consistent with the observation that MIS11c, 9e, 5e, and 1 are stronger than other interglacials over the last 800 kyr (Fig. 2e). Predictions of MIS15e and MIS11c diverge from observations in Model 1 but are well-reproduced in Models 2 and 3, which include I_S. We also tested simpler models without the δ¹⁸O_max term, but they give poor fits, particularly for the strength of MIS7e, 7c, and 5e compared to the original models (Fig. S1 and Table S1 in the Supplement). The contribution of δ¹⁸O_max is also evident from the BIC (the BICs for the original models are substantially lower than those for the simpler models). In Appendix A, we also explored the effect of using the first crossing point between F_N(t) and $\stackrel{\mathrm{‾}}{F}$ for the starting time t_s of Termination III on our model results: the predicted interglacial intensity δ¹⁸O_min of MIS7e is stronger than the observation (Fig. A1) and the prediction skills slightly decrease (Table A1), but the main features of interglacial intensities are reproduced in Models 2 and 3. Taken together, our results indicate that interglacial intensity is determined by the strength of the preceding glacial and the summer energy received during deglaciation in both polar regions.

3.2 Glacial intensities

We now attempt to explain glacial intensity δ¹⁸O_max on the premise that we know the timing of the glacial maximum t_max, the timing of the previous interglacial peak t_min, and its intensity δ¹⁸O_min from observations.

Inspection of Fig. 4a shows that a longer glacial duration ($T=t_{\text{max}}-t_{\text{min}}$) typically results in a larger increase in δ¹⁸O (larger ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}$). For example, the two weakest glacials, MIS15b and MIS7b, are also the shortest, while strong glacials are longer (T ≳ 60 kyr) (Fig. 4e). However, the relation between ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}$ and T is nonlinear. Figure 5a shows the time evolution of δ¹⁸O during each glacial period relative to the time elapsed since the previous interglacial peak δ¹⁸O_min. While the time evolution varies with each glacial, the rate of increase in δ¹⁸O is typically higher in the beginning (T ≲ 20 kyr) (reflecting in large part rapid deep-ocean cooling, Elderfield et al., 2012) and decreases as time elapses (or as δ¹⁸O increases). As a zeroth-order approximation, the increase in δ¹⁸O may be expressed by a linear relaxation process $\mathrm{d}x/\mathrm{d}t=(A-x)/\mathit{\tau }$, x(0)=0. Its solution $x\left(t\right)=A(\mathrm{1}-e^{-t/\mathit{\tau }})$ roughly fits the observed δ¹⁸O changes in Fig. 5a (magenta solid line) for τ∼25 kyr and A∼1.3 ‰ (the same functional form is assumed for an ice sheet response in Imbrie et al., 1993). In the following model, x(T) represents the basic time dependence of δ¹⁸O change in a glacial period with duration T (kyr).

Figure 5 (a) Time evolution of δ¹⁸O during each glacial period from its previous interglacial δ¹⁸O_min. The magenta lines are baseline increase profiles $x\left(t\right)=A(\mathrm{1}-e^{-t/\mathit{\tau }})$ for $\mathit{\tau }=\mathrm{10},\mathrm{25}$ and 50 kyr and A=1.3 ‰. (b) Growth rate of δ¹⁸O ($\mathrm{d}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/\mathrm{d}t$) vs. the caloric summer insolation at 65^∘ N only during a glacial period [t_min, t_max]. $\mathrm{d}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/\mathrm{d}t$ is calculated after 15-point Gaussian smoothing on 1 kyr re-sampled data. The magenta circle with a 1-sigma range shows the mean in a bin of the caloric summer insolation of size 0.1 GJ m⁻².

Table 2 Coefficients and statistics of the regression models for glacial intensities (corresponding to Fig. 4). The model is given as ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}(\mathrm{1}-e^{-T/\mathrm{25}})+{\mathit{\beta }}_{\mathrm{2}}L$. The overall F test provides a p value of less than 0.001 in each model, which rejects the null hypothesis that none of the variables in the model is significant. The asterisks indicate the significance of each coefficient: ^∗ for $p\in (\mathrm{0.01},\mathrm{0.05}]$, ${}^{\ast \ast }$ for $p\in (\mathrm{0.001},\mathrm{0.01}]$, and ${}^{\ast \ast \ast }$ for $p\in [\mathrm{0},\mathrm{0.001}]$.

We then consider the effect of insolation. As observed in Fig. 5b, the growth rate of δ¹⁸O ($\mathrm{d}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/\mathrm{d}t$) is high when the caloric summer insolation at 65^∘ N is low, while the average growth rate is close to zero for high insolation values. Specifically, δ¹⁸O almost exclusively grows (at ∼10 kyr timescales) for values of insolation below ∼5.7–5.8 GJ m⁻². We therefore introduce a measure for low insolation periods (Fig. 4d): the total time L (kyr) within a glacial period [t_min, t_max] during which the caloric summer insolation at 65^∘ N is below an empirical threshold 5.7 GJ m⁻² (see below for this specific value). We obtain a similar result if we use the integral of the insolation below a threshold (Fig. S2) instead of the total time, while the model fit using the total time L is slightly better than this alternative. In either case, a measure of insolation deficit contributes to the prediction of glacial intensities. Thresholding is the simplest way of capturing insolation deficits.

We assume that the δ¹⁸O increase between an interglacial minimum and the ensuing glacial maximum (${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}$) is decomposed into the basic time dependence $\mathrm{1}-e^{-T/\mathrm{25}}$ and the total low-insolation period L (kyr) during the period of glaciation (Fig. 4). That is, we suppose the relation ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}(\mathrm{1}-e^{-T/\mathrm{25}})+{\mathit{\beta }}_{\mathrm{2}}L$. The result of the regression analysis is given in Table 2. The estimated intercept β₀ is quite small −0.033 ‰, and the null hypothesis β₀=0 cannot be rejected (p=0.88). If this model is compared to the one without an intercept, i.e. β₀≡0 (Table 2), the BIC positively supports the latter (ΔBIC>2). Therefore, we select the model without an intercept, β₀, that is, ${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}-{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}=\mathrm{1.4}(\mathrm{1}-e^{-T/\mathrm{25}})+\mathrm{0.025}L$. The predicted δ¹⁸O_max, shown in Fig. 4g, is strongly correlated with observations (R=0.90). The p values for each coefficient are both less than 0.005, suggesting that both explanatory variables are important. What emerges is that the length of the glacial between the interglacial peak and the glacial maximum, the total time during which insolation is below a threshold, and the δ¹⁸O_min of the preceding interglacial all have a large impact on the intensity of the ensuing glacial maximum.

Figure 6 (a) Contributions of different factors to the predicted δ¹⁸O_min anomaly estimated by Model 3 using δ¹⁸O_max and I_AV. The height of each bar indicates the anomaly of each term from its average over 11 interglacials. The labels on the top of the bars show corresponding Marine Isotope Stages (MISs) for δ¹⁸O_min. The labels below the bars show those for δ¹⁸O_max. (b) Contributions of different factors to the predicted δ¹⁸O_max anomaly: the previous interglacial peak value δ¹⁸O_min in the actual data (yellow), the contribution of the basic time dependence, $\mathrm{1.4}(\mathrm{1}-e^{-t/\mathrm{25}})$ (green), and the contribution of the total time of the low insolation spells during the glacial period, 0.025L (dark blue). The height of each bar indicates the anomaly of each term from its average over 11 glacials. The labels on the top of the bars show corresponding MISs for δ¹⁸O_max. The labels below the bars show those for δ¹⁸O_min.

Figure 7 (a) LR04 benthic δ¹⁸O stack (Lisiecki and Raymo, 2005). (b) Obliquity (Laskar et al., 2004). (c) Shackleton05 composite benthic δ¹⁸O record from the eastern equatorial Pacific (Tzedakis et al., 2017).

The final model has four parameters: regression coefficients β_1,2, time constant τ (=25 kyr), and the threshold defining the low-insolation period (5.7 GJ m⁻²). The latter two values are suggested from the observations (Figs. 4 and 5) and are adopted because they provide a good prediction. The result is however relatively insensitive to those parameters; we obtain a good prediction with a correlation coefficient close to 0.9 in a range of τ and threshold values, as shown in Fig. S3. The model with four parameters might be considered complex for predicting 11 glacial intensities. If the model is severely overfitted to the data, the model would not possess prediction ability for new data. In Appendix B, we have shown by the leave-one-out cross-validation method that the model actually has prediction ability for unseen data. It should also be noted that, in this model, the glacial duration T is taken from the LR04 benthic δ¹⁸O record assuming that its orbital tuning is right.

4 Summary and discussion

We have inspected the predictability of interglacial and glacial intensities over the last 800 kyr. Interglacial strength can be well-predicted at the previous glacial maximum with the caveat that, in two cases, the amplitude was achieved only in the second or third insolation peak. Glacial strength is well-predicted at the previous interglacial peak, with the caveat that the length of the glacial is currently taken from the data. While the models contain three to four parameters, they are still simple explanatory frameworks. These models show that interglacial intensity over the last 800 kyr can be described as a function of the strength of the previous glacial and the summer insolation at high latitudes in both hemispheres during the deglaciation and that glacial intensity is linked to the strength of the previous interglacial, the time elapsed from it, and the evolution of boreal summer insolation.

While previous studies (e.g. Lang and Wolff, 2011) have underlined the influence of large ice sheets on the extent of deglaciation, our analysis provides a quantitative description of its contribution. Figure 6a shows the decomposition of the predicted δ¹⁸O_min anomaly into the contribution of different factors: the weak MIS15e deglaciation is attributed to the low average insolation integral I_AV, even though it follows one of the strongest glaciations, MIS16; on the other hand, the weak deglaciations of MIS15a, 13a, 7c, and 7e are attributed to the weak MIS15b, 14, 7d, and 8 glacials respectively. Large ice sheets are more unstable as a result of ice sheet physics, glacio-isostatic adjustments, extensions to lower latitude, and changes in ice sheet albedo (MacAyeal, 1979; Birchfield et al., 1981; Marshall and Clark, 2002; Ganopolski et al., 2010; Abe-Ouchi et al., 2013) and are therefore more sensitive to insolation increases. In addition, the effect of glacial intensity on deglaciation might be operating partly through its influence on atmospheric CO₂. Larger NH ice sheets have a capacity to produce larger amounts of freshwater or a longer period of freshwater discharge into the North Atlantic, weakening the Atlantic Meridional Overturning Circulation and leading to activation of the bipolar seesaw (Knutti et al., 2004; Denton et al., 2010; Ganopolski and Brovkin, 2017; Tzedakis et al., 2022). Large ice sheets may also promote stronger deep-ocean salinity stratification, which stabilizes relatively warm waters in the glacial deep ocean and then amplifies the rate of Antarctic warming during the activation of the bipolar seesaw of ensuing deglaciation (Knorr et al., 2021). These changes in turn lead to more CO₂ outgassing from the Southern Ocean (Watson and Garabato, 2006; Stephens and Keeling, 2000), accelerating ice sheet ablation, and are consistent with the combined role of insolation and CO₂ proposed by Yin and Berger (2012).

Our analysis also developed an empirical relation linking glacial intensity to the length of the preceding glacial period, the total length of low-insolation periods during that glacial, and the preceding interglacial intensity. The decomposition of the predicted δ¹⁸O_max anomaly into the contributions from these factors (Fig. 6b) shows that all pre-MBE glacials (and also MIS7d and 6) were strengthened by the extent of the δ¹⁸O_min. This represents remnant ice at the end of the preceding interglacial that contributes by reducing albedo and as a seed for ice growth that occurs during periods of low insolation. The decomposition also reveals that MIS12, the largest Quaternary glaciation, is the only occasion when all three factors make a positive contribution.

The MBE. A striking feature of our results is the change in the time integrals of the averaged boreal and austral insolation, I_AV, across the MBE (Fig. 2c). The decomposition of the predicted δ¹⁸O_min anomaly (Fig. 6a) clearly shows that while δ¹⁸O_max has positive contributions before and after the MBE, a shift from negative to positive contributions is observed in I_AV after the MBE. As previously discussed, I_AV is dominated by obliquity, which has a strong influence on the insolation received over the local summer season and over the year at high latitudes of both hemispheres (Milankovitch, 1941). The influence of obliquity on high latitudes is further amplified by sea ice and snow albedo feedbacks, while the contribution of precession may be stronger at northern high latitudes (Yin and Berger, 2012). The EPICA Dome C sea salt Na flux record, a proxy for sea ice extent (Wolff et al., 2006), shows a reduction in interglacial sea ice after the MBE (Fig. 1e), which may have allowed more CO₂ outgassing from the Southern Ocean (Watson and Garabato, 2006; Stephens and Keeling, 2000), leading to increased global annual mean temperatures and ice ablation. In addition, the stronger post-MBE component of I_AV and also I_S during deglaciations, enhancing the melting of sea ice and warming of the Southern Ocean, may have destabilized the AIS, increasing its contribution to sea-level highstands. This is consistent with offshore sediment and geochemical data on provenance changes, suggesting increased ice loss from the Wilkes Subglacial Basin during post-MBE interglacials (MIS11c, 5e, and 9e), although the pre-MBE record remains weakly constrained (Wilson et al., 2021). Recently, Mitsui and Boers (2022) developed an artificial neural network (ANN) model that performs a skilful 21 kyr ahead prediction of δ¹⁸O on the basis of the past δ¹⁸O history and the insolation evolution. Through the sensitivity analysis of the ANN model, they concluded that the intensification of interglacials across the MBE is attributed to the amplitude increase in the obliquity forcing. While this is consistent with our conclusions, our present regression model is more physically interpretable than the ANN model and even more precise in predicting δ¹⁸O_min.

Thus, the shift in interglacial intensities at the MBE may be ultimately related to the amplitude modulation of obliquity with a duration of ∼1.2 Myr (Lourens and Hilgen, 1997), which led to higher obliquity variations after 430 kyr BP. If this conjecture is correct, then a similar shift from weaker to stronger interglacials should have occurred about 1.6 Myr BP. The LR04 stack (Fig. 7) hints at an interval of weaker interglacials ∼1.9–1.6 Myr BP, but the averaging of several records to create the stack tends to smooth δ¹⁸O variability. By comparison, inspection of the Shackleton05 composite benthic δ¹⁸O record from the eastern equatorial Pacific (see Tzedakis et al., 2017, for references) shows more clearly the occurrence of weaker interglacials from around 1.9 to 1.6 Myr BP and a shift to stronger interglacials after that. Although it would be interesting to explore whether our model can reproduce a shift in interglacial intensities at ∼1.6 Myr BP, uncertainties about the size of ice sheets and the deep-water temperature components of δ¹⁸O complicate such an undertaking at present.

While we remain some distance from a fully predictive model of temperature, ice volume, and sea level over the entire sequence of glacial cycles, our analysis lays out some of the key predictors that need to be understood in physical models and coupled together and underlines the importance of ice volume legacies and the time course of insolation on the amplitude of glacial cycles.

Appendix A: Sensitivity to the starting time of the termination

In the models of interglacial intensities in Sect. 3.1, the starting time t_s of the termination is chosen as the first crossing point between F_N(t) and $\stackrel{\mathrm{‾}}{F}$, except for Termination III (Fig. 2). Termination III started after the local insolation minimum at 254 kyr BP (orange dotted line) behind the crossing point by about 11 kyr, responding to the second rise in the insolation. This is why we set t_s=254 kyr BP for Termination III in our main models in Fig. 2. Here we show the effect of using the first crossing point between F_N(t) and $\stackrel{\mathrm{‾}}{F}$ for the starting time of Termination III on our model results: the predicted interglacial intensity δ¹⁸O_min of MIS7e is stronger than the observation (Fig. A1) and the prediction skills slightly decrease (Table A1), but the main features of interglacial intensities are reproduced in Models 2 and 3.

Figure A1 Same as Fig. 2 but without the exception for the start timing of Termination III leading to MIS7e. The final correlation with the observed δ¹⁸O_min (red circle) decreases from 0.83 to 0.81 in Model 1, from 0.94 to 0.86 in Model 2, and from 0.92 to 0.83 in Model 3 (see also Table A1). However, the main features of interglacial intensities are reproduced in Models 2 and 3. (a) LR04 δ¹⁸O. The red circles indicate the minima of δ¹⁸O (δ¹⁸O_min) at each interglacial and the blue triangles the maxima (δ¹⁸O_max) at glacials. See below for the grey strips and the dashed lines. (b) Caloric summer half-year insolation at 65^∘ N (F_N, black) and 65^∘ S (F_S, green). The average of the two (magenta) is also shown. The blue dashed lines show timings t_s at which the caloric summer half-year insolation at 65^∘ N exceeds the average 5.845 GJ m⁻² (black horizontal line), and the red dashed lines show timings t_e at which the insolation falls back below the average. Each termination starts roughly around t_s, and it is completed around t_e. The grey strips show the termination intervals [t_s, t_e] based on the insolation curve. (c) Integral of the caloric summer insolation anomaly between t_s and t_e at 65^∘ N, I_N (black cross), the integral at 65^∘ S for the same period, I_S (green diamond), and the average $I_{\text{AV}}=\frac{\mathrm{1}}{\mathrm{2}}(I_{\mathrm{N}}+I_{\mathrm{S}})$ (magenta square). (d) δ¹⁸O_max. (e) Predictions by linear regression models with explanatory variables in panels (c) and (d): Model 1 with I_N (black cross); Model 2 with both I_N and I_S with their own coefficients (blue diamond with cross); Model 3 with I_AV (magenta squares).

Table A1 Coefficients and statistics of the regression models for the data corresponding to Fig. A1 (without the exception for the start timing of Termination III leading to MIS7e). Model 1 (${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\mathrm{N}}$), Model 2 (${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\mathrm{N}}+{\mathit{\beta }}_{\mathrm{3}}{I}_{\mathrm{S}}$), and Model 3 (${{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{min}}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{max}}+{\mathit{\beta }}_{\mathrm{2}}{I}_{\text{AV}}$). The overall F test provides a p value of less than 0.05 in each model, which rejects the null hypothesis that none of the variables in the model is significant. The asterisks indicate the significance of each coefficient: ^∗ for $p\in (\mathrm{0.01},\mathrm{0.05}]$, ${}^{\ast \ast }$ for $p\in (\mathrm{0.001},\mathrm{0.01}]$, and ${}^{\ast \ast \ast }$ for $p\in [\mathrm{0},\mathrm{0.001}]$.

Appendix B: Prediction ability of the models for unseen data

Models that contain more unknown parameters than can be justified by the data are called overfitted (Everitt and Skrondal, 2010). Overfitted models can have low prediction ability for new data, even if they appear skilful for training data. The leave-one-out cross-validation (LOOCV) approach is one of the methods to assess the prediction ability of a machine learning model for new data: a model is trained with N−1 data points, removing one data point from the entire dataset with N points (here N=11). The trained model is then used to predict the removed data point. This procedure is applied for every data point. Based on the average prediction error in LOCCV, we can infer the prediction ability of the model for unseen data. If a regression model is severely overfitted to data, the correlation coefficient between the data and the prediction in LOCCV would substantially decrease from the correlation coefficient obtained by the usual model fitting.

Model 1 for the interglacial intensities gives a correlation coefficient of R=0.70 in LOOCV, while the correlation coefficient in the usual fitting is R=0.83 (Table 1). Model 2 for the interglacial intensities gives R=0.82 in LOOCV and R=0.94 in the usual fitting (Table 1). Model 3 for the interglacial intensities gives R=0.85 in LOOCV and R=0.92 in the usual fitting (Table 1). The model for glacial intensities gives R=0.84 in LOOCV and R=0.90 in the usual fitting. In summary, the correlation in LOOCV is slightly lower than the correlation in the usual fitting, but the difference is not substantially large in each model. That is, our models with three to four parameters have prediction ability for unseen data, even when trained with 10 data points. Therefore we consider that the models are not severely overfitted. Although the number of parameters (three to four) is large in comparison to the number of data points N=11, the models are reasonably simple compared to the complexity of ice age cycles arising from various feedbacks.