Articles | Volume 18, issue 9
Clim. Past, 18, 1983–1996, 2022
https://doi.org/10.5194/cp-18-1983-2022
Clim. Past, 18, 1983–1996, 2022
https://doi.org/10.5194/cp-18-1983-2022
Research article
01 Sep 2022
Research article | 01 Sep 2022

Insolation evolution and ice volume legacies determine interglacial and glacial intensity

Insolation evolution and ice volume legacies determine interglacial and glacial intensity
Takahito Mitsui1,2, Polychronis C. Tzedakis3, and Eric W. Wolff4 Takahito Mitsui et al.
  • 1Earth System Modelling, School of Engineering & Design, Technical University of Munich, Munich, Germany
  • 2Potsdam Institute for Climate Impact Research, Potsdam, Germany
  • 3Environmental Change Research Centre, Department of Geography, University College London, London, UK
  • 4Department of Earth Sciences, University of Cambridge, Cambridge, UK

Correspondence: Takahito Mitsui (takahito321@gmail.com)

Abstract

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.

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f01

Figure 1Changes 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 CO2 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 δ18O 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 (δ18O) 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 δ18O 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 δ18O, atmospheric CO2 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 CO2 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 δ18O 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 δ18O 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 δ18O from the insolation, the only external driver of the climate system. Atmospheric CO2 or other climate feedbacks are considered agents in between insolation and δ18O 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 δ18O of seawater (e.g. Elderfield et al., 2012) exceed differences in interglacial highstands (Dutton et al., 2015), we use the LR04 benthic δ18O stack record (Lisiecki and Raymo, 2005) (mean standard error of 0.06 ‰) as an integrated metric of interglacial intensities. Although benthic δ18O 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 δ18O records, therefore, reflects the slow ice sheet buildup and rapid decay. By extension, differences in δ18Omax between glacial maxima in the LR04 record largely reflect differences in ice volume. We assume that the orbital tuning in the LR04 δ18O 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 δ18O 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).

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f02

Figure 2Modelling interglacial intensities. (a) LR04 δ18O. The red circles indicate the minima of δ18O (δ18Omin) at each interglacial and the blue triangles the maxima (δ18Omax) at glacials. See below for the grey strips and the dashed lines. (b) Caloric summer half-year insolation at 65 N (FN, black) and 65 S (FS, green). The average of the two (magenta) is also shown. The blue dashed lines show timings ts at which the caloric summer half-year insolation at 65 N exceeds the average 5845 GJ m−2 (black horizontal line) and the red dashed lines show timings te at which the insolation falls back below the average. Each termination starts roughly around ts, and it is completed around te. 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 [ts, te] based on the insolation curve. (c) Integral of the caloric summer insolation anomaly between ts and te at 65 N, IN (black cross), the integral at 65 S for the same period, IS (green diamond), and the average IAV=12(IN+IS) (magenta square). (d) δ18Omax. (e) Predictions by linear regression models with explanatory variables in panels (c) and (d): Model 1 with IN (black cross); Model 2 with both IN and IS with their own coefficients (blue diamond with cross); Model 3 with IAV (magenta squares).

Download

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f03

Figure 3Scatter plot in the space: (a) (δ18Omax, IN, δ18Omin), (b) (δ18Omax, IAV, δ18Omin), and (c) (δ18Omax, β2IN+β3IS, δ18Omin). (a) We can separate the four strong interglacials (MIS11c, 9e, 5e, and 1) from the weak interglacials in the plane of explanatory variables (δ18Omax, IN). In this case, however, the weak interglacial MIS15e is close to the strong interglacials in the δ18Omax-IN plane. Replacing IN with a mixture of IN and IS (b, c) leads to a clear separation of the four strongest interglacials (MIS11c, 9e, 5e, and 1) from the rest.

Download

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 F=5.845GJ m−2 over the last 1×106 years (Myr). The solar constant used is 1368 W m−2.

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 (ΔBIC=BICk-BICj) 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 0<ΔBIC<2, positive if 2<ΔBIC<6, strong if 6<ΔBIC<10, and very strong if ΔBIC>10.

