Predicting trends in atmospheric CO₂ across the Mid-Pleistocene Transition using existing climate archives

During the Mid-Pleistocene Transition (MPT), ca. 1200–800 000 years ago (ka), the Earth's glacial cycles changed from 41 to 100 kyr periodicity. The emergence of this longer ice age periodicity was accompanied by higher global ice volume in glacial periods and lower global ice volume in interglacial periods. Since there is no known change in external orbital forcing across the MPT, it is generally agreed that the cause of this transition is internal to the Earth system. Resolving the climate, carbon cycle, and cryosphere processes responsible for the MPT remains a major challenge in Earth and palaeoclimate science. To address this challenge, the international ice core community has prioritised recovery of an ice core record spanning the MPT interval.

Here we present results from a simple generalised least-squares (GLS) model that predicts atmospheric CO₂ out to 1.8 Myr. Our prediction utilises existing records of atmospheric carbon dioxide (CO₂) from Antarctic ice cores spanning the past 800 kyr along with the existing LR04 benthic $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{calcite}}$ stack (Lisiecki and Raymo, 2005; hereafter “benthic δ¹⁸O stack”) from marine sediment cores. Our predictions assume that the relationship between CO₂ and benthic δ¹⁸O over the past 800 000 years can be extended over the last 1.5 million years. The implicit null hypothesis is that there has been no fundamental change in feedbacks between atmospheric CO₂ and the climate parameters represented by benthic δ¹⁸O, global ice volume, and ocean temperature.

We test the GLS-model-predicted CO₂ concentrations against observed blue ice CO₂ concentrations, δ¹¹B-based CO₂ reconstructions from marine sediment cores, and δ¹³C of leaf-wax-based CO₂ reconstructions (Higgins et al., 2015; Yan et al., 2019; Yamamoto et al., 2022). We show that there is no clear evidence from the existing blue ice or proxy CO₂ data to reject our predictions or our associated null hypothesis. A definitive test and/or rejection of the null hypothesis may be provided following recovery and analysis of continuous oldest ice core records from Antarctica, which are still several years away. The record presented here should provide a useful comparison for the oldest ice core records and an opportunity for further constraints on the processes involved in the MPT.

1 Introduction

Ice core records from Antarctica provide comprehensive and continuous records of many climate parameters over the last 800 000 years, e.g. the ice cores from Vostok (Petit et al., 1999) and the European Project for Ice Coring in Antarctica's Dome-C (EDC) (Jouzel et al., 2007). One of the major challenges in climate science lies beyond the current threshold of the ice core record. The Mid-Pleistocene Transition (MPT) spans from ca. 1200–800 000 years ago (ka) (Chalk et al., 2017) and is characterised by a change from regularly paced 40 000-year (kyr) glacial cycles with thinner glacial ice sheets to quasi-periodic 100 kyr glacial cycles in which ice sheets are more persistent and thicker (Clark et al., 2006; Chalk et al., 2017). To resolve the forcings and feedbacks involved in this transition, multiple nations are targeting recovery of continuous ice cores spanning the MPT under the framework of the International Partnerships in Ice Core Science (IPICS) oldest ice core challenge (IPICS, 2020).

The purpose of the current study is to make a simple prediction of atmospheric CO₂ across the MPT. Cross-comparison of our and other predicted CO₂ records against observed MPT CO₂ data will aid in testing competing hypotheses on the cause of the transition, in particular the role of carbon cycle changes.

The MPT occurred in the absence of any changes to orbital insolation forcing; therefore, the mechanisms behind the MPT must be internal to the Earth system (Raymo, 1997; Ruddiman et al., 1989). Multiple hypotheses have been put forward to explain the transition. A common element in many of these is internal climate–Earth system changes which allow for the development of thicker, more extensive ice sheets that could endure insolation peaks corresponding to the 23 kyr precession and 41 kyr obliquity cycles, i.e. an increase in the threshold for deglaciation and altered sensitivity to orbital forcings (McClymont et al., 2013; Tzedakis et al., 2017). Indeed, the skipped obliquity cycle hypothesis proposes that the 100 kyr signal seen in spectral analysis of the post-MPT benthic δ¹⁸O stack (e.g. Fig. 1a) may be comprised of alternating 80 and 120 kyr signals, i.e. in which the intervening obliquity cycles are skipped. Among the prominent hypotheses to explain an increased threshold for deglaciation are the following three.

1. Long-term decrease in radiative forcing due to a secular reduction in atmospheric CO₂ across the transition (e.g. Berger et al., 1999; Hönisch et al., 2009; Raymo et al., 1988). According to this view, reduced radiative forcing drives the formation of larger and more stable ice sheets.

2. Progressive removal of sub-glacial regolith during the 41 kyr glacial cycles. Clark and Pollard (1998) proposed that ice sheet basal sliding prior to the MPT was enhanced by the presence of a low-friction sedimentary regolith layer between the Laurentide ice sheet and the crystalline bedrock. According to this view, progressive removal of this sedimentary layer then favoured the development of larger and more persistent post-MPT ice sheets.

