Causes of ice age intensification across the Mid-Pleistocene Transition.

During the Mid-Pleistocene Transition (MPT; 1,200-800 kya), Earth's orbitally paced ice age cycles intensified, lengthened from ∼40,000 (∼40 ky) to ∼100 ky, and became distinctly asymmetrical. Testing hypotheses that implicate changing atmospheric CO2 levels as a driver of the MPT has proven difficult with available observations. Here, we use orbitally resolved, boron isotope CO2 data to show that the glacial to interglacial CO2 difference increased from ∼43 to ∼75 μatm across the MPT, mainly because of lower glacial CO2 levels. Through carbon cycle modeling, we attribute this decline primarily to the initiation of substantive dust-borne iron fertilization of the Southern Ocean during peak glacial stages. We also observe a twofold steepening of the relationship between sea level and CO2-related climate forcing that is suggestive of a change in the dynamics that govern ice sheet stability, such as that expected from the removal of subglacial regolith or interhemispheric ice sheet phase-locking. We argue that neither ice sheet dynamics nor CO2 change in isolation can explain the MPT. Instead, we infer that the MPT was initiated by a change in ice sheet dynamics and that longer and deeper post-MPT ice ages were sustained by carbon cycle feedbacks related to dust fertilization of the Southern Ocean as a consequence of larger ice sheets.

The Mid-Pleistocene Transition (MPT) marks a major shift in the response of Earth’s climate system to orbital forcing. During the Early Pleistocene, glacial–interglacial (G-IG) climate cycles were paced by ∼40,000 y (40 ky) obliquity cycles, whereas G-IG cycles after the MPT gradually intensified over multiple obliquity cycles (i.e., 80- to 120-ky periodicity) (1, 2) and acquired a distinctively asymmetric character with gradual glacial growth and abrupt glacial terminations that were paced by a combination of obliquity and precession (1). These changes gave rise to longer, colder, and dustier Late Pleistocene ice ages with larger continental ice sheets and lower global sea level (SL) (3–5) (Fig. 1). The MPT occurred in the absence of any significant change in the pacing or amplitude of orbital forcing, indicating that it arose from an internal change in the response of the climate system rather than a change in external forcing (1, 6, 7).

Fig. 1.

Climate records across the MPT. (A) CO2 records are shown as follows: black line, ice core compilation (14); blue, our δ11B-based LP260 data; red, our δ11B-based eMPT data; and purple squares, low-resolution MPT δ11B record of ref. 11 (all with 2σ error bars/envelopes). The range of ice core CO2 measurements (17) from stratigraphically disturbed blue ice and their approximate ages are indicated. (B) SL records, where orange indicates the Red Sea record (21), dark blue represents Mg/Ca-based deconvolution of deep sea benthic foraminiferal oxygen isotope data (3), and pink shows a record from the Mediterranean Sea (4). (C) Dust mass accumulation rate (MAR) in a sub-Antarctic site ODP 1090 on the southern flank of the Agulhas Ridge (24). (D) LR04 benthic foraminiferal oxygen isotope stack (26). Warm intervals are highlighted by gray bars.

Proposed explanations for the MPT fall into two primary groups: those that invoke a change in ice sheet dynamics and those that call on some subtle change in the climate system’s global energy budget. Two prominent hypotheses posit that either removal of the subglacial regolith beginning at about 1,200 ky (8, 9) or phase-locking of Northern and Southern Hemisphere ice sheets at about 1,000 ky (10) gave rise to deeper and ultimately longer G-IG climate cycles by allowing for a greater buildup of ice independent of a change in CO2 radiative climate forcing (scenario 1 in Fig. 2). Alternatively, it has been argued that an underlying change in the global carbon cycle could have triggered the MPT through a decline in ΔRCO2 [i.e., the radiative climate forcing exerted by CO2 decline (11–13) (scenario 2 in Fig. 2)]. The continuous 800-ky-long ice core record of atmospheric CO2 (i.e., compiled by ref. 14) is well-correlated to and shares spectral power with orbital-scale changes in temperature, ice volume, SL, and the oxygen isotopic composition of benthic foraminifera (Figs. 1 and 3). State of the art coupled climate–ice sheet models can simulate climate cycles that are longer than single obliquity cycles, provided that mean CO2 concentrations are within certain model-dependent bounds (15, 16) (e.g., 200–260 μatm). These studies suggest that the absolute CO2 level attained during rising obliquity (i.e., during increasing high-latitude Northern Hemisphere summer insolation) may be a critical control that determines whether ice sheets are strictly locked to the ∼40-ky beat of obliquity or survive for longer periods. Recent work has provided some evidence for an overall CO2 decline since the MPT (11, 17), supporting this view. The study by Hönisch et al. (11), in particular, provides evidence that CO2 decline was most pronounced during glacial stages. Here, we build on that work with the aim to resolve the coupling of CO2 and climate on orbital timescales to address major unanswered questions regarding the role of CO2 change in the MPT.