Table 1Coefficients and statistics of the regression models for interglacial intensity δ18Omin (corresponding to Fig. 2). Model 1 (δ18Omin=β0+β1δ18Omax+β2IN), Model 2 (δ18Omin=β0+β1δ18Omax+β2IN+β3IS), and Model 3 (δ18Omin=β0+β1δ18Omax+β2IAV). 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(0.01,0.05], for p(0.001,0.01], and for p[0,0.001].

Download Print Version | Download XLSX

3 Results

3.1 Interglacial intensities

Figure 2a shows the LR04 benthic δ18O record over the last 800 kyr, where interglacial peaks (δ18Omin) are indicated by red circles and glacial peaks (δ18Omax) by blue triangles. In Fig. 2b, the caloric summer half-year insolation at 65 N (FN(t), black) and that at 65 S (FS(t), green) are shown respectively. The average of the two insolations (12(FN+FS), 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 FN(t) exceeds its average F=5.845GJ m−2. Exceptionally, the start of Termination III leading to MIS7e is delayed relative to the crossing point between FN(t) and F by about 11 kyr, possibly because the rise in the insolation FN(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 FN(t) and F associated with terminations by ts (including an exception ts=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 (te) when the insolation FN(t) falls below the average 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 δ18Omin. We postulate that the magnitude of each termination is related to the boreal summer insolation anomaly integrated during the insolation-based termination period [ts, te] at 65 N and the austral summer insolation (at 65 S) anomaly integrated during the same period:

IN=t=tste(FN(t)-F),IS=t=tste(FS(t)-F).

These quantities shown in Fig. 2c are used as explanatory variables for the interglacial intensities. In addition to these, we consider the average, IAV=12(IN+IS), as an alternative explanatory variable. We then include the previous glacial intensity δ18Omax 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 (ts and δ18Omax) are taken from the palaeodata, while the end of complete deglaciations, te, 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 δ18Omin with equation δ18Omin=β0+β1δ18Omax+β2IN+β3IS. We refer to this case in which β2 and β3 are chosen independently as Model 2. Model 1 does not use the term in IS (i.e. β3≡0). We also consider Model 3 defined by the regression equation δ18Omin=β0+β1δ18Omax+β2IAV, which is a special case of Model 2, in which the coefficients of IN and IS are equal. In all three models, the actual and predicted δ18Omin 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 IAV 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.

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f04

Figure 4Modelling glacial intensities. (a) LR04 δ18O. The red circles indicate δ18Omin at each interglacial and the blue triangles the maxima δ18Omin 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) δ18Omin for each interglacial. (d) Total time L during which the caloric summer insolation is below a threshold 5.7 GJ m−2 between the interglacial peak and the glacial peak. (e) Time span T between the interglacial peak and the glacial peak. (f) 1-e-T/25. (g) Prediction of δ18Omax from the linear regression relation with explanatory variables in panels (d) and (f) (R=0.90).

Download

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 IS. We also tested simpler models without the δ18Omax 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 δ18Omax 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 FN(t) and F for the starting time ts of Termination III on our model results: the predicted interglacial intensity δ18Omin 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 δ18Omax on the premise that we know the timing of the glacial maximum tmax, the timing of the previous interglacial peak tmin, and its intensity δ18Omin from observations.

Inspection of Fig. 4a shows that a longer glacial duration (T=tmax-tmin) typically results in a larger increase in δ18O (larger δ18Omax-δ18Omin). 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 δ18Omax-δ18Omin and T is nonlinear. Figure 5a shows the time evolution of δ18O during each glacial period relative to the time elapsed since the previous interglacial peak δ18Omin. While the time evolution varies with each glacial, the rate of increase in δ18O 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 δ18O increases). As a zeroth-order approximation, the increase in δ18O may be expressed by a linear relaxation process dx/dt=(A-x)/τ, x(0)=0. Its solution x(t)=A(1-e-t/τ) roughly fits the observed δ18O 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 δ18O change in a glacial period with duration T (kyr).

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f05

Figure 5(a) Time evolution of δ18O during each glacial period from its previous interglacial δ18Omin. The magenta lines are baseline increase profiles x(t)=A(1-e-t/τ) for τ=10,25 and 50 kyr and A=1.3 ‰. (b) Growth rate of δ18O (dδ18O/dt) vs. the caloric summer insolation at 65 N only during a glacial period [tmin, tmax]. dδ18O/dt 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−2.

Download

Table 2Coefficients and statistics of the regression models for glacial intensities (corresponding to Fig. 4). The model is given as δ18Omax-δ18Omin=β0+β1(1-e-T/25)+β2L. 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(0.01,0.05], for p(0.001,0.01], and for p[0,0.001].

Download Print Version | Download XLSX

We then consider the effect of insolation. As observed in Fig. 5b, the growth rate of δ18O (dδ18O/dt) 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, δ18O almost exclusively grows (at ∼10 kyr timescales) for values of insolation below ∼5.7–5.8 GJ m−2. We therefore introduce a measure for low insolation periods (Fig. 4d): the total time L (kyr) within a glacial period [tmin, tmax] during which the caloric summer insolation at 65 N is below an empirical threshold 5.7GJ m−2 (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 δ18O increase between an interglacial minimum and the ensuing glacial maximum (δ18Omin-δ18Omax) is decomposed into the basic time dependence 1-e-T/25 and the total low-insolation period L (kyr) during the period of glaciation (Fig. 4). That is, we suppose the relation δ18Omax-δ18Omin=β0+β1(1-e-T/25)+β2L. The result of the regression analysis is given in Table 2. The estimated intercept β0 is quite small 0.033 ‰, and the null hypothesis β0=0 cannot be rejected (p=0.88). If this model is compared to the one without an intercept, i.e. β0≡0 (Table 2), the BIC positively supports the latter (ΔBIC>2). Therefore, we select the model without an intercept, β0, that is, δ18Omax-δ18Omin=1.4(1-e-T/25)+0.025L. The predicted δ18Omax, 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 δ18Omin of the preceding interglacial all have a large impact on the intensity of the ensuing glacial maximum.

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f06

Figure 6(a) Contributions of different factors to the predicted δ18Omin anomaly estimated by Model 3 using δ18Omax and IAV. 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 δ18Omin. The labels below the bars show those for δ18Omax. (b) Contributions of different factors to the predicted δ18Omax anomaly: the previous interglacial peak value δ18Omin in the actual data (yellow), the contribution of the basic time dependence, 1.4(1-e-t/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 δ18Omax. The labels below the bars show those for δ18Omin.

Download

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f07

Figure 7(a) LR04 benthic δ18O stack (Lisiecki and Raymo, 2005). (b) Obliquity (Laskar et al., 2004). (c) Shackleton05 composite benthic δ18O 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−2). 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 δ18O 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 δ18Omin anomaly into the contribution of different factors: the weak MIS15e deglaciation is attributed to the low average insolation integral IAV, 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 CO2. 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 CO2 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 CO2 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 δ18Omax 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 δ18Omin. 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, IAV, across the MBE (Fig. 2c). The decomposition of the predicted δ18Omin anomaly (Fig. 6a) clearly shows that while δ18Omax has positive contributions before and after the MBE, a shift from negative to positive contributions is observed in IAV after the MBE. As previously discussed, IAV 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 CO2 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 IAV and also IS 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 δ18O on the basis of the past δ18O 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 δ18Omin.

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 δ18O variability. By comparison, inspection of the Shackleton05 composite benthic δ18O 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 δ18O 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 ts of the termination is chosen as the first crossing point between FN(t) and 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 ts=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 FN(t) and F for the starting time of Termination III on our model results: the predicted interglacial intensity δ18Omin 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.

https://cp.copernicus.org/articles/18/1983/2022/cp-18-1983-2022-f08

Figure A1Same as Fig. 2 but without the exception for the start timing of Termination III leading to MIS7e. The final correlation with the observed δ18Omin (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 δ18O. The red circles indicate the minima of δ18O (δ18Omin) at each interglacial and the blue triangles the maxima (δ18Omax) at glacials. See below for the grey strips and the dashed lines. (b) Caloric summer half-year insolation at 65 N (FN, black) and 65 S (FS, green). The average of the two (magenta) is also shown. The blue dashed lines show timings ts at which the caloric summer half-year insolation at 65 N exceeds the average 5.845 GJ m−2 (black horizontal line), and the red dashed lines show timings te at which the insolation falls back below the average. Each termination starts roughly around ts, and it is completed around te. The grey strips show the termination intervals [ts, te] based on the insolation curve. (c) Integral of the caloric summer insolation anomaly between ts and te at 65 N, IN (black cross), the integral at 65 S for the same period, IS (green diamond), and the average IAV=12(IN+IS) (magenta square). (d) δ18Omax. (e) Predictions by linear regression models with explanatory variables in panels (c) and (d): Model 1 with IN (black cross); Model 2 with both IN and IS with their own coefficients (blue diamond with cross); Model 3 with IAV (magenta squares).

Download

Table A1Coefficients 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 (δ18Omin=β0+β1δ18Omax+β2IN), Model 2 (δ18Omin=β0+β1δ18Omax+β2IN+β3IS), and Model 3 (δ18Omin=β0+β1δ18Omax+β2IAV). 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(0.01,0.05], for p(0.001,0.01], and for p[0,0.001].

Download Print Version | Download XLSX

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.

Code and data availability

R package palinsol version 0.97 (Crucifix, 2016) used for calculating the caloric summer half-year insolation at 65 N and at 65 S is available from (last access: 26 August 2022). The post-processed data used for Figs. 2 and 4 are provided as the Supplement. The R codes used in this study are available from the corresponding author (Takahito Mitsui) upon request.

Supplement

The post-processed data used for Figs. 2 and 4 are provided as the Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/cp-18-1983-2022-supplement.

Author contributions

PCT and TM conceived the study. TM conducted the analyses with contributions from PCT and EWW. All three authors interpreted and discussed the results and wrote the manuscript.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Climate of the Past. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer

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

Acknowledgements

The authors thank Michel Crucifix, who modified R package palinsol in order to compute the Southern Hemisphere caloric summer half-year insolation. Takahito Mitsui was funded by the Volkswagen Foundation. Polychronis C. Tzedakis was funded by the UK Natural Environment Research Council (NE/V001620/1). Eric W. Wolff is supported by a Royal Society Professorship. This is a contribution to the PAGES Working Group on Quaternary Interglacials (QUIGS).

Financial support

This work was supported by the German Research Foundation (DFG) and the Technical University of Munich (TUM) in the framework of the Open Access Publishing Program.

Review statement

This paper was edited by Erin McClymont and reviewed by two anonymous referees.

References

Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J. I., Takahashi, K., and Blatter, H.: Insolation-driven 100 000 year glacial cycles and hysteresis of ice-sheet volume, Nature, 500, 190–193, 2013. 

Barth, A. M., Clark, P. U., Bill, N. S., He, F., and Pisias, N. G.: Climate evolution across the Mid-Brunhes Transition, Clim. Past, 14, 2071–2087, https://doi.org/10.5194/cp-14-2071-2018, 2018. 

Batchelor, C. L. , Margold, M., Krapp, M., Murton, D. K., Dalton, A. S., Gibbard, P. L., Stokes, C. R. , Murton, J. B., and Manica, A.: The configuration of Northern Hemisphere ice sheets through the Quaternary, Nat. Commun., 10, 1–10, 2019. 

Bauska, T. K., Marcott, S. A., and Brook, E. J.: Abrupt changes in the global carbon cycle during the last glacial period, Nat. Geosci., 14, 91–96, 2021. 

Bazin, L., Landais, A., Lemieux-Dudon, B., Toyé Mahamadou Kele, H., Veres, D., Parrenin, F., Martinerie, P., Ritz, C., Capron, E., Lipenkov, V., Loutre, M.-F., Raynaud, D., Vinther, B., Svensson, A., Rasmussen, S. O., Severi, M., Blunier, T., Leuenberger, M., Fischer, H., Masson-Delmotte, V., Chappellaz, J., and Wolff, E.: An optimized multi-proxy, multi-site Antarctic ice and gas orbital chronology (AICC2012): 120–800 ka, Clim. Past, 9, 1715–1731, https://doi.org/10.5194/cp-9-1715-2013, 2013. 

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J.: Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present, Geophys. Res. Lett., 42, 542–549, 2015. 

Berger, A., Li, X. S., and Loutre, M. F.: Modelling northern hemisphere ice volume over the last 3 Ma, Quaternary Sci. Rev., 18, 1–11, 1999. 

Berger, A. L.: Long-Term Variations of Caloric Insolation Resulting from the Earth's Orbital Elements, Quaternary Res., 9, 139–167, 1978. 

Berger, W. H. and Wefer, G.: On the Dynamics of The Ice Ages: Stage-11 Paradox, Mid-Brunhes Climate Shift, and 100-Ky Cycle, in: Earth's Climate and Orbital Eccentricity: The Marine Isotope Stage 11 Question, edited by: Droxler, A. W., Poore, R. Z., and Burckle, L. H., https://doi.org/10.1029/137GM04, 2003. 

Birchfield, G. E., Weertman, J., and Lunde, A. T.: A paleoclimate model of Northern Hemisphere ice sheets, Quaternary Res., 15, 126–142, 1981. 

Crucifix, M.: Palinsol: insolation for palaeoclimate studies, R package version 0.93, 2016. 

Denton, G. H., Anderson, R. F., Toggweiler, J. R., Edwards, R. L., Schaefer, J. M., and Putnam, A. E.: The last glacial termination, Science, 328, 1652–1656, 2010. 

Dutton, A., Carlson, A. E., Long, A. J., Milne, G. A., Clark, P. U., DeConto, R., Horton, B. P., Rahmstorf, S., and Raymo, M. E.: Sea-level rise due to polar ice-sheet mass loss during past warm periods, Science, 349, aaa4019, https://doi.org/10.1126/science.aaa4019, 2015. 

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. 

EPICA community members: Eight glacial cycles from an Antarctic ice core, Nature, 429, 623–628, 2004. 

Everitt, B. S. and Skrondal, A.: The Cambridge dictionary of statistics, Cambridge University Press, ISBN-13 978-0-521-76699-9, 2010. 

Ganopolski, A. and Brovkin, V.: Simulation of climate, ice sheets and CO2 evolution during the last four glacial cycles with an Earth system model of intermediate complexity, Clim. Past, 13, 1695–1716, https://doi.org/10.5194/cp-13-1695-2017, 2017. 

Ganopolski, A., Calov, R., and Claussen, M.: Simulation of the last glacial cycle with a coupled climate ice-sheet model of intermediate complexity, Clim. Past, 6, 229–244, https://doi.org/10.5194/cp-6-229-2010, 2010. 

Ganopolski, A., Winkelmann, R., and Schellnhuber, H. J.: Critical insolation–CO2 relation for diagnosing past and future glacial inception, Nature, 529, 200–203, 2016. 

Grant, K. M., Rohling, E. J., Ramsey, C. B., Cheng, H., Edwards, R. L., Florindo, F., Heslop, D., Marra, F., Roberts, A. P., Tamisiea, M. E., and Williams, F.: Sea-level variability over five glacial cycles, Nat. Commun., 5, 1–9, 2014. 

Holden, P. B., Edwards, N. R., Wolff, E. W., Valdes, P. J., and Singarayer, J. S.: The Mid-Brunhes event and West Antarctic ice sheet stability, J. Quaternary Sci., 26, 474–477, 2011. 

Hughes, P. D. and Gibbard, P. L.: Global glacier dynamics during 100 ka Pleistocene glacial cycles, Quaternary Res., 90, 222–243, 2018. 

Hughes, P. D., Gibbard, P. L., and Ehlers, J.: The “missing glaciations” of the Middle Pleistocene, Quaternary Res., 96, 161–183, 2020. 

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

Imbrie, J., Berger, A., Boyle, E. A., Clemens, S. C., Duffy, A., Howard, W. R., Kukla, G., Kutzbach, J., Martinson, D. G., McIntyre, A., Mix, A. C., Molfino, B., Morley, J. J., Peterson, L. C., Pisias, N. G., Prell, W. L., Raymo, M. E., Shackleton, N. J., and Toggweiler, J. R.: On the structure and origin of major glaciation cycles 2. The 100 000 year cycle, Paleoceanography, 8, 699–735, 1993. 

Knorr, G., Barker, S., Zhang, X., Lohmann, G., Gong, X., Gierz, P., Stepanek, C., and Stap, L. B.: A salty deep ocean as a prerequisite for glacial termination, Nat. Geosci., 14, 930–936, 2021. 

Knutti, R., Flückiger, J., Stocker, T. F., and Timmermann, A.: Strong hemispheric coupling of glacial climate through freshwater discharge and ocean circulation, Nature, 430, 851–856, 2004. 

Lang, N. and Wolff, E. W.: Interglacial and glacial variability from the last 800 ka in marine, ice and terrestrial archives, Clim. Past, 7, 361–380, https://doi.org/10.5194/cp-7-361-2011, 2011. 

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M., and Levrard, B.: A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285, 2004. 

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records, Paleoceanography, 20, https://doi.org/10.1029/2004PA001071, 2005. 

Lourens, L. J. and Hilgen, F. J.: Long-periodic variations in the Earth's obliquity and their relation to third-order eustatic cycles and late Neogene glaciations, Quatern. Int., 40, 43–52, 1997. 

MacAyeal, D. R.: A catastrophe model of the paleoglimate, J. Glaciol., 24, 245–257, 1979. 

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, 2214, https://doi.org/10.1029/2002GL015192, 2002. 

Milankovitch, M.: Kanon der Erdbestrahlung und seine Anwendung auf des Eiszeitenproblem. Special Publication 132, Section of Mathematical and Natural Sciences, vol. 33, p. 633, Royal Serbian Academy of Sciences, Belgrade, Serbia (“Canon of Insolation and the Ice Age Problem”) (trans. Israel Program for the US Department of Commerce and the National Science Foundation, Washington, D.C., USA, 1969, and by Zavod za udzbenike i nastavna sredstva in cooperation with Muzej nauke i tehnike Srpske akademije nauka i umetnosti, Belgrade, Serbia, 1998), 1941. 

Mitsui, T. and Boers, N.: Machine learning approach reveals strong link between obliquity amplitude increase and the Mid-Brunhes transition, Quaternary Sci. Rev., 277, 107344, https://doi.org/10.1016/j.quascirev.2021.107344, 2022. 

Nehrbass-Ahles, C., Shin J., Schmitt, J., Bereiter, B., Joos, F., Schilt, A., Schmidely, L., Silva, L., Teste, G., Grilli, R., Chappellaz, J., Hodell, D., Fischer, H., and Stocker, T. F.: Abrupt CO2 release to the atmosphere under glacial and early interglacial climate conditions, Science, 369, 1000–1005, 2020. 

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

Parrenin, F. and Paillard, D.: Amplitude and phase of glacial cycles from a conceptual model, Earth Planet. Sc. Lett., 214, 243–250, 2003. 

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

Past Interglacials Working Group of PAGES: Interglacials of the last 800 000 years, Rev. Geophys., 54, 162–219, 2016. 

Obase, T., Abe-Ouchi, A., and Saito, F.: Abrupt climate changes in the last two deglaciations simulated with different Northern ice sheet discharge and insolation, Sci. Rep., 11, 1–11, 2021. 

Raftery, A. E.: Bayesian model selection in social research, Sociol. Methodol., 25, 111–163, 1995. 

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

Shakun, J. D., Lea, D. W., Lisiecki, L. E., and Raymo, M. E.: An 800 kyr record of global surface ocean δ18O and implications for ice volume-temperature coupling, Earth Planet. Sc. Lett., 426, 58–68, 2015. 

Shin, J., Nehrbass-Ahles, C., Grilli, R., Chowdhry Beeman, J., Parrenin, F., Teste, G., Landais, A., Schmidely, L., Silva, L., Schmitt, J., Bereiter, B., Stocker, T. F., Fischer, H., and Chappellaz, J.: Millennial-scale atmospheric CO2 variations during the Marine Isotope Stage 6 period (190–135 ka), Clim. Past, 16, 2203–2219, https://doi.org/10.5194/cp-16-2203-2020, 2020. 

Snyder, C. W.: Evolution of global temperature over the past two million years, Nature, 538, 226–228, 2016. 

Spratt, R. M. and Lisiecki, L. E.: A Late Pleistocene sea level stack, Clim. Past, 12, 1079–1092, https://doi.org/10.5194/cp-12-1079-2016, 2016. 

Stephens, B. B. and Keeling, R. F.: The influence of Antarctic sea ice on glacial–interglacial CO2 variations, Nature, 404, 171–174, 2000. 

Tzedakis, P. C., Raynaud, D., McManus, J. F., Berger, A., Brovkin, V., and Kiefer, T.: Interglacial diversity, Nat. Geosci., 2, 751–755, 2009. 

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

Tzedakis, P. C., Hodell, D. A., Nehrbass-Ahles, C., Mitsui, T., and Wolff, E. W.: Marine Isotope Stage 11c: An unusual interglacial, Quaternary Sci. Rev., 284, 107493, https://doi.org/10.1016/j.quascirev.2022.107493, 2022. 

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

Watson, A. J. and Naveira Garabato, A. C.: The role of Southern Ocean mixing and upwelling in glacial-interglacial atmospheric CO2 change, Tellus B, 58, 73–87, 2006. 

Willeit, M., Ganopolski, A., Calov, R., and Brovkin, V.: Mid-Pleistocene transition in glacial cycles explained by declining CO2 and regolith removal, Sci. Adv., 5, eaav7337, https://doi.org/10.1126/sciadv.aav7337, 2019. 

Wilson, D. J., van de Flierdt, T., McKay, R. M., and Naish, T. R.: Pleistocene Antarctic climate variability: ice sheet, ocean and climate interactions, in: Antarctic Climate Evolution, 2nd edn., edited by: Florindo, F., Siegert, M., de Santis, L., and Naish, T., Elsevier, 523–621, https://doi.org/10.1016/B978-0-12-819109-5.00001-3, 2021. 

Wolff, E. W., Fischer, H., Fundel, F., Ruth, U., Twarloh, B., Littot, G., Mulvaney, R., Röthlisberger, R., de Angelis, M., Boutron, C., Hansson, M., Jonsell, U., Hutterli, M., Bigler, M., Lambert, F., Kaufmann, P., Stauffer, B., Steffensen, J. P., Siggaard-Andersen, M. L., Udisti, R., Becagli, S., Castellano, E., Severi, M., Wagenbach, D., Barbante, C., Gabrielli, P., and Gaspari, V.: Southern Ocean sea-ice extent, productivity and iron flux over the past eight glacial cycles, Nature, 440, 491–496, 2006. 

Yin, Q.: Insolation-induced mid-Brunhes transition in Southern Ocean ventilation and deep-ocean temperature, Nature, 494, 222–225, 2013. 

Yin, Q. Z. and Berger, A.: Individual contribution of insolation and CO2 to the interglacial climates of the past 800 000 years, Clim. Dynam., 38, 709–724, https://doi.org/10.1007/s00382-011-1013-5, 2012. 

Download
Short summary
We provide simple quantitative models for the interglacial and glacial intensities over the last 800 000 years. Our results suggest that the memory of previous climate states and the time course of the insolation in both hemispheres are crucial for understanding interglacial and glacial intensities. In our model, the shift in interglacial intensities at the Mid-Brunhes Event (~430 ka) is ultimately attributed to the amplitude modulation of obliquity.