# Atmospheric CO2 during the Mid-Piacenzian Warm Period and the M2 glaciation

Elwyn de la Vega, Thomas Chalk, Paul Wilson, Ratna Priya Bysani, Gavin Foster

► To cite this version:

Elwyn de la Vega, Thomas Chalk, Paul Wilson, Ratna Priya Bysani, Gavin Foster. Atmospheric CO2 during the Mid-Piacenzian Warm Period and the M2 glaciation. Scientific Reports, 2020, 10 (1), pp.11002. 10.1038/s41598-020-67154-8. hal-04699109

HAL Id: hal-04699109

https://hal.science/hal-04699109v1

Submitted on 19 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Atmospheric CO₂ during the Mid-Piacenzian Warm Period and the M2 glaciation

Elwyn de la Vega✉, Thomas B. Chalk, Paul A. Wilson, Ratna Priya Bysani & Gavin L. Foster

The Piacenzian stage of the Pliocene (2.6 to 3.6 Ma) is the most recent past interval of sustained global warmth with mean global temperatures markedly higher (by ~2–3 °C) than today. Quantifying CO₂ levels during the mid-Piacenzian Warm Period (mPWP) provides a means, therefore, to deepen our understanding of Earth System behaviour in a warm climate state. Here we present a new high-resolution record of atmospheric CO₂ using the δ¹¹B-pH proxy from 3.35 to 3.15 million years ago (Ma) at a temporal resolution of 1 sample per 3–6 thousand years (kyrs). Our study interval covers both the coolest marine isotope stage of the mPWP, M2 (~3.3 Ma) and the transition into its warmest phase including interglacial KM5c (centered on ~3.205 Ma) which has a similar orbital configuration to present. We find that CO₂ ranged from 394±³⁴ ppm to 330±¹⁴ ppm: with CO₂ during the KM5c interglacial being 391±³⁰ ppm (at 95% confidence). Our findings corroborate the idea that changes in atmospheric CO₂ levels played a distinct role in climate variability during the mPWP. They also facilitate ongoing data-model comparisons and suggest that, at present rates of human emissions, there will be more CO₂ in Earth's atmosphere by 2025 than at any time in at least the last 3.3 million years.

The Pliocene Epoch (2.588 to 5.3 Ma) was a time when global temperatures were ~3 °C warmer than the pre-industrial¹ and sea level was ~20 m higher than present²,³, largely due to the presence of smaller Greenland and Antarctic ice sheets². Given that many other tectonic boundary conditions were similar (but not identical⁴) to today, the Pliocene Epoch, and the mid-Piacenzian warm period (mPWP; 3 to 3.3 Ma) in particular, are useful targets for climate model validation studies (e.g. refs. ⁴⁻⁷). Phase 1 of the Pliocene Model Intercomparison Project (PlioMIP), found an overall agreement between climate model simulations of mPWP surface temperatures and the available data when run with a CO₂ of 405 ppm¹. The North Atlantic and Pacific Oceans, however, are areas of consistently poor data-model agreement¹,⁶. Haywood et al.⁴ noted that the experiment design of PlioMIP simulations precluded a determination of whether this result is attributed to poor model performance or poor data quality. This deficiency in experimental design, in part, resulted from the data collection time interval being from 3 to 3.3 Ma which spans a number of orbital climate cycles, whereas the models were run for less than 1,000 years with an invariant modern orbit⁵.

To address this weakness, PlioMIP2 (phase 2)⁵ part of the model evaluations⁸ feeding into the 6th Assessment Report (AR) for the Intergovernmental Panel on Climate Change IPCC AR6, will focus on the KM5c interglacial interval at ~3.205 Ma which has a close-to-modern orbital configuration⁴. Other data compilation efforts with future data-model comparisons in mind, have targeted the interval from 3.3 Ma to 3.205 Ma¹,⁹ because this also includes marine isotope stage M2, a prominent anomalously cold marine isotope stage (MIS) that provides a cold-to-warm transition within the overall warm background climate state of the late Pliocene (Fig. 1). A major limitation of these latest efforts, however, is the lack of data on atmospheric CO₂ during these rather narrow time intervals⁹⁻¹¹. For instance, although ~40 δ¹¹B-based determinations of CO₂ are available for the 3 to 3.3 Ma window¹²⁻¹⁴, the M2 to KM5 transition remains poorly sampled and disagreement between datasets from different analytical techniques persist (Fig. 1). Furthermore, the available δ¹¹B-CO₂ data exhibit relatively large short-term variability, which could signal orbital cyclicity (41 kyrs) aliased by low sampling resolution. Data from other CO₂ proxy systems such as stomata and palaeosol δ¹³C, are lacking in this interval and it has recently been shown that the marine-based alkenone-δ¹³C-CO₂ proxy underestimates CO₂ levels in the Pliocene¹⁵,¹⁶ therefore limiting its usefulness to providing a minimum CO₂ during the mPWP of >270 ppm¹⁶.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top panel: Current CO₂ estimates from boron isotopes across the Plio-Pleistocene boundary. *G. ruber* data in red circles from site 999 (Martinez-Boti et al. ¹³), *T. sacculifer* in red squares from site 999 (Seki et al. ¹⁴), blue squares from site 926 (Sosdian et al. ²⁰) and pale red squares from site 999 (Bartoli et al. ¹²). Bottom Panel: δ¹⁸O from benthic foraminifer *Cibicidoides wuellerstorfi* at ODP Site 999 (blue circles) with a 5 point running mean (this study and ref. ⁵³) compared to the benthic isotope stack of ref. ¹⁹. M2 glacial and early Pleistocene strong glacials are highlighted in blue for context. Interglacial KM5 and KM5c are highlighted in yellow and orange, respectively. Note that there are no estimates for the M2 glacial and very few across the mid-Piacenzian warm period (mPWP), the low resolution of previous studies makes pin-pointing individual interglacials such as the KM5c future analogue difficult. These studies also differ in their estimates of Mid-Piacenzian CO₂ ¹²⁻¹⁴,²⁰.</figcaption></figure>