Fig. 2.

Changing relationship between CO2 climate forcing and ice sheet size. Three scenarios (A–C) for the MPT intensification of glacial cycles compared with observations (D). Reconstructed SL is taken here to reflect continental ice sheet size in relationship to CO2 climate forcing (∆RCO2) calculated (33) from our orbitally resolved CO2 data. In all panels, red and blue represent conditions during our two sampling intervals before and after the MPT (i.e., eMPT and LP260), respectively. The end member scenarios posit (A) a change in ice sheet dynamics, causing ice volume to become more sensitive to unchanged G-IG climate forcing, and (B) an unchanged sensitivity of ice sheet size to forcing, with glacial intensification driven by additional CO2 drawdown. Neither one of these two scenarios adequately describes both observed changes of increased ice sheet sensitivity (greater slope) and additional glacial CO2 drawdown (more negative climate forcing). Here, we argue for a hybrid scenario with a change in ice sheet dynamics (possibly caused by regolith removal of ref. 8 or ice sheet phase-locking of ref. 10), allowing ice sheets to grow larger and to trigger a positive ice–dust–CO2 feedback that promotes additional glacial intensification. In D, the regression confidence intervals account for uncertainty in both SL and ΔRCO2 (), but to avoid clutter, we only display the regression based on the Mediterranean SL reconstruction (4) and the uncertainty on the slope rather than the individual data points. We refer the reader to and Fig. S7 for other SL records and full treatment of data uncertainties.

Fig. 3.

Reconstructed ice age CO2 cycles before and after MPT. (A) Boron isotope data from ODP 999 (Fig. S1) shown in blue (LP260) and red (eMPT) along with the LR04 deep sea benthic foraminiferal oxygen isotope stack (black) (26). (B) CO2 levels calculated from boron isotopes (same colors as above) compared with ice core (black) (14) and previous low-resolution boron-derived CO2 data (purple) (11). Probabilistic assessments are shown as the colored bands, with the probability maximum shown within a dark band that represents its 95% probability envelope (∼±6 ppm) and a lighter band that represents the full 95% envelope of the sampled distribution. As illustrated by B, Inset, comparison between our (red) eMPT and (blue) LP260 records reveals that glacials on average experienced higher CO2 levels during eMPT than LP260 (eMPT: 241 ± 21 μatm vs. LP260: 203 ± 14 μatm; 2σ), whereas interglacial levels were indistinguishable between the two time slices (eMPT: 284 ± 17 μatm vs. LP260: 277 ± 18 μatm; 2σ).

To better quantify the role of CO2 during the MPT, we present two orbitally resolved, boron isotope-based CO2 records generated using the calcite tests of surface-dwelling planktonic foraminifera from Ocean Drilling Program (ODP) Site 999 in the Caribbean (Fig. 3 and Figs. S1 and S2). Boron isotopes (δ11B) in foraminifera have proven to be a reliable indicator of past ocean pH (18, 19) and with appropriate assumptions regarding a second carbonate system parameter ( and Fig. S3), allow reconstruction of atmospheric CO2 levels. Site 999 likely remained near air–sea CO2 equilibrium through time (20), and this is further supported by agreement of our data (blue and red in Figs. 1 and 3) with published low-resolution δ11B-derived CO2 data from ODP Site 668 in the equatorial Atlantic (11) (purple squares in Figs. 1 and 3) and with the ice core CO2 compilation (14).

Results