3. Phase locking of the Northern and Southern Hemisphere ice sheets. In frequency spectra of the global marine benthic δ¹⁸O record (Fig. 1) there is no evidence of the precession (23 kyr) component of Northern Hemisphere insolation prior to the MPT; the spectra are dominated by the obliquity (41 kyr) component (Fig. 1c). Emergence of significant precession and 100 kyr signals occurs across the MPT (Fig. 1b), and all three components are clearly present after the MPT (Fig. 1a). Raymo et al. (2006) suggested that precession-paced changes in Northern and Southern Hemisphere ice volumes may have occurred prior to the MPT but are cancelled due to out-of-phase ice volume changes between the two hemispheres. According to this view, during the MPT the precession-paced changes fall into phase between the two hemispheres such that the precession signal emerges (Raymo et al., 2006). In this view the global synchronisation of ice volume drives the formation of larger and more stable ice sheets.

Figure 1 Thomson multi-taper method (MTM) spectral analysis representing the relative power of signal periodicity for (a) the benthic δ¹⁸O stack after (0–800 ka) the Mid–Pleistocene Transition (MPT), (b) benthic δ¹⁸O across the MPT (800–1200 ka), and (c) benthic δ¹⁸O prior to the onset of the MPT (1200–1500 ka). Each has a robust AR(1) 95 % confidence interval (dashed red line). Benthic δ¹⁸O stack data are from Lisiecki and Raymo (2005).

These hypotheses are not mutually exclusive. For a recent review on the cause of the MPT see Berends et al. (2021b).

For a long-term decrease in radiative forcing by atmospheric CO₂ to be the cause of the MPT, the reduction in CO₂ might be expected in both glacial and interglacial stages (Chalk et al., 2017). However, low-resolution boron-isotope-based CO₂ reconstructions by Hönisch et al. (2009) and Chalk et al. (2017) suggest that glacial-stage CO₂ drawdown occurred over the MPT in the absence of interglacial CO₂ drawdown. Glacial-stage CO₂ drawdown across the MPT may be a positive climate–carbon cycle feedback to changes in ice sheet dynamics, including CO₂ drawdown by enhanced iron fertilisation of the Southern Ocean in response to exposed continental shelves due to lower sea level, as well as planetary drying associated with colder climate conditions (Chalk et al., 2017). Colder glacial temperatures that enhance the solubility of CO₂ in the oceans and reduced abyssal ocean ventilation have also been implicated in enhanced glacial-stage ocean storage of CO₂ (McClymont et al., 2013; Hasenfratz et al., 2019).

Testing of hypotheses on the cause of the MPT is currently limited by the lack of a continuous ice core that spans its duration. The International Partnership in Ice Core Sciences (IPICS) has nominated recovery of such a record as a key priority in ice core research (IPICS, 2020). Multiple national and international projects have commenced, or are soon to commence, drilling for “oldest ice” (see e.g. Shugi, 2022). In this project, we take inspiration from the “EPICA challenge” in which the palaeoclimate and modelling community was challenged to predict the global atmospheric carbon dioxide and methane concentrations from 800–400 ka based on the existing 400 kyr Vostok ice core record (Wolff et al., 2004). Here, we use a generalised least-squares (GLS) model trained on continuous climate archives to predict a CO₂ record out to 1.8 Ma. We utilise two primary datasets for the GLS model: the existing 800 kyr ice core composite record of atmospheric CO₂ (Bereiter et al., 2015) and the LR04 benthic stack of 57 globally distributed records of the ¹⁸O to ¹⁶O ratio of fossil benthic foraminifera calcite (hereafter referred to as the LR04 δ¹⁸O benthic stack). The δ¹⁸O ratios in the LR04 benthic stack are governed primarily by deep-ocean temperature and global ice volume at the time the foraminifera lived, with higher values indicating both increased ice volume and a colder climate. The relationship between the ice volume and ocean temperature components contributing to the δ¹⁸O benthic stack is not linear. Separating the two signals remains challenging and has been attempted elsewhere using a range of approaches from comparison with paired deep-ocean temperature proxies (Elderfield et al., 2012) and inverse modelling (Berends et al., 2021a) to spectral analysis (e.g. Huybers and Wunsch, 2005).

Figure 2 shows a scatterplot of the LR04 δ¹⁸O benthic stack versus observed ice core CO₂ over the past 800 kyr. Both datasets are binned to equivalent 3 kyr time steps (Methods). The Pearson's correlation coefficient (r) between the datasets is −0.82 (p<0.05), indicating that ∼68 % of the variance in observed CO₂ is shared with the LR04 δ¹⁸O benthic stack. This strong relationship provides an initial rationale for using the LR04 δ¹⁸O benthic stack as an input parameter to predict CO₂ beyond 800 kyr. Mechanistically, multiple processes are expected to contribute to the shared variance. A first-order factor is the dependency of CO₂ solubility on ocean temperature (e.g. Millero, 1995). From the simple solubility perspective, colder climate states with increased ice volume and colder ocean temperatures will drive increased ocean uptake of CO₂ (Berends et al., 2021b). However, the solubility effect only accounts for a portion of observed glacial CO₂ drawdown (Archer et al., 2000). Multiple additional contributors to the shared variance are proposed in the literature. These include (not exhaustively) direct radiative forcing of ice volume changes by CO₂ (e.g. Shackleton and Pisias, 1985), the impact of ice volume and/or sea level changes on atmospheric CO₂ via ocean productivity and carbonate chemistry changes (e.g. Broecker, 1982; Archer et al., 2000; Ushie and Matsumoto, 2012), CO₂ drawdown during periods of high ice volume by increased iron fertilisation (e.g. Röthlisberger et al., 2004; Martinez-Garcia et al., 2014), and enhanced sea ice extent during periods of high ice volume capping the ventilation of CO₂ from the ocean interior at high latitudes (Stephens and Keeling, 2000).

Figure 2 Scatterplot of the composite observed atmospheric CO₂ record (Bereiter et al., 2015) against the LR04 benthic stack of marine δ¹⁸O records (Lisiecki and Raymo, 2005). The red line is a linear line of best fit (R²=0.68; p<0.05).

Quantitative separation and attribution of the processes linking global ice volume, ocean temperature, and atmospheric CO₂ on millennial to orbital timescales are not currently available (e.g. Archer et al., 2000; Sigman et al., 2010; Gottschalk et al., 2019) and will not be attempted here. Rather, we make the simple assumption that the relationships between the LR04 benthic δ¹⁸O stack and CO₂ can be extended beyond 800 ka and use generalised least-squares (GLS) regression modelling between benthic δ¹⁸O and CO₂ to make a prediction of CO₂ spanning 800–1500 ka. The deliberately simple implicit assumption, and null hypothesis, is that there is no change to the feedback processes linking benthic δ¹⁸O and CO₂ before and after the MPT.

This approach differs from previous more complex model studies that have attempted to reconstruct CO₂ using the LR04 benthic δ¹⁸O stack as an input variable (van de Wal et al., 2011; Stap et al., 2016; Berends et al., 2021a). Previous studies have used an inverse forward modelling approach, in which climate and ice sheet models of various complexities are used to capture physical relations between CO₂, global temperature, and ice volume. For example, in Berends et al. (2021a) the offset between modelled and observed benthic δ¹⁸O is used to calculate a value for atmospheric CO₂ that is iterated back to the inverse model. The CO₂ record which minimises the difference between the modelled and observed benthic stack is then taken as an estimate of how atmospheric CO₂ may have evolved to force coupled climate, deep-ocean temperature, and land ice volume changes that reproduce the observed benthic δ¹⁸O signal. The accuracy of the reconstructions in the inverse modelling approach depends on the ability of the climate and ice sheet models used to capture the correct climate dynamics across the MPT. Our GLS method is a simpler statistical approach, designed with the specific null hypothesis in mind, that does not attempt to simulate the physics linking benthic δ¹⁸O signal, land ice volume, global temperature, and CO₂. A range of approaches to reconstructing CO₂ have been called for and are of value in the context of forthcoming continuous ice core records across the MPT from oldest ice projects currently underway in Antarctica (IPICS, 2020).

To test our null hypothesis in advance of the recovery of a continuous ice core, we compare our predicted CO₂ record to two sets of low-resolution ice core data that exist outside the current 800 kyr observed CO₂. These data come from direct CO₂ measurements from ancient “blue ice” from the Allan Hills in East Antarctica (hereafter referred to as BI-CO₂) from ca. 1 Ma (Higgins et al., 2015) and 1.5 Ma (Yan et al., 2022). We use the term blue ice to describe deep, ancient glacial ice that has been brought nearer to the surface of an ice sheet by ice flow. Blue ice is sampled by cutting trenches or shallow drilling of up to several hundred metres (e.g. Higgins et al., 2015). The vertical migration of blue ice is associated with high deformation, making the ice samples stratigraphically complex and hard to date (Higgins et al., 2015). As a result, blue ice records alone do not provide a continuous CO₂ record across the MPT. In the Discussion section, we also compare our predicted record to existing proxy CO₂ reconstructions from boron isotope analysis of benthic foraminifera in marine sediment records (Chalk, et al., 2017; Dyez et al., 2018; Guillermic et al., 2022), leaf wax δ¹³C carbon isotope ratios (Yamamoto et al., 2022), and predictions from previous models of various complexities (van de Wal et al., 2011; Willeit et al., 2019; Berends et al., 2021a). We conclude with a discussion of the implications of our results and data comparisons for the understanding MPT dynamics.

2 Methods

We use a generalised least-squares (GLS) model with an autoregressive (AR) factor 1 to predict atmospheric CO₂ from the LR04 benthic δ¹⁸O stack (Fig. 3a and b). We use GLS because the assumptions of ordinary least squares (OLS) are violated by the presence of autocorrelation and heteroskedasticity in the regression errors. We selected the AR(1) correlation factor as it yielded the lowest Akaike information criterion (AIC) value from a test of multiple correlation factors. The AR(1) process assumes and accounts for the dependence of error at a given point in time on the previous error term. In practise this makes the model assumptions more realistic and improves parameter estimation where, as in the climate system, observations are dependent on past values.

Figure 3 (a) The LR04 benthic stack of 57 globally distributed δ¹⁸O records (Lisiecki and Raymo, 2005). (b) Comparison of our PRED-CO₂ (ppm) record to the current continuous composite record (0–800 ka) and to direct CO₂ measurements from Allan Hills blue ice cores (BI-CO₂) at ca. 1 Ma (±89 kyr) (Higgins et al., 2015) and ca. 1.5 Ma (±213 kyr) (Yan et al., 2022). Age uncertainty boundaries for the BI-CO₂ data are represented by dashed box boundaries. Marine isotope stages 5, 9, and 11 are numbered on the plot according to Lisiecki and Raymo (2005). Blue shading around PRED-CO₂ is the 95 % CI from bootstrap analysis. (c) Mean concentrations of the PRED-CO₂ and observed composite CO₂ records over the range of the observed composite record (offset for clarity) and the mean concentrations of the PRED-CO₂ and BI-CO₂ data at 1 Ma and again at 1.5 Ma averaged over the age uncertainty range of each BI-CO₂ dataset. (d) As for (c) but filtered by the upper and lower 25th and 75th percentiles to estimate glacial and interglacial periods.

To obtain common time steps and resolution between the predictor (LR04 benthic δ¹⁸O stack) and response (CO₂) variables, we re-grid the LR04 benthic stack and Bereiter et al. (2015) CO₂ data into time bins with a resolution of 3 kyr. The GLS regression model was then applied over the 0–800 kyr range of the predictor and response variables as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathrm{CO}}_{\mathrm{2}}& {\displaystyle }=-\mathrm{33.37}\times \mathit{\delta }{}^{\mathrm{18}}\mathrm{O}+\mathrm{365.15},\\ {\displaystyle }& {\displaystyle }\text{autoregressive (AR) factor:}\mathrm{1}.\end{array}\end{array}$$

Based on the regression model, the δ¹⁸O values of the LR04 benthic stack from 800–1500 ka were used to predict CO₂ concentration over this range (hereafter referred to as PRED-CO₂). To gauge the GLS model stability we took a bootstrap approach, selecting a random 50 % subset of our data (with replacement) and re-running the model 1000 times to determine 95 % confidence intervals (CIs) for the predictions. While the GLS method itself addresses autocorrelation, the bootstrap method introduces variability such that each iteration of the model has different combinations of the original data points (including repeated ones); this variability helps assess the robustness and sensitivity of the model e.g. to variable data and dating uncertainty.

Uncertainties in the independent age scales of both the LR04 stack and the compiled CO₂ record are inherited by our GLS model and its predictions. The LR04 stack includes 57 globally distributed benthic δ¹⁸O sediment core records. The age models for these cores are constructed by alignment of their δ¹⁸O signals, followed by tuning of the stack to a simple ice model based on 21 June insolation at 65° N in a way which maintains relatively stable global mean sedimentation rates (Lisiecki and Raymo, 2005). The authors estimate uncertainty of 6 kyr from 1.5–1.0 Ma and 4 kyr from 1–0 Ma (Lisiecki and Raymo, 2005). The observed CO₂ composite ice core record for the past 800 kyr (Bereiter et al., 2015) uses six independent dating methods for various core locations both spatially across Antarctica and stratigraphically for different sections of the same core. The age uncertainty in the gas timescale has a median over the 0–800 ka interval of 2 kyr, but individual uncertainties can reach up to 5 kyr (Veres et al., 2013; Bazin et al., 2013). The relative age uncertainties between these input variables may diminish the regression or in some instances lead to spurious correlation. However, we expect any such effects are minor on the basis that our predictions show little sensitivity to the bootstrap analysis, with a median 2σ error of 5.8 ppm from 0 to 1.8 Ma (see Fig. 3b and c and the Discussion section).

3 Results

Figure 3b shows the time series of our LR04-benthic-δ¹⁸O-stack-based GLS model predictions of atmospheric CO₂ (PRED-CO₂) over the past 800 kyr in comparison to the observed ice core CO₂ record from Bereiter et al. (2015). The correlation coefficient (R²) between the predicted and observed records is 0.68 (p≪0.01). Our PRED-CO₂ record out to 1.8 Ma with shaded 95 % CIs from the bootstrap analysis is also shown, overlain with observed Allan Hills blue ice CO₂ (BI-CO₂) datasets of age 1000 ± 89 ka (Higgins et al., 2015) and 1.5 Ma ± 213 kyr (Yan et al., 2022).

We evaluate the PRED-CO₂ record against the observed CO₂ data according to criteria of mean concentrations across the common intervals and mean concentrations in the glacial and interglacial subsets of the data. First is the mean CO₂ concentration over the common intervals (Fig. 3c). From 0–800 ka the mean concentrations in observed (Bereiter et al., 2015) and PRED-CO₂ data are in close agreement (225.2 ± 3.03 ppm versus the predicted 225.2 ± 2.5 ppm, respectively; uncertainties are 95 % confidence intervals, i.e. 1.96σ). In the 1000 ± 89 ka interval (i.e. averaged across the age uncertainty of the Higgins et al., 2015, blue ice data) the BI-CO₂ concentration is ∼11 ppm higher than PRED-CO₂ (246.7 ± 8.4 ppm versus the predicted 235.3 ± 3.9 ppm); this difference is not significant at the 95 % confidence level. For the 1.5 Ma ± 213 kyr interval, the mean BI-CO₂ concentration is ∼9 ppm lower than PRED-CO₂ (231.9 ± 5.6 ppm versus the predicted 240.7 ± 2.1 ppm), which is marginally significant at the 95 % level. Comparisons of mean levels across intervals spanning multiple glacial and interglacial cycles may be biased if (as is likely) the blue ice data are not sampling glacial and interglacial values with the same uniformity as a continuous record.

To address this, we define the glacial and interglacial thresholds of PRED-CO₂ to be the lower and upper 25th percentiles, respectively, of the LR04 δ¹⁸O predictor variable (following Chalk et al., 2017). Filtering the observed (Bereiter et al., 2015) CO₂ record and our predicted CO₂ record according to these definitions, we find a very close match for glacial (202.0 ± 3.2 versus the predicted 199.7 ± 1.6 ppm) and interglacial intervals (253.9 ± 4.1 ppm versus the predicted 253.1 ± 2.3 ppm) over the past 800 kyr (see Fig. 3d for these comparisons). For blue ice (BI-CO₂) data, a corresponding LR04 isotope signal could not be confidently applied to the measured CO₂ concentration due to the uncertainties associated with blue ice dating; therefore, we defined the glacial and interglacial thresholds of blue ice data according to the top (interglacial) and bottom (glacial) 25th percentiles of actual CO₂. Applying this to the 1000 ± 89 ka interval finds that observed BI-CO₂ data are ∼9 ppm higher than PRED-CO₂ during the glacial stages (226.2 ± 4.0 ppm versus the predicted 217.6 ± 2.3 ppm) and ∼15 ppm higher than PRED-CO₂ during the interglacial stages (271.3 ± 4.5 versus the predicted 256.3 ± 3.8 ppm). These differences are significant with respect to the constrained uncertainties. During the 1.5 Ma ± 213 kyr interval, the mean BI-CO₂ concentration did not show any significant difference from PRED-CO₂ in interglacial stages (254.1 ± 10.3 versus the predicted 257.2 ± 1.7 ppm). During glacial stages there is a small 2.9 ppm difference between the upper estimate of BI-CO₂ and the lower estimate of PRED-CO₂ (218.4 ± 1.3 and 224 ± 1.4 ppm, respectively; see Fig. 3d). In our view these results, notwithstanding the 2.9 ppm difference at 1.5 Ma, do not give sufficient cause to reject the GLS model. Furthermore, the comparison indicates that PRED-CO₂ is not drifting systematically away from the existing observed BI-CO₂ data (Fig. 3d). The differences could of course represent a failing in the model, potential biases in the blue ice data, dating uncertainty, and/or other unconstrained uncertainties (see the Discussion section for blue ice caveats).

Figure 4 (a) Predicted CO₂ (this work) compared to observed proxy CO₂ estimates from a range of other sources: δ¹¹B-based pCO₂ reconstructions and measurements by Dyez et al. (2018), Guillermic et al. (2022), and Chalk et al. (2017); blue ice CO₂ measurements by Yan et al. (2019) and Higgins et al. (2015); and δ¹³C leaf wax proxy reconstructions by Yamamoto et al. (2022). The dashed boxes indicate the dating uncertainty and range of the respective BI-CO₂ records. (b) Our predicted record compared to various model simulations: a regolith removal hypothesis simulation by Willeit et al. (2019) and inverse-model-based CO₂ reconstructions by van de Wal et al. (2011) and Berends et al. (2021a).

We now consider long-term trends in interglacial and (separately) glacial CO₂ levels across the past 1.8 Myr in PRED-CO₂ and in the existing ice core CO₂ data. For PRED-CO₂ there is no significant difference between CO₂ concentrations in the interglacial stages of the 1.5 Ma ± 213 ka, 1000 ± 89 ka, and 0–800 ka windows (Fig. 4d, blue bars). In the ice core observations, interglacial levels at 1.5 Ma in BI-CO₂ are also within the uncertainties of those in the 0–800 ka interval. Notably, the BI-CO₂ concentrations in the 1000 ± 89 ka interval appear elevated with respect to the 0–800 kyr and 1.5 Ma ± 213 ka intervals; however, this elevated (ca. 271 ppm) level is consistent with the observed interglacial CO₂ concentration during interglacials 5, 9, and 11 (Fig. 3b). Overall, there is no indication in the observed ice core CO₂ data or in PRED-CO₂ of a long-term trend in interglacial CO₂ levels across the past 1.8 Myr.

In comparison, there are significant declines in glacial CO₂ levels across the MPT in PRED-CO₂ and the observed ice core data. For PRED-CO₂, glacial CO₂ concentrations are not significantly different during the 1.5 Ma ± 213 ka and 1000 ± 89 ka windows. However, across the MPT, PRED-CO₂ glacial concentrations drop by ∼18 ppm (Fig. 3d). This pattern is similar to the observed BI-CO₂ data, where glacial CO₂ levels show no decline between the 1.5 Ma ± 213 ka and 1000 ± 89 ka windows (indeed there is a marginal increase from 218.4 ± 1.3 to 226.2 ± 4.0 ppm, respectively), before falling by 24 ppm to the 0–800 kyr observed glacial mean of 202.0 ± 3.2 ppm (Fig. 3d). Glacial-stage drawdown of CO₂ across the MPT in the absence of interglacial drawdown is consistent with previous observations of boron-isotope-based CO₂ reconstructions (e.g. Chalk et al., 2017; Hönisch et al., 2009; see the Discussion section). In the following section we also compare PRED-CO₂ data to boron-isotope-based and other CO₂ proxy records covering the 0 to 1.8 Myr interval.

4 Discussion

Our objective with this paper was to generate the simplest reasonable model to predict CO₂ from the LR04 δ¹⁸O benthic stack and to test the predictions against available observations. It is possible that the fit between observed and our predicted CO₂ data could be further improved using a non-linear approach. However, we refrain from a non-linear approach for several key reasons. First, a scatterplot of the LR04 δ¹⁸O benthic stack versus observed ice core CO₂ over the past 800 kyr yields a Pearson's correlation coefficient (R) of −0.82 (Fig. 2), indicating that ∼68 % of the variance in observed CO₂ is shared with the benthic stack. This is similar to what was reported in ordinary linear least-squares regression (R²=0.70) by Berends et al. (2021a). Importantly, there is no evidence in this scatterplot of departure from the linear relationship at high or low CO₂ or benthic δ¹⁸O levels. Second, following the approach of Chalk et al. (2017) and interpreting the upper 25th percentile of CO₂ data as representing mean interglacial-stage CO₂ and the lower 25th percentile of CO₂ data as representing mean glacial-stage CO₂ levels, we see that our predicted interglacial mean value for the past 800 kyr (253.1 ± 2.3 ppm) closely overlaps with the observed interglacial mean value (253.9 ± 4.1 ppm), and similarly, the predicted glacial-stage mean (199.7 ± 1.7 ppm) closely overlaps with the observed glacial-stage mean (202.0 ± 3.2 ppm). Third, the predictions are remarkably insensitive to bootstrap analysis in which 50 % of the data are omitted with each iteration of the GLS model. Such insensitivity to the bootstrap analysis and accurate prediction of glacial and interglacial CO₂ values would be unlikely in the case of major non-linear dependencies between the LR04 predictor and CO₂ response variables. Fourth, non-linear approaches would risk generating an improved fit due to statistical artefacts that do not meaningfully relate to any dependence between benthic δ¹⁸O and CO₂. Finally, the specific causes and sources and sinks involved in glacial to interglacial and millennial-scale CO₂ variations remain poorly constrained (e.g. Archer et al., 2000; Sigman et al., 2010; Gottschalk et al., 2019). Given this process uncertainty, the GLS model fits our criteria of the simplest reasonable model. Further, the use of benthic δ¹⁸O to predict atmospheric CO₂ has precedence; in response to the EPICA challenge (Wolff et al., 2004) Nicholas  Shackleton predicted atmospheric CO₂ out to 800 kyr based on a number of benthic δ¹⁸O records from the eastern Pacific (Wolff et al., 2005).

There are several caveats with blue ice data that may affect their use in evaluating our GLS model predictions. The blue ice data may have been subject to diffusional smoothing of CO₂ (e.g. Yan et al., 2019), which would act in the direction of elevating the (lower 25th percentile) assumed glacial concentrations above the glacial atmospheric values and reducing the (upper 25th percentile) assumed interglacial concentrations. There is also the potential for artificially elevated CO₂ concentrations in blue ice due to in situ respiration of CO₂ from microbial activity in detrital matter. Respiration effects are screened for by measurements of δ¹³C of CO₂, but it is difficult to demonstrate that all samples are unaffected (Yan et al., 2019). These uncertainties support our argument that the GLS model predictions are not rejected by the available observed BI-CO₂ data.

We consider the BI-CO₂ data to provide the most reliable measurements of CO₂ concentration in the absence of a continuous ice core record across the MPT. However, further comparison of our CO₂ predictions can also be made against CO₂ proxy data from non-ice-core archives (Fig. 4a). We consider δ¹¹B-based atmospheric CO₂ reconstructions (Chalk et al., 2017; Dyez et al., 2018; Guillermic et al., 2022) and a recent atmospheric CO₂ reconstruction from δ¹³C of leaf wax (Yamamoto et al., 2022). The continuous δ¹¹B-based reconstructions of Dyez et al. (2018) overlap PRED-CO₂ from ∼1.38–1.5 Ma, while the Chalk et al. (2017) reconstruction overlaps PRED-CO₂ from 1.09–1.43 Ma. Discrete reconstructions from Guillermic et al. (2022) are distributed non-uniformly across the ∼800 to 1.5 Ma interval. For the two continuous δ¹¹B-based reconstructions (Chalk et al., 2017; Dyez et al., 2018) the glacial CO₂ levels appear to be consistent with the PRED-CO₂ record within their reported 30–60 ppm uncertainties. However, δ¹¹B-based interglacial stages in these reconstructions exceed those of the PRED-CO₂ record (Fig. 4a). The Guillermic et al. (2022) reconstructions suggest a larger range of CO₂ concentrations than the overlapping intervals of PRED-CO₂ and of the two continuous δ¹¹B-based reconstructions (Fig. 4a). The large range of the Guillermic et al. (2022) data and the high interglacial maxima in the Chalk et al. (2017) and Dyez et al. (2018) data all significantly exceed the range and interglacial maxima from the BI-CO₂ estimates. These discrepancies internally between different δ¹¹B-based CO₂ reconstructions and between the δ¹¹B-based reconstructions and the BI-CO₂ data may be due to uncertainties associated with the δ¹¹B proxy transfer function. The δ¹¹B-based CO₂ reconstructions are dependent on assumptions about multiple components of the carbonate system, including local marine carbon chemistry and the CO₂ saturation state in the past (Hönisch et al., 2009). Evidence that δ¹¹B-based reconstructions may overestimate interglacial-stage CO₂ is also seen in data from Chalk et al. (2017) spanning ca. 0–250 ka, where the δ¹¹B-based interglacial CO₂ levels exceed the continuous ice core CO₂ record by up to ca. 30 ppm.

By comparison, the δ¹³C of leaf wax data (Yamamoto et al., 2022) have a similar glacial to interglacial range as PRED-CO₂ but a ca. 20 ppm lower mean concentration than our predictions (Fig. 4a). Hence, our PRED-CO₂ data fall lower than interglacial δ¹¹B-based interglacial levels but are higher than the δ¹³C of leaf-wax-based estimates. The strong spread between these different proxies and the large associated uncertainty of the alternative marine and leaf wax proxy CO₂ reconstructions mean that we do not find cause from the existing CO₂ proxy data to reject our predictions or our associated null hypothesis.

We also compare our predictions to existing more complex model simulations (Fig. 4b), first against a transient simulation using an intermediate-complexity Earth system model (CLIMBER-2) by Willeit et al. (2019). This study suggests that a combination of gradual regolith removal and atmospheric CO₂ decline can explain the long-term climate variability over the past 3 Myr. The second comparison is against a longer-term reconstruction by van de Wal et al. (2011), which uses benthic δ¹⁸O that utilises deep-sea benthic isotope records to reconstruct a continuous CO₂ record over the past 20 Myr. The third comparison is to a CO₂ reconstruction based on an inverse forward modelling approach forced by the LR04 benthic stack, in which the forward model is incrementally updated through interaction with general circulation model snapshots and the ANICE 3-D ice sheet–shelf model (Berends et al., 2021a). Our simple GLS model demonstrates a similar long-term trend and timing of glacial–interglacial signals as well as an atmospheric CO₂ level that sits approximately mid-way between the van de Wal et al. (2011) and Willeit et al. (2019) models and is remarkably similar to the Berends et al. (2021a) reconstruction, despite their different approach. Notably the Berends et al. reconstruction shows greater glacial to interglacial amplitude in the CO₂ signal compared to our GLS model. The decreasing linear trend in CO₂ in Willeit et al. (2019), which is not seen in the other reconstructions, was directly prescribed in that study to induce Northern Hemisphere glaciation at 2.6 Myr ago.

A complete and critical test of our and other CO₂ predictions awaits the upcoming analysis of the continuous oldest ice core records. We now discuss some potential applications of the PRED-CO₂ record for hypothesis testing on the cause of the MPT.

PRED-CO₂ shows a long-term decline in glacial CO₂ across the MPT but no long-term decrease in interglacial CO₂. This pattern is consistent with the boron-isotope-based CO₂ reconstructions shown earlier, where it is often described as an increase in the interglacial to glacial CO₂ difference (e.g. Chalk et al., 2017; Hönisch et al., 2009). Chalk et al. (2017) conclude that the MPT was initiated by a change in ice sheet dynamics and that longer and higher-ice-volume post-MPT ice ages are sustained by carbon cycle feedbacks, in particular dust fertilisation of the Southern Ocean. The fact that our LR04-based prediction of CO₂ captures this same trend, with predicted glacial CO₂ fairly constant from 1.5 to ca. 1.0 Ma before declining from 1.0 to 0.6 Ma, reflects the fact that the LR04 benthic stack also features an increase in the interglacial to glacial benthic δ¹⁸O difference across this same interval, which is dominated by the glacial-stage changes (Fig. 3a). Here, a comparison of PRED-CO₂ to a realised continuous oldest ice core record will be of value. The agreement or disagreement would provide information on the proportionality of the CO₂ coupling with ice volume; if there were a major new or non-linear process across the MPT that changed the nature of coupling between CO₂ and ice volume, the PRED-CO₂ and observed CO₂ records would be expected to diverge.

Another avenue to use the PRED-CO₂ record for hypothesis testing on the cause of the MPT concerns the phase locking hypothesis. The phase locking hypothesis is proposed to explain the absence of precession-related (23 kyr) periods in the LR04 benthic stack prior to the MPT (Fig. 1), despite the strong precession cycle in insolation (Raymo et al., 2006; Morée et al., 2021). The key concept is that prior to the MPT the Northern Hemisphere and Antarctic ice sheets were responsive (in ice volume) to insolation changes in the precession band, but because precession forcing is out of phase between the hemispheres, the ice volume changes were opposing between the hemispheres and therefore cancelled in the benthic stack. This cancellation of the precession signal left insolation forcing in the 41 kyr obliquity band to dominate globally integrated ice volume changes expressed in the benthic stack. A transition from a smaller and more dynamic terrestrial-terminating Antarctic ice sheet to a larger and more stable marine-terminating ice sheet with cooling climate across the MPT (e.g. Elderfield et al., 2012) is then proposed to remove sensitivity of Antarctic ice volume to local precession forcing in favour of quasi-100 kyr ice volume changes that are in phase between the hemispheres (Raymo et al., 2006).

Recently presented data from Yan et al. (2022) lend some support to the phase locking hypothesis, specifically with evidence that pre-MPT Antarctic temperature (and by extension ice volume) is positively correlated with a local precession band insolation proxy based on the oxygen to nitrogen ratio of trapped air (Yan et al., 2022), whereas the correlation becomes negative in the blue ice and continuous ice core data in the post-MPT record. If Yan et al. (2022) are correct and the phase locking hypothesis holds, then an implication is that prior to the MPT, Antarctic climate, Antarctic ice volume, and by extension Southern Ocean climate conditions would fall out of phase with the LR04 benthic stack. To now extend the argument to potential impacts on CO₂ exchange, if the phase locking hypothesis holds, then prior to the MPT the Antarctic and Southern Ocean climate conditions and by extension the Southern Ocean mechanisms of CO₂ exchange described earlier would also be expected to fall out of phase with the benthic stack. Since our regression model assumes continuation of the in-phase relationship between the benthic stack and Antarctic and Southern Ocean climate conditions (as inherited from the post-MPT training data) we would expect to see major disagreement between our pre-MPT CO₂ predictions and a realised oldest ice continuous ice core CO₂ record.

5 Summary and conclusions

In this study we have used a simple generalised least-squares (GLS) model to predict atmospheric CO₂ from the LR04 benthic δ¹⁸O stack for the period spanning the Mid-Pleistocene transition: 800–1800 kyr. Our CO₂ prediction is therefore based on the assumption that the physical processes linking CO₂, sea level, global ice volume, and ocean temperature over the past 800 kyr do not fundamentally change across the 800–1800 ka time period. The null hypothesis is deliberately simplistic on the basis that differences between our predictions and observed or proxy CO₂ records may be revealing of the physical processes involved in the Mid-Pleistocene Transition.

We made initial tests of the null hypothesis by comparing our predicted CO₂ record to existing discrete blue ice CO₂ records and other non-ice-core proxy CO₂ records from the 800–1800 kyr interval. Our predicted CO₂ concentrations do not show any systematic departure from observed blue ice CO₂ concentrations. The predictions are marginally lower (during glacial and interglacial stages) than those observed in blue ice from 1000 ± 89 ka and marginally higher than observed in blue ice data from 1.5 Ma ± 213 kyr. Our predictions were generally lower than interglacial δ¹¹B-based CO₂ reconstructions but higher than recent δ¹³C of leaf-wax-based CO₂ reconstructions. Overall, we do not find clear evidence from the existing blue ice or proxy CO₂ data to reject our predictions or our associated null hypothesis. The definitive test of our and other CO₂ predictions therefore awaits the future analysis of the upcoming continuous oldest ice core records. The PRED-CO₂ record presented here should provide a useful comparison to forthcoming oldest ice core records and an opportunity for further constraints on the processes involved in the MPT.