To address this data deficiency, we developed δ¹¹B-based CO₂ estimates from Ocean Drilling Program (ODP) Site 999 in the Caribbean (Supplementary Fig. 1) at a resolution of 1 sample per 3–6 kyr for the time interval 3.15 to 3.35 Ma, encompassing the M2 glaciation and the KM5 interglacial (including KM5c). Although focusing predominantly on the mixed layer dwelling planktic species *Globigerinoides ruber* (45 new data points, 63 in total), we also present new measurements of *Trilobus sacculifer* (2 new data points, 5 total) on the same samples to provide a check on the consistency of the δ¹¹B-pH calibration for *G. ruber* which has been recently called into question¹⁷.

## Results

Our new high-resolution CO₂ record is shown in Fig. 2 (and Supplementary Fig. 2 and Supplementary Table 2) and is consistent with earlier studies¹²⁻¹⁴ in showing that CO₂ was higher than the pre-industrial during the mPWP (mean=360 ppm). Our more detailed record reveals that CO₂ variations ranged from 330±¹⁴ to 394±⁹ ppm (based on the mean and distribution of CO₂ in the <25% and >75% interquartile range during the whole studied period; ref. ¹⁸), with a peak-to-trough range of 56–75 ppm determined by a Welch T-test of the data within the upper and lower quartiles (at 95% confidence; p < 0.01). From a comparison with benthic δ¹⁸O data from ODP 999 (Supplementary Table 2) and the δ¹⁸O stack¹⁹ (LR04) we observe that cold marine isotope stages (e.g. KM2; Fig. 2) are typically closely associated with low CO₂ and warm stages with high CO₂ levels. However, during the prominent M2 cold stage and through the warming out of M2, CO₂ appears to lag benthic δ¹⁸O by ~10 kyr (Supplementary Fig. 5). This lag is not attributable to age model uncertainty because it is present when comparing δ¹¹B-derived CO₂ and benthic δ¹⁸O from the same samples (Fig. 2). CO₂ during the interglacial KM5c, determined using the mean of the δ¹¹B of the five samples in this interglacial, is estimated at 391±³⁰ ppm (at 95% confidence). Using a broader time window (10 to 15 ky) for the KM5 interglacial moderately alters the estimate (see Supplementary Table 1).

Our estimates of borate ion, pH and CO₂ derived from the δ¹¹B of *T. sacculifer* (without sac-like final chamber) from ODP 999 (our new data, and ref. ¹⁴) and ODP 926²⁰ overlap well with those based on *G. ruber* from ODP 999 in the mPWP (Fig. 2 and Supplementary Fig. 3). Calculated CO₂ does not exhibit substantive inter-species offset with a mean difference of 10 ± 29 ppm with no consistent bias towards higher or lower pH/CO₂.

**Accuracy of δ¹¹B-CO₂ from G. ruber in the Plio-Pleistocene.** Ref. ¹⁷ suggested that the δ¹¹B-pH calibration of *G. ruber* may have evolved through time and that *G. ruber* may suffer from morphotype-differences in δ¹¹B-derived pH estimates for the surface Pliocene ocean resulting in underestimates of the true pH and corresponding overestimates of true CO₂. The principal evidence presented for this assertion was the disagreement between Pliocene CO₂ calculated using the *T. sacculifer* data of Bartoli *et al.*¹² and the *G. ruber* data from Martinez-Boti *et al.*¹³ (Fig. 1). As shown here, when *T. sacculifer* and *G. ruber* are compared from the same samples and measured with the same analytical technique (in this case MC-ICPMS) there is no significant offset between the methods in terms of reconstructed borate ion δ¹¹B, pH or CO₂ (Fig. 2 and Supplementary Fig. 3). This finding suggests that the *G. ruber* δ¹¹B-pH calibration has not evolved through time and therefore that the pH (and hence CO₂) we reconstruct here is an accurate record of the surface water carbonate system parameter. This finding indicates that disagreements between published Pliocene δ¹¹B-based datasets¹²,¹³,²⁰,²¹ are likely attributable to either: (i) sampling driven aliasing due to the relatively large short-term CO₂ variability in the mPWP (e.g.
<figure><img src="image_4.png" /></figure>

**Figure 2.** Top panel: Red circles and lines show δ¹¹B-derived CO₂ data from *Globigerinoides ruber* at ODP Site 999 (this study and Martinez-Boti *et al.*¹³, Chalk *et al.*¹⁸), red squares are *Trilobatus sacculifer* at ODP 999 (this study and Seki *et al.*¹⁴), purple squares are *T. sacculifer* from ODP 668 (Honisch *et al.*²³) and blue squares are *T. sacculifer* from ODP 926 (Sosdian *et al.*²⁰). Black solid line shows ice core-derived CO₂ from ref. ⁵⁸. Left; Late Pleistocene CO₂ from boron isotopes¹⁴,¹⁸,²³,⁵² and ice core data. Also shown are CO₂ projections in line with RCP8.5 at current emission rates to the year 2040 (black broken line). Middle column; MPT CO₂ estimates¹⁸,²³ including disturbed ice estimates²⁴,²⁵ (Note: age adjusted for scale). Right; mPWP estimates of CO₂ (this study combined with Martinez-Boti *et al.*¹³), new data from *T. sacculifer* is shown in red squares and shows no offset from *G. ruber* estimates. Second panel: Time periods as above, LR04 and ODP 999 δ¹⁸O from *C. wuellerstorfi*¹⁸,¹⁹,⁵³. Third panel: Iron mass accumulation rate from the Southern Atlantic ODP Site 1090²⁸. Fourth panel: % Northern Component Water (NCW) estimated from δ¹³C in benthic foraminifera (grey) and εₐ₉₉ₐ₉₉₄ from fish debris (dark green) in the deep North Atlantic (core U1313³¹). Note the lag of ocean circulation and CO₂ relative to the M2 glaciation³¹.

Figure 2), (ii) differences in core site location between studies (with possible different local CO₂ disequilibrium), (iii) the lack of a comparison between species on exactly the same sample, or (iv) the well-documented analytical offset between MC-ICPMS and negative-ion thermal ionization mass spectrometry (NTIMS; see refs. ²¹,²²). We note that, if the offset is attributable to analytical issues, the agreement between data generated by both methodologies for the younger Pleistocene time slices examined here (Fig. 2) confirms the suggestion of Sosdian *et al.*²⁰ that the *T. sacculifer* dataset of Bartoli *et al.*¹² (measured with NTIMS) requires an additional analytical correction (see ref. ¹⁶ for details), beyond that currently used for the NTIMS δ¹¹B datasets of refs. ¹⁷,²³.

**mPWP CO₂ cycles – variability and causes.** Our new high-resolution data set clearly demonstrates that the mPWP is an interval of relatively high CO₂, with a mean of 367 ppm whereas lower values are observed during the mid- and late-Pleistocene (262 ppm and 241 ppm¹⁸, respectively Fig. 2). Our dense sampling permits, for
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Distribution plots showing the distribution of absolute CO₂ (left), climate forcing from CO₂ (middle), and δ¹⁸O (right) across three time intervals (where CO₂ cycles are consistent); Late Pleistocene 0–250 ka (LP250, bottom), Mid-Pleistocene (MPT, 1050–1500 ka) and Middle Piacenzian including M2 (mPWP, 3–3.3 Ma). The white dots represent the median values with the quartile ranges shown in the thick and thin black lines. The width of the distribution shows the density of data in each range. Climate forcing (ΔR<sub>CO₂</sub>) is expressed relative to pre-industrial CO₂ value of 278 ppm. The effect of excluding M2 in the mPWP interval has a negligible effect on the distribution of CO₂ and CO₂ forcing but a visible impact on δ¹⁸O.</figcaption></figure>

the first time, an assessment of CO₂ variability during mPWP on orbital timescales, although the length of our record still precludes reliable time series analysis. The CO₂ cycles (and derived climate CO₂ forcing) we observe in the mPWP (including the M2 glaciation) are similar, yet smaller in amplitude, than those evident in δ¹¹B-CO₂ data from the Late Pleistocene (LP = 0–250 kyr; ref. ¹⁸) and early mid-Pleistocene Transition (eMPT = 1050–1250 kyr; Figs. 2 and 3). The variation in climate forcing between these different intervals correlates with the magnitude of δ¹⁸O variability, although there is an apparent increased response in δ¹⁸O during the late Pleistocene (δ¹⁸O range ~ 2‰) relative to the mPWP (δ¹⁸O range ~ 0.5–1‰ Fig. 3; ref. ¹³) likely reflecting the increased influence of ice-volume change on δ¹⁸O following the northern hemisphere glaciation¹³.

We also compare the δ¹¹B-derived CO₂ data with ice-core based CO₂ for the Pleistocene intervals¹⁸,²³ (Supplementary Fig. 4). We expect the range of CO₂ variability from the various Antarctic ice cores to be narrower than our marine-based records because of the greater precision associated with the ice core records (± 6 ppm vs. ±20–30 ppm) and the CO₂ data from blue ice from the mid-Pleistocene may not capture a full climate cycle²⁴,²⁵. Despite these caveats, and as discussed elsewhere¹⁸,²⁵, there is good agreement between the ice-core and δ¹¹B-derived CO₂, providing confidence in the accuracy of the distribution (and absolute values) we determine here for the mPWP.

Three phenomena are suggested to exert a major control over glacial-interglacial CO₂ change in the Late Pleistocene record: (i) changes in stratification south of the Antarctic Polar Front influencing the ventilation of CO₂-rich circumpolar deep water;²⁶ (ii) variations in the magnitude of dust-borne Fe fertilization in the sub-Antarctic that influences the strength of the biological pump in this region;²⁷,²⁸ (iii) changes in the extent to which southern component deep water fills the North Atlantic increasing deep ocean carbon storage²⁹. Ref. ³⁰ proposed that polar waters in the North Pacific and the entire Southern Ocean (SO) were destratified prior to the onset of Northern Hemisphere Glaciation at 2.7 Ma. This mechanism may contribute to elevated CO₂ during the mPWP³⁰, but it would also rule out water mass ventilation in the high-latitude SO as an important driver of the CO₂ cycles we observe in the mPWP. Both the accumulation and variability of dust-borne Fe at ODP Site 1090 (Fig. 2) were reduced during the mPWP to a fraction of that observed in the other time intervals examined (~factor of four reduction during the mPWP, compared to the late Pleistocene), suggesting that dust supply did not play a major role in generating the observed CO₂ cycles in the mPWP.

Chemical stratification of the North Atlantic Ocean, due to incursions of southern component water (SCW), during the mPWP is suggested by gradients in ε<sub>Nd</sub> and δ<sup>13</sup>C from ref. ³¹, albeit at reduced magnitude compared to the latest Pleistocene and eMPT. Ref. ³¹ showed that the SCW-signal was particularly marked after the M2 glacial maximum as shown by reduced % NCW (northern component water, Fig. 2), and also observed that ε<sub>Nd</sub> lagged δ<sup>18</sup>O by 10 kyr during the onset of the M2 glaciation, as is also evident in our new CO<sub>2</sub> data (Fig. 2, and supplementary Fig. 5). In general, those intervals of the mPWP where SCW contributed significantly to the waters of the deep North Atlantic were characterized by low CO<sub>2</sub> (and vice versa; Fig. 2), yet with variable leads and lags (supplementary Fig. 5). This association is consistent with the suggestion³¹,³² that glacial expansion of a CO<sub>2</sub>-charged SCW reservoir plays at least a first order role in driving orbital CO<sub>2</sub> cycles, even prior to the onset of northern hemisphere glaciation at 2.7 Ma.

**M2 Glaciation – the role of CO₂.** Marine Isotope Stage M2 (at 3.3 Ma) is an anomalously cold stage clearly evident in the LR04 δ¹⁸O stack¹⁹ and many other temperature records, such as arctic air temperature³³ and sea surface temperature³⁴,³⁵ (including site ODP 999, supplementary Fig. 6) It is also often described as a failed attempt at Northern Hemisphere Glaciation³⁴,³⁶,³⁷, that eventually initiated ~600 kyr later at ~2.6 Ma. Refs. ¹³,¹⁴ used δ¹¹B-CO₂ to suggest that CO₂ dropped below 280 ppm for the first time during the intensification of Northern Hemisphere
Glaciation (iNHG), consistent with the suggestion that 280 ppm CO₂ is an important threshold below which extensive glaciation of continents in the northern hemisphere develop³⁸. Our new data reveal that while the lower bound of the error envelop in our smoothed CO₂ record approaches this value, atmospheric CO₂ during M2 is unlikely to have crossed this threshold. Furthermore, in our records, CO₂ variability associated with M2 lags δ¹⁸O by 10 kyr, which also means that minimum CO₂ is not coincident with minimum northern hemisphere orbital forcing (Supplementary Fig. 5). Additional support for CO₂ playing a secondary role in M2 is that other periods of low CO₂ are evident in the mPWP (Fig. 2). For instance, the KM2 glaciation is clearly evident in the benthic δ¹⁸O data, Mg/Ca-SST at site 999 (Fig. 2 and Supplementary Fig. 6) and our CO₂ record, but is not considered to be a major glacial interval¹⁹ (Fig. 2). This therefore suggests that other boundary conditions, such as orbital con-figuration³⁹ (Supplementary Fig. 5) were perhaps dominant in triggering the M2 cold stage.

**M2 Glaciation - CO₂ lags δ¹⁸O.** Minor dephasing between εNd, δ¹³C and benthic δ¹⁸O has been observed in the late Pleistocene⁴⁰, with carbon budget change lagging δ¹⁸O and preceding the change in ocean circulation, but with smaller lags than we observe during M2 (~2.5 ky vs. ~10 ky).

Similarly, while variations in atmospheric CO₂ are not entirely in phase with ice volume changes over the late Pleistocene, the leads-lags are small⁴¹ (<3 ky) and are not readily observed when comparing δ¹¹B-derived CO₂ and benthic foraminiferal δ¹⁸O time series (e.g. Figure 2 left panels).

This therefore either requires the operation of subtly distinct carbon-cycle dynamics during the Pliocene, and M2 in particular, or implies some other driver operated to explain the observed ~10 kyr lag of CO₂ during M2.

One possible non-carbon cycle driver for the observed lag at M2 could be a preservation bias in our data because partial dissolution of planktic foraminiferal tests drives δ¹¹B towards more negative isotopic composition (lower pH, higher CO₂) in some species¹⁴,⁴². However, while our chosen species for pH/CO₂ reconstruction, *G. ruber*, is known to be relatively susceptible to partial dissolution on the seafloor⁴³, its δ¹¹B signal has been observed to be robust¹⁴,⁴⁴. Furthermore, a recent study⁴⁵ of test weight and fragmentation at ODP 999 showed that tests become better preserved during M2 (Supplementary Fig. 7), as observed during glacial periods of the late Pleistocene⁴⁶. This is inconsistent with the high CO₂ observed during the descent into M2 maximum (especially given fragmentation and δ¹⁸O are in phase, Supplementary Fig. 7), ruling out an effect of partial dissolution on our CO₂ reconstruction and observed lag.

An alternative explanation for the delayed pH change observed at ODP 999 during M2 lies in local water mass changes during this glacial episode. Although the Central American Seaway was closed to deep water by this time³⁶, it is possible that limited exchange of surface water was still occurring around M2³⁴ and the early Pleistocene⁴⁷. Temperature data from Mg/Ca in *T. sacculifer*³⁴ and *G. ruber* at site 999 (this study, Supplementary Fig. 6, data from the East Equatorial Pacific site 1241⁴⁸ is also shown for comparison) show a cooling prior to the M2 maximum. A connection between the East Equatorial Pacific and the Caribbean, or local intensification in upwelling could have brought cold, nutrient- and carbon-rich waters (low Mg/Ca, low pH) to Site 999, potentially explaining the apparent high CO₂ at the inception of M2. However, because no noted decline in the abundance of *G. ruber* occurred, given this is a species known to favour oligotrophic conditions⁴⁹, it seems unlikely this was accompanied by a significant influx of such water. Also, a mechanism for increased influx of EEP water during sea level regression is lacking. Importantly, during KM5c (3205ky), the CAS likely remained closed impeding the influx of Pacific waters to site 999, as shown by elevated temperature (relatively to the M2 interval) in the Caribbean (Supplementary Fig. 6), hence local changes in hydrography are unlikely to have unduly impacted our CO₂ estimates for this central interval.

The cause of the apparent lag of CO₂ compared to δ¹⁸O during the inception of M2 therefore remains uncertain, but because this lag is also present in North Atlantic εNd records³¹ we favour a carbon cycle-based interpretation. A Southern Ocean-driven explanation for the CO₂ lag during M2 is perhaps indicated by the observation that the tail of the CO₂ decline during M2 is out of phase with decline in northern hemisphere insolation, but apparently in-phase with insolation decline at 65°S (i.e. offset by one half precession cycle, Supplementary Fig. 5).

**Past to future – when will we exceed mPWP CO₂ levels?** Atmospheric CO₂ has been increasing since the industrial revolution from a background of ~280 ppm, reaching 411 ppm in 2019. The mPWP is often cited as being the last time CO₂ levels were as high as today¹³, although we note methane and other greenhouse gases remain poorly constrained and may contribute to extra radiative forcing⁵⁰.

Our refined view of CO₂ during this interval however reveals that in terms of the mean (367 ppm), mPWP values were exceeded in the mid-1990s. Our upper quartile range (394±³⁴ ppm) suggests that CO₂ during the mPWP is very likely to have been ≤ 427 ppm (using the distinctions of the IPCC). Atmospheric CO₂ rose by 2.5 ppm from 2017 to 2018, if this rate is sustained, our new data indicate that CO₂ will surpass even the highest mPWP values within the next 5 to 6 years (2025–2026).

**Conclusions**

Our new δ¹¹B-data for the mPWP provide a tightly constrained and robust estimate of CO₂ during KM5c of 391±³⁰ to 391±³⁴ ppm (at 95% confidence), and documents CO₂ variability during the mPWP from 330±¹⁴ to 391±³⁴ ppm. This extended view suggests that changes in ocean carbon storage may play an important role in natural variability in CO₂ on orbital timescales in warmer than present climate states. By better constraining the upper levels of CO₂ during the mPWP we conclude that, given current rates of CO₂ increase, we will very likely surpass mPWP values within 4 to 5 years, meaning that, by 2024/2025 levels of atmospheric CO₂ will be higher than any point in the last 3.3 million years.
## Methods

ODP Site 999 is located in the Caribbean Sea and has been reliably used to reconstruct past atmospheric CO₂ in a number of previous studies¹⁸,⁵¹,⁵². Air-sea disequilibrium for CO₂ in the surface waters in the region of ODP Site 999 are close to 0 (+20 ppm) today and ref. ¹³ suggest this remained the case for at least the last 3.3 million years. An age model for the interval 3.15 to 3.35 Ma studied here was constructed by analyzing the epibenthic foraminifera *Cibicidoides wuellerstori* for δ¹⁸O in each sample, combining these data with similar data from ref. ⁵³, and tuning the combined record to the LR04 benthic δ¹⁸O stack¹⁹ (Fig. 1). From each sample, ~120 individual tests of the mixed layer dwelling species *Globigerinoides ruber* white sensu stricto (300–355 µm) were hand separated, clay removed, oxidatively cleaned and analyzed for boron isotopic (δ¹¹B) and trace element composition (e.g. Mg/Ca, Al/Ca) at the University of Southampton using well-established procedures²¹,⁵². These data were converted to pH using the *G. ruber* δ¹¹B-pH calibration of ref. ⁵¹ and a modern δ¹¹B of seawater⁵⁴ (39.6‰). Sea surface temperatures (SST) were derived from each samples Mg/Ca ratio using the calibrations of ref. ⁵⁵ corrected for the changing Mg/Ca for seawater following ref. ²⁰ based on ref. ⁵⁶.

It was recently suggested¹⁷ that the δ¹¹B-pH calibration for *G. ruber* may have varied over the last 3.5 million years. We therefore also analyse *T. sacculifer* (300–355 µm) from two samples from ODP 999 to combine with existing *T. sacculifer* data¹⁴,²⁰,⁵² to provide a test of the *G. ruber* δ¹¹B-pH calibration of ref. ⁵¹. Following previous studies¹³, uncertainty relating to the following factors, including uncertainty in the δ¹¹B-pH calibration, were propagated using a Monte Carlo approach (n = 10,000): temperature (±1.5 °C; 2σ), salinity (± 3 psu, 2σ), δ¹¹Bsw (± 0.2, 2σ; refs. ²⁰,⁵⁷), analytical uncertainty (± 0.12–0.3‰, 2σ; see Supplementary Table 2).

In order to calculate atmospheric CO₂ from the reconstructed surface water pH we use dissolved inorganic carbon (DIC) from the reconstructions of ref. ²⁰. Uncertainty in this term is also propagated into CO₂ uncertainty using a Monte Carlo approach (n = 10,000), but with a uniform distribution encompassing the range of the reconstructed DIC for our study interval (1765 to 1840 µmol/kg). CO₂ was then determined by the maximum probability of all 10,000 realisations.

To obtain CO₂ during the KM5c interglacial the mean and associated uncertainty of the reconstructed borate ion δ¹¹B and SST of the data lying within +/- 7 kyr of the peak of KM5c (n = 5) was determined. In our age model we defined that the peak of KM5c was at 3212 kyr, although expanding the window to ±15 kyr, or uncertainty in the identification of the peak of KM5c to within 10 kyr, did not have a significant effect on our calculated mean (Supplementary Table 1).

This average was then used to calculate pH and CO₂, but only considering the uncertainty in δ¹¹Bsw and DIC detailed above. This approach was chosen because it allows the random uncertainties arising from SST and δ¹¹B measurement error to be reduced through replication, whilst still realistically accounting for the systematic uncertainties in δ¹¹Bsw and DIC.

In order to place our new data into a Plio-Pleistocene context, we recalculated the δ¹¹B data of ref. ¹³ from 3.0 to 3.3 Ma using the same methodology outlined here and combined it with our new data. The difference in the method is the choice of the second carbonate parameter where constant modern alkalinity is used in ref. ¹³ whereas we use DIC²⁰ in this study. These new and published data are then combined into a single record with an average temporal resolution of 1 sample every 4 kyr. This was then interpolated to 1 kyr resolution and a 6-point running mean was used to reduce the influence of analytical uncertainty on our reconstructed CO₂ record. Uncertainty in the smoothed record was determined using the output of the original Monte Carlo simulation.

CO₂ forcing was calculated relatively to preindustrial CO₂ (278 ppm) and is defined as:

$$ (R_{CO_2}) = \ln \frac{CO_2}{CO_{2preind}} $$

Received: 30 December 2019; Accepted: 18 May 2020;
Published: 9 July 2020

## References

1. Dowsett, H. J. et al. Sea surface temperature of the mid-Piacenzian ocean: a data-model comparison. *Scientific reports* 3, https://doi.org/10.1038/srep02013 (2013).

2. Dutton, A. et al. Sea-level rise due to polar ice-sheet mass loss during past warm periods. *science* 349, aaa4019 (2015).

3. Dumitru, O. A. et al. Constraints on global mean sea level during Pliocene warmth. *Nature* 574, 233–236 (2019).

4. Haywood, A. et al. Large-scale features of Pliocene climate: results from the Pliocene Model Intercomparison Project. *Climate of the Past* 9, 191–209 (2013).

5. Haywood, A. et al. The Pliocene model intercomparison project (PlioMIP) phase 2: scientific objectives and experimental design. *Climate of the Past* 12, 663–675 (2016).

6. Haywood, A. M., Dowsett, H. J. & Dolan, A. M. Integrating geological archives and climate models for the mid-Pliocene warm period. *Nature communications* 7, 10646 (2016).

7. DeConto, R. M. & Pollard, D. Contribution of Antarctica to past and future sea-level rise. *Nature* 531, 591 (2016).

8. Eyring, V. et al. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. *Geoscientific Model Development (Online)* 9 (2016).

9. McClymont, E. et al. Lessons from a high CO₂ world: an ocean view from ~3 million years ago. *Climate of the Past Discussion*, in review (2020).

10. Tan, N. et al. Exploring the MIS M2 glaciation occurring during a warm and high atmospheric CO₂ Pliocene background climate. *Earth and Planetary Science Letters* 472, 266–276 (2017).

11. Tan, N. et al. Modeling a modern-like pCO₂ warm period (Marine Isotope Stage KM5c) with two versions of an Institut Pierre Simon Laplace atmosphere-ocean coupled general circulation model. *Climate of the Past* 16, 1–16 (2020).

12. Bartoli, G., Hönisch, B. & Zeebe, R. E. Atmospheric CO₂ decline during the Pliocene intensification of Northern Hemisphere glaciations. *Paleoceanography* 26, PA4213 (2011).

13. Martínez-Boti, M. et al. Plio-Pleistocene climate sensitivity evaluated using high-resolution CO₂ records. *Nature* 518, 49 (2015).
14. Seki, O. et al. Alkenone and boron-based Pliocene pCO2 records. *Earth and Planetary Science Letters* **292**, 201–211 (2010).

15. Badger, M. P., Schmidt, D. N., Mackensen, A. & Pancost, R. D. High-resolution alkenone palaeobarometry indicates relatively stable p CO2 during the Pliocene (3.3–2.8 Ma). *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **371**, 20130094 (2013).

16. Badger, M. P. et al. Insensitivity of alkenone carbon isotopes to atmospheric CO2 at low to moderate CO2 levels. *Climate of the Past* **15**, 539–554 (2019).

17. Dyez, K. A., Hönisch, B. & Schmidt, G. A. Early Pleistocene obliquity-scale pCO2 variability at~ 1.5 million years ago. *Paleoceanography and Paleoclimatology* **33**, 1270–1291 (2018).

18. Chalk, T. B. et al. Causes of ice age intensification across the Mid-Pleistocene Transition. *Proceedings of the National Academy of Sciences* **114**, 13114–13119 (2017).

19. Lisiecki, L. E. & Raymo, M. E. A Plio-Pleistocene stack of 57 globally distributed benthic δ 18 O records. *Paleoceanography* **20**, 1–17 (2005).

20. Sosdian, S. M. et al. Constraining the evolution of Neogene ocean carbonate chemistry using the boron isotope pH proxy. *Earth and Planetary Science Letters* **498**, 362–376 (2018).

21. Foster, G. L. et al. Interlaboratory comparison of boron isotope analyses of boric acid, seawater and marine CaCO3 by MC-ICPMS and NTIMS. *Chemical Geology* **358**, 1–14 (2013).

22. Farmer, J. R., Hönisch, B. & Uchikawa, J. Single laboratory comparison of MC-ICP-MS and N-TIMS boron isotope analyses in marine carbonates. *Chemical Geology* **447**, 173–182 (2016).

23. Hönisch, B., Hemming, N. G., Archer, D., Siddall, M. & McManus, J. F. Atmospheric carbon dioxide concentration across the mid-Pleistocene transition. *Science* **324**, 1551–1554 (2009).

24. Higgins, J. A. et al. Atmospheric composition 1 million years ago from blue ice in the Allan Hills, Antarctica. *Proceedings of the National Academy of Sciences* **112**, 6887–6891 (2015).

25. Yan, Y. et al. Two-million-year-old snapshots of atmospheric gases from Antarctic ice. *Nature* **574**, 663–666 (2019).

26. Jaccard, S. et al. Two modes of change in Southern Ocean productivity over the past million years. *Science* **339**, 1419–1423 (2013).

27. Hain, M. P., Sigman, D. M. & Haug, G. H. Carbon dioxide effects of Antarctic stratification, North Atlantic Intermediate Water formation, and subantarctic nutrient drawdown during the last ice age: Diagnosis and synthesis in a geochemical box model. *Global Biogeochemical Cycles* **24**, GB4023 (2010).

28. Martínez-Garcia, A. et al. Southern Ocean dust-climate coupling over the past four million years. *Nature* **476**, 312 (2011).

29. Curry, W. B. & Oppo, D. W. Glacial water mass geometry and the distribution of δ13C of ΣCO2 in the western Atlantic Ocean. *Paleoceanography* **20** (2005).

30. Sigman, D. M., Jaccard, S. L. & Haug, G. H. Polar ocean stratification in a cold climate. *Nature* **428**, 59 (2004).

31. Lang, D. C. et al. Incursions of southern-sourced water into the deep North Atlantic during late Pliocene glacial intensification. *Nature Geoscience* **9**, 375 (2016).

32. Sigman, D. M., Hain, M. P. & Haug, G. H. The polar ocean and glacial cycles in atmospheric CO2 concentration. *Nature* **466**, 47–55 (2010).

33. Brigham-Grette, J. et al. Pliocene warmth, polar amplification, and stepped Pleistocene cooling recorded in NE Arctic Russia. *Science* **340**, 1421–1427 (2013).

34. De Schepper, S. et al. Northern hemisphere glaciation during the globally warm early late Pliocene. *PloS one* **8**, e81508 (2013).

35. Herbert, T. D. et al. Late Miocene global cooling and the rise of modern ecosystems. *Nature Geoscience* **9**, 843 (2016).

36. Haug, G. H. & Tiedemann, R. Effect of the formation of the Isthmus of Panama on Atlantic Ocean thermohaline circulation. *Nature* **393**, 673 (1998).

37. Mudelsee, M. & Raymo, M. E. Slow dynamics of the Northern Hemisphere glaciation. *Paleoceanography* **20**, PA4022 (2005).

38. DeConto, R. M. et al. Thresholds for Cenozoic bipolar glaciation. *Nature* **455**, 652 (2008).

39. Laskar, J. et al. A long-term numerical solution for the insolation quantities of the Earth. *Astronomy & Astrophysics* **428**, 261–285 (2004).

40. Piotrowski, A. M., Goldstein, S. L., Hemming, S. R. & Fairbanks, R. G. Temporal relationships of carbon cycling and ocean circulation at glacial boundaries. *Science* **307**, 1933–1938 (2005).

41. Mudelsee, M. The phase relations among atmospheric CO2 content, temperature and global ice volume over the past 420 ka. *Quaternary Science Reviews* **20**, 583–589 (2001).

42. Hönisch, B. & Hemming, N. G. Surface ocean pH response to variations in pCO2 through two full glacial cycles. *Earth and Planetary Science Letters* **236**, 305–314 (2005).

43. Berger, W. H. Planktonic foraminifera: selective solution and the lysocline. *Marine Geology* **8**, 111–138 (1970).

44. Ni, Y. et al. A core top assessment of proxies for the ocean carbonate system in surface-dwelling foraminifers. *Paleoceanography* **22**, PA3212 (2007).

45. Todd, C., Schmidt, D., Robinson, M. & De Schepper, S. Planktic foraminiferal test size and weight response to the late Pliocene environment. *Paleoceanography and Paleoclimatology*, e20831 (2019).

46. Schmidt, M. W., Vautravers, M. J. & Spero, H. J. Western Caribbean sea surface temperatures during the late Quaternary. *Geochemistry, Geophysics, Geosystems* **7** (2006).

47. Groeneveld, J. et al. Glacial induced closure of the Panamanian Gateway during Marine Isotope Stages (MIS) 95–100 (~2.5 Ma). *Earth and Planetary Science Letters* **404**, 296–306 (2014).

48. Groeneveld, J; Tiedemann, R. Mg/Ca and sea surface temperature data for Site 202-1241. PANGAEA, https://doi.org/10.1594/PANGAEA.315654 (2005).

49. Lombard, F. et al. Modelling planktic foraminifer growth and distribution using an ecophysiological multi-species approach. *Biogeosciences Discussions* **8**, 853–873 (2011).

50. Hopcroft, P. O. et al. Polar amplification of Pliocene climate by elevated trace gas radiative forcing. *Proceedings of the National Academy of Sciences*, in press (2020).

51. Henehan, M. J. et al. Calibration of the boron isotope proxy in the planktonic foraminifera Globigerinoides ruber for use in palaeo-CO2 reconstruction. *Earth and Planetary Science Letters* **364**, 111–122 (2013).

52. Foster, G. Seawater pH, pCO2 and [CO2 – 3] variations in the Caribbean Sea over the last 130 kyr: A boron isotope and B/Ca study of planktic foraminifera. *Earth and Planetary Science Letters* **271**, 254–266 (2008).

53. Haug, G. H., Tiedemann, R., Zahn, R. & Ravelo, A. C. Role of Panama uplift on oceanic freshwater balance. *Geology* **29**, 207–210 (2001).

54. Foster, G., Pogge von Strandmann, P. A. & Rae, J. Boron and magnesium isotopic composition of seawater. *Geochemistry, Geophysics, Geosystems* **11** (2010).

55. Delaney, M. L., Bé, A. W. & Boyle, E. A. Li, Sr, Mg, and Na in foraminiferal calcite shells from laboratory culture, sediment traps, and sediment cores. *Geochimica et Cosmochimica Acta* **49**, 1327–1341 (1985).

56. Evans, D. & Müller, W. Deep time foraminifera Mg/Ca paleothermometry: Nonlinear correction for secular change in seawater Mg/Ca. *Paleoceanography* **27**, PA4205 (2012).

57. Greenop, R. et al. A record of Neogene seawater δ11B reconstructed from paired δ11B analyses on benthic and planktic foraminifera. *Climate of the Past* **13**, 149–170 (2017).

58. Bereiter, B. et al. Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present. *Geophysical Research Letters* **42**, 542–549 (2015).
**Acknowledgements**

This research was supported by the National Environmental Research Council grant NE/H006273/1 and Royal Society Wolfson Merit awards to P. A. W. and G. L. F. We thank J. Andy Milton for plasma mass spectrometry support, Agnes Michalik and Matthew Cooper for clean laboratory assistance, and Megan Spencer and Bastian Hambach for support in stable isotope analysis.

**Author contributions**

E.d.l.V., G.L.F. and T.B.C. wrote the manuscript. E.d.l.V. generated boron isotope and element data. R.P.B. and T.B.C. generated oxygen and carbon isotope data. P.A.W. provided input on the manuscript. All authors contributed to the interpretation of the data.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41598-020-67154-8.

**Correspondence** and requests for materials should be addressed to E.d.l.V.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_6.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020, corrected publication 2021