Our two datasets span an early portion of the Mid-Pleistocene Transition (eMPT) from 1,080 to 1,250 kya (n = 51) and for validation against the ice core CO2 record, the Pleistocene interval from 0 to 260 kya (LP260; n = 59, including 32 recalculated data points from ref. 18), yielding a similar median sampling interval of ∼3.5–4.5 ky for both records. Our LP260 CO2 dataset has a confidence interval of ±20 μatm (2σ) and is offset by a mean of +7 μatm from the ice core CO2 data when accounting for both CO2 and age uncertainties (21) (Fig. 3 and ). Comparison between our two CO2 records reveals that eMPT glacials on average were associated with higher CO2 levels than LP260 glacials (eMPT: 241 ± 21 μatm vs. LP260: 203 ± 14 μatm; 2σ), whereas interglacial levels were indistinguishable (eMPT: 284 ± 17 μatm vs. LP260: 277 ± 18 μatm; 2σ). This analysis uses highest and lowest 25th percentiles of δ18O values to define “glacial” and “interglacial” subsets of the data, although this pattern is independent of the thresholds that we define (Fig. 4 and Fig. S4). Our analysis reproduces the glacial stage-specific decline in CO2 levels found in ref. 11, leading to similar reconstructed increases in the glacial to interglacial CO2 difference since the MPT (40 ± 47 and 32 ± 35 μatm based on ref. 11 and our data, respectively) (Fig. 4). The higher resolution of these datasets allows this approach to yield useful data about our timespans, despite the relatively large uncertainty on each individual data point. When analyzed in a similar way, recent direct measurements of CO2 from a stratigraphically disturbed section of ∼1-My-old “blue ice” (17) offer a fully independent test for the two δ11B-based reconstructions and are consistent with these findings (Fig. 4 and Fig. S4). Thus, all available evidence suggests that the MPT was associated with a transition in the global carbon cycle characterized mainly by enhanced glacial-stage drawdown of CO2.

Fig. 4.

CO2 change since the MPT. Quantified from different datasets: boron isotope data from ODP 999 (this study) and ODP 668 (11), CO2 directly measured on stratigraphically disturbed ∼1-My-old blue ice from the Allan Hills (17), and CYCLOPS model inversion (this study). For each dataset, we quantify the change in (Top) interglacial and (Middle) glacial CO2 level as well as (Bottom) the change in the magnitude of interglacial–glacial CO2 cycles. For this analysis, we define glacial and interglacial subsets of the datasets based on a 25% cutoff criterion, subsampling the data with the 25% lowest/highest δ18O (marine records) or CO2 (ice core; model). As further discussed in , the results are robust for a wide range of cutoff values (Fig. S4). Thick black bars denote 1σ uncertainty of the estimated CO2 change, while thin black bars denote the one-sided test of the sign of CO2 change at 95% significance level. We note that the ODP 668 uncertainties do not encompass the underlying alkalinity and seawater boron isotope composition assumptions, which are included in the uncertainty propagation for our ODP 999 data. The Allan Hills ice may not capture the full range of CO2 levels (17).

We evaluate the reconstructed G-IG CO2 change across our study interval with a carbon cycle model inversion of Southern Ocean and Atlantic mechanisms thought to have contributed to the most recent Late Pleistocene G-IG CO2 cycles (22). For this, we force the CYCLOPS carbon cycle model (23) with ODP 1090 sedimentary iron mass accumulation rates (24), ODP 1094 Ba/Fe ratios (25), and ODP 982/U1313 (Fig. S1) benthic ∆δ13C variations (26, 27) to represent, respectively, (i) sub-Antarctic dust-borne iron fertilization; (ii) combined changes in polar Antarctic stratification, nutrient drawdown, and export production; and (iii) transitions in the geometry and depth structure of the Atlantic Meridional Overturning Circulation (AMOC) (Fig. S5). These mechanisms and their model sensitivities have been documented elsewhere (23). Here, we invert the model and the forcing to minimize the mismatch between simulated atmospheric CO2 levels and the ice core CO2 record of the last 800 ky (residual rms error of 12.3 μatm) () and then, to predict atmospheric CO2 levels back to 1,500 ky (Fig. S5) for comparison with our data.

We find that changes in the periodicity of simulated CO2 levels closely match those in the ice core CO2 record, in the benthic foraminiferal oxygen isotope record, and in our δ11B-based CO2 reconstruction (Fig. S6). Within the relative age uncertainty between the model forcing and our δ11B record, we find that the model explains more than 60% of the variance observed in our eMPT CO2 reconstruction, in line with model and reconstruction uncertainties. The model inversion does not include any secular change in the silicate weathering cycle (11) (), so that simulated CO2 change is exclusively related to carbon redistribution within the ocean–atmosphere system and associated CaCO3 compensation dynamics (22, 23).

In good agreement with the δ11B-based CO2 reconstructions and the ice core CO2 measurements, the model inversion yields (i) insignificant (−1 ± 3 μatm; 2σ) eMPT to LP260 interglacial CO2 change and (ii) a −22 ± 5 μatm (2σ) eMPT to LP260 decline in glacial-stage CO2 levels (Fig. 4 and Fig. S4). In the model, we can attribute most of the additional glacial CO2 drawdown to MPT intensification of glacial dust-borne iron fertilization of biological productivity and nutrient utilization in the Sub-Antarctic Zone of the Southern Ocean (24, 28–30) (Fig. S5). AMOC shoaling also seems to have become more prevalent after ∼1,200 ky but contributes less to simulated CO2 change (23). The model reproduces relatively low reconstructed interglacial CO2 from 400 to 800 ky, because use of ODP 1094 Ba/Fe in the model inversion results in persistent polar Southern Ocean stratification as suggested previously (25). Through our eMPT sample interval, the model reproduces the ∼80-ky CO2 periodicity that is evident in our eMPT δ11B data (Fig. S6), mainly because of an ∼80-ky periodicity in eMPT polar Antarctic stratification and nutrient cycling recorded in ODP 1094 Ba/Fe (25). While all three forcings (iron fertilization, Atlantic circulation, coupled polar Antarctic changes) contribute to the simulated changes in CO2 periodicities that are highly coherent with the MPT change in rhythm of the climate system, the iron fertilization influence dominates the MPT intensification of ice age CO2 drawdown (Fig. S5).

Discussion

MPT intensification of glacial-stage CO2 drawdown is consistent with stabilization of continental ice sheets during increasing orbital obliquity by reduced greenhouse gas forcing, thereby helping ice sheets to grow larger and for periods longer than one obliquity cycle (scenario 2 in Fig. 2). However, when we directly compare changes in SL as a measure for ice volume against CO2 climate forcing (ΔRCO2) from our records (Fig. 2), we find that, between eMPT and LP260, ice sheet mass increased progressively more per CO2 lowering, thereby increasing the SL–ΔRCO2 slope in Fig. 2. This suggests an increase in ice sheet sensitivity to CO2 forcing across the MPT, with the caveat that eMPT may not fully capture pre-MPT conditions, although it agrees with the longer-term record of Hönisch et al. (11). This finding is robust, regardless of which SL reconstruction is used (Fig. S7); in all cases, the SL to ΔRCO2 relationships appear to be linear, with increasing slopes from eMPT to LP260. The steepening relationship is also evident when regressing δ11B to δ18O relationships, with both isotope ratios measured on the same sample material (Fig. S8). Using the SL record with the best coverage of both intervals, relative SL from the Mediterranean Sea (4), we estimate 25 ± 3 and 45 ± 5 m of SL lowering for each 1-Wm−2 reduction in radiative forcing during eMPT and LP260, respectively. Such a pronounced increase in sensitivity implicates a change in ice sheet dynamics as predicted by the regolith hypothesis (8, 9) or the establishment of marine-based ice sheet margins in East Antarctica (10) (scenario 1 in Fig. 2).

The observed changes in the SL to ΔRCO2 relationships contain elements of both end member scenarios shown in Fig. 2 , in which a greater slope is possibly related to changes internal to the ice sheets (scenario 1) and amplified glacial to interglacial CO2 climate forcing is linked (this study) to increased glacial dustiness that causes enhanced Southern Ocean iron fertilization (scenario 2). Therefore, we propose a hybrid scenario (Fig. 2) that incorporates both heightened ice sheet sensitivity to CO2 forcing and dust-driven ocean sequestration of CO2 to represent the observed climate system change across the MPT.

First, we propose that—independent of orbital and CO2 forcing—a process internal to the climate system yielded greater glacial buildup of ice sheets [e.g., regolith removal (8) or ice sheet phase-locking (10)]. Second, we infer that larger ice sheets led to increased glacial atmospheric dustiness (31, 32), either directly through SL lowering or indirectly because of atmospheric cooling, drying, and/or changes in surface winds. This, in turn, induced glacial iron fertilization of the Sub-Antarctic Zone of the Southern Ocean, thereby effecting the 20- to 40-μatm increase in the amplitude of the G-IG CO2 cycles documented here (Fig. 4) (11). In our hybrid scenario, the positive climate–dust–CO2 feedback is required to (i) drive additional ice sheet growth and (ii) stabilize those ice sheets during the critical orbital phase of rising obliquity, ensuring the survival of ice sheets beyond single obliquity cycles. Therefore, regardless of the mechanism that served as the initial MPT trigger, our findings further illustrate the exquisite coupling that exists in the Earth System between climate change, ice sheet mass, and the polar ocean mechanisms that regulate G-IG CO2 change.

Materials and Methods

Globigerinoides ruber white sensu stricto (300–355 μm) were picked from sediments from ODP 999A (Fig. S1), and the age model was constructed by benthic oxygen isotopes from the same samples and X-ray fluorescence scanning data. Samples were measured for boron isotope composition using a Thermo Scientific Neptune multicollector inductively coupled plasma mass spectrometer at the University of Southampton according to methods described elsewhere (18). Analytical uncertainty is given by the external reproducibility of repeat analyses of Japanese Geological Survey Porites coral standard at the University of Southampton and is typically <0.2‰ (at 95% confidence). Metal element–calcium ratios (Mg, B, Al) were analyzed using Thermo Element 2XR inductively coupled plasma mass spectrometer at the University of Southampton. Here, these data are used to assess adequacy of clay removal (Al/Ca < 100 μmol/mol) and to generate down core temperature estimates. CO2 was calculated using a Monte Carlo approach (10,000 replicates) with estimates of salinity and alkalinity using a flat probability spanning a generous range (34–37 psu and 2,100–2,500 μmol/kg, respectively). A normal distribution around proxy data was used for all other input variables (temperature, pH, δ11Bsw, δ11Bforam) ( has full details). The CO2 record was then probabilistically assessed using a Monte Carlo approach that considers uncertainties in both age and CO2 values and that preserves the stratigraphy of the record, which minimizes age uncertainty in a relative sense between samples (shown as an envelope in Figs. 1 and 3). Each of 2,000 Monte Carlo iterations involved independent random resampling of each sample within its x and y uncertainty distributions. The stratigraphic constraint prevents age reversals in this resampling procedure. Linear interpolation was performed between resampled points, and the distribution of values thus generated was analyzed per time step for the modal value and its 95% probability interval as well as the 95% probability envelope of data in the sampled distribution (using the 2.5th and 97.5th percentiles). Because uncertainties in both x and y directions are considered, the record of probability maxima (modes) gives a smoothed representation of the record, with quantified uncertainties ().

Inverse carbon cycle modeling was carried out using the CYCLOPS model (23), with the forward model forcing derived from pertinent paleoceanographic records (25–27) and the forcing scaling parameters inverted to minimize model misfit with respect to the ice core CO2 record of the last 800 ky. Significant linear correlation with and matching spectral content to our boron isotope-based CO2 data confirm the skill of the model inversion (Fig. S6). Detailed statistical analysis is carried out to identify and quantify changes in absolute glacial and interglacial CO2 as well as the G-IG CO2 range from the model inversion results, our high-resolution CO2 data, and some previous datasets (11, 17) that are not well dated or lack the required temporal resolution for comparison in the time and/or frequency domains. This analysis is based on estimation of the population means of cumulative probability density of glacial and interglacial subsamples, which were selected based on either available benthic foraminiferal δ18O or CO2 rank (Fig. 4). Factorial analysis of the validated model allows for the mechanistic attribution to sub-Antarctic iron fertilization of glacial stage-specific CO2 reduction associated with the MPT interval (Fig. S5, Bottom), which is the pattern that we identified as common between model and all three empirical datasets. More detailed descriptions of inverse modeling and model/data cross-validation and statistical quantification of CO2 change can be found in and , respectively.