Eccentricity-paced atmospheric carbon-dioxide variations across the middle Miocene climate transition

The middle Miocene climate transition ~14 Ma marks a fundamental step towards the current “icehouse” climate, with a ~1 ‰ δO increase and a ~1 ‰ transient δC rise in the deep ocean, indicating rapid expansion of the East Antarctic Ice Sheet associated with a change in the operation of the global carbon cycle. The variation of atmospheric CO2 across the carbon-cycle perturbation has been intensely debated as proxy records of pCO2 for this time interval are sparse and partly contradictory. Using boron isotopes (δB) in planktonic foraminifers from drill site ODP 1092 in the South Atlantic, we show that long-term pCO2 variations between ~14.3 and 13.2 Ma were paced by 400 k.y. eccentricity cycles, with decreasing pCO2 at high eccentricity and vice versa. Our data support results from a carbon-cycle model study, according to which increased monsoon intensity at high eccentricity enhanced weathering and river fluxes in the tropics, resulting in increasing carbonate and organic carbon burial and hence decreasing atmospheric CO2. In this scenario, a combination of the eccentricity-driven climatic cycle and enhanced meridional deep-ocean circulation during Antarctic ice-sheet expansion may have both contributed to the pCO2 rise following Antarctic glaciation, acting as a negative feedback on the progressing glaciation and helping to stabilize the climate system on its way to the late Cenozoic “icehouse” world.

The nature of the carbon cycle perturbation could be better constrained if we knew the evolution of atmospheric CO2 across the MMCT. Understanding the role of the carbon cycle in this cooling step is a key to assess Earth-system sensitivity to CO 2 forcing and the long-term stability of the Antarctic ice sheet under rising CO2 concentrations. However, most proxy records for the history of pCO2 across the MMCT are incomplete or at low resolution, thus prohibiting resolution of the CM events (Pagani et al., 1999;Kürschner et al., 2008;Foster et al., 2012;Ji et al., 2018;Sosdian et al., 2018;Super et al., 2018) and making it difficult to identify the mechanism responsible for the major step into the "icehouse" world. The only highresolution pCO2 reconstruction across CM6 is based on alkenone δ 13 C data from a Miocene outcrop on Malta (Badger et al., 2013), but it does not reveal a significant change.
Therefore, to better understand atmospheric CO2 evolution over the period of EAIS growth, we generated a high-resolution reconstruction based on δ 11 B measurements of fossil planktonic foraminifers from the South Atlantic (Figs. 1, 2). The boron isotopic composition of biogenic carbonates is a reliable recorder of ambient pH, which in turn is closely linked to atmospheric pCO2 in oligotrophic surface waters. Our record fully captures the carbon isotope excursions CM5b and CM6, which include the major expansion of the EAIS.

Sampling strategy
This research used samples from Ocean Drilling Program (ODP) Site 1092, located at 46.41 S and 7.08 E (Fig. 1) in a water depth of 1973 m. The interval studied (~178-184 mcd, see Table S1) consists of nannofossil ooze with excellent carbonate preservation, as shown by SEM-imaged undissolved coccoliths, which are susceptible to corrosion (Fig. 2). The shown example is from ~13.8 Ma, where pH was low, but shell preservation is similarly good throughout the entire record.
Approximately 200 specimens of the cold-water dwelling planktonic foraminifer G. bulloides were picked for δ 11 B analysis from 35 processed 10 cm 3 sediment samples. In addition, about 5 to 8 specimens of the benthic foraminifer Cibicidoides wuellerstorfi were picked from 14 intervals to reconstruct deep-sea pH (Table S1). Today, surface-water pCO2 at Site 1092 is close to equilibrium with the atmosphere (ΔpCO2sea-air< -30 µatm (Takahashi et al., 1993(Takahashi et al., , 2009). Plankton tow data revealed that G. bulloides generally lives within the upper 100 m of the water column (Mortyn and Charles, 2003). Based on DIC (dissolved inorganic carbon) and TA (total alkalinity) from seasonal TCO2+TALK (Goyet et al., 2000), and temperature and salinity from WOA13 data sets, modern gradients in pH and pCO2 are small at this site within the upper 200 m of the water column. During the austral spring season, when the shell flux of G. bulloides is highest (Jonkers and Kučera, 2015;Raitzsch et al., 2018), the differences in pH and pCO2 within this depth range are -0.03 and 17 µatm, respectively.
Hence, we expect G. bulloides to be a reliable recorder of subsurface pH and pCO2.

Age model
The age model used in this study is adopted from Kuhnert et al. (2009), but was revised to bring it in line with the astronomically tuned stable isotope record of the reference IODP Site U1338 (Holbourn et al., 2014) (Fig. 3, Table S2). In addition, chron C5ADn (top) identified at 182.51 mcd and originally dated to 14.178 Ma (Censarek and Gersonde, 2002) was tied to 14.081 Ma from the age model of Site U1338 (Holbourn et al., 2014. For direct comparison of our pCO2 record with the boron-based reconstructions available from the literature, the age models of ODP 761 (Holbourn et al., 2004, Foster et al., 2012, Sosdian et al., 2018 and the Blue Clay Formation of Ras-il-Pellegrin on Malta (Abels et al., 2005;Badger et al., 2013) were also re-tuned to match the reference curve U1338 (Table S2). Based on these revised age models, the δ 18 O and δ 13 C profiles of Sites 1092, ODP 791 and the Blue Clay Formation generally show a good agreement with those of Site U1338 (Fig. 3).

Boron isotope analysis
Foraminifer shells were cleaned following the protocol of Barker et al. (2003). Trace metal and boron isotope (δ 11 B) were measured following Raitzsch et al. (2018). Briefly, the cleaned samples were dissolved in 60 µL of 1 N HNO3 and microdistilled on a hotplate to separate boron from the carbonate matrix. The microdistillation method has been proven to yield a B recovery of ~100 %, a low procedural blank, and accurate results, even at low B concentrations (Gaillardet et al., 2001;Wang et al., 2010;Misra et al., 2014;Raitzsch et al., 2018). The distillate containing only boron was diluted with 2 % HNO 3 and analyzed for isotopes in triplicate using a Nu Plasma II multi-collector ICPMS at AWI (Bremerhaven, Germany) that is equipped with a customized detector array of 16 Faraday cups and 6 secondary electron multipliers (SEM), also termed ion counters (IC). 11 B and 10 B were collected in IC5 and IC0, respectively, at a boron concentration of ~3 ppb. As three of the SEMs were later replaced by Daly detectors, 11 B and 10 B of C. wuellerstorfi samples were collected in D5 and D0, respectively, at a boron concentration of ~2 ppb. 11 B/ 10 B was standardized against concentration-matched NBS 951 using the standard-sample-standard bracketing technique, and frequent analysis of control standard AE121 with an isotopic composition similar to that of foraminifers was monitored to ensure accuracy of measurement. Measurement uncertainties are reported as 2 standard deviations (2V) derived from triplicate measurements or as ±0.30 ‰ (for SEM-analyzed samples) and ±0.25 ‰ (for Daly-analyzed samples) determined from the long-term reproducibility (2V) of the control standard, whichever is larger (see Table S1).

Reconstruction of pH and pCO2 from δ 11 B
Planktonic δ 11 B was converted to sea-surface pH (Fig. S1) using the δ 11 Bforam/δ 11 Bborate relationship for G. bulloides (Raitzsch et al., 2018), where δ 11 Bborate is determined by the dissociation constant pKB and δ 11 Bseawater. pKB is dependent on salinity, temperature and pressure. Mg/Ca was taken from Kuhnert et al. (2009) for temperature estimates, and pressure was chosen from 50 m water depth, the assumed average calcification depth of G. bulloides. To account for seawater Mg/Ca, which was approximately 3.2 mol/mol ~14 Ma ago (Stanley and Hardie, 1999;Horita et al., 2002), Mg/Ca temperatures were corrected by applying a correction factor of 0.825 to the pre-exponential constant of the Mg/Ca-temperature equation of Mashiotta et al. (1999). This correction factor was previously determined for G. sacculifer by Evans and Müller (2012), but we assume that the effect of seawater Mg/Ca on shell Mg/Ca is similar for G. bulloides, since both have similar Mg/Ca ratios. However, the difference between Mg/Casw-corrected and non-corrected Mg/Ca absolute temperatures is in the order of 1.7 °C, but the relative changes are nearly identical.
This minimum gradient is required to keep the upper ocean density difference across the Subantarctic Front to enable the formation of Antarctic Intermediate Water (Kuhnert et al., 2009), which existed as Southern Component Intermediate Water (SCIW) since at least ~16 Ma  and spread into all oceans adjacent to the Southern Ocean (Wright et al., 1992). However, since the computed absolute salinity values are unrealistically low, an offset of 17.2 (psu) was added to the entire salinity record to achieve post-glaciation values similar to today. The reconstructed relative salinity change across the MMCT at Site 1092 is slightly more than 1 (psu), which is equaivalent to the salinity gradient across the Subantarctic front today. Deep-sea pH from C. wuellerstorfi δ 11 Bforam shown in Fig. S1 was calculated using a δ 11 Bforam:δ 11 Bborate relationship of 1:1 (Rae et al., 2011), deep-sea temperatures derived from Mg/Ca (data courtesy of H. Kuhnert) using the species-specific calibration of Raitzsch et al. (2008), a salinity of 34 (psu), and a paleo-water depth of 1794 m.
The most critical parameter for correct conversion to pH is δ 11 Bseawater, which is fortunately well constrained for the middle Miocene. We used a δ 11 Bseawater value of 37.80 ‰, which is the mean value from different independent studies and in close agreement with each other to within 0.1 ‰ (1V) (Pearson and Palmer, 2000;Foster et al., 2012;Raitzsch and Hönisch, 2013;Greenopet al., 2017), and also comparable to the 38.5 ‰ modeled by Lemarchand et al. (2000).
To calculate pCO2 from sea-surface pH, a second carbonate system parameter needs to be constrained. For this study, we used a sea-surface TA of 2000 µmol/kg, which is in line with various carbon cycle model results (Tyrrell and Zeebe, 2004;Ridgwell, 2005;Caves et al., 2016;Sosdian et al., 2018, Zeebe andTyrrell, 2019), and what we applied to the entire record.
However, to verify the validity of our approach we tested the effect of varying TA on pCO2 calculations (see section 3.2).
To assess whether a pH effect on Mg/Ca-derived temperatures and, in turn, pH-corrected temperature influences estimated pCO2 as shown by Gray and Evans (2019), we used a slightly modifed version of their 'MgCaRB' package written in 'R'.
Briefly, this program corrects Mg/Ca for the effects of salinity and pH and the resulting change in the temperature effect on pKB and hence on pH, until the change in both parameters between two iterations fall below certain threshold values. The original R code from Gray and Evans (2019) uses sediment age and associated error to determine past salinity and its uncertainty, based on the secular seawater salinity model from Spratt and Lisiecki (2016). Hence, the model is limited to the last ~800 ka. To adopt the code for our Miocene study, we slightly modified 'MgCaRB.d11B.R' from Gray and Evans (2019) by allowing for directly use estimated salinities and associated uncertainties. In addition, the seawater B isotopic composition was changed from 39.61 ‰ to the Miocene value of 37.80 ‰ (see code in the supplement). While the Miocene relative salinity variations can be constrained from changes in benthic foraminiferal δ 18 O, the absolute salinity levels are difficult to determine. Hence this approach only allows for comparing relative differences between the conventionally and the iteratively corrected pCO2 estimates. For this reason, the temperature and pCO2 values of the two approaches were brought in coincidence at the start of the record (Fig. S2).
Uncertainties of pCO2 estimates were fully propagated from individual uncertainties in pH (converted from 2V measurement uncertainty of δ 11 B), TA (±100 µmol/kg), temperature (±1 °C), and salinity (±1 psu) using the 'seacarb' package (Lavigne et al., 2011) programmed in 'R'. The applied temperature uncertainty is similar to the mean difference of ~2 °C between the non-corrected and pH-corrected Mg/Ca temperatures (Fig. S2), while the TA uncertainty is similar to the mean standard deviation of ±~130 µmol/kg between different TA models and clearly encompassing the standard deviation of ±~50 µmol/kg for each model across the MMCT (Sosdian et al., 2018). The applied salinity uncertainty encompasses the potential salinity change across the MMCT. We did not apply uncertainties for Mg/Casw and δ 11 Bsw, since both are systematic errors and would just shift the pCO2 record in either direction.

Carbon dioxide record
Our boron-based record indicates a steady pCO2 increase from ~390 to 520 ± 60 µatm between 14.25 and 14.08 Ma and a subsequent decrease to ~360 µatm until 13.87 Ma, culminating after the onset of the major ice-sheet expansion (Fig. 4E). It should be noted that this ~160 µatm drop in pCO2 occurred at a time of northward shifting Southern Ocean fronts, which is only related to a step-like drop in sea-surface temperature (SST) and salinity around 14 Ma at this site ( Fig. 4C) (Kuhnert et al., 2009). At 13.82 Ma, when SST had already reached a lower stable level and after the inception of the ~1 ‰ rise in benthic δ 18 O (Fig. 4A), pCO2 increased rapidly by ~120 µatm and displayed high-amplitude variations of more than 50 µatm until 13.57 Ma. This transient rise in pCO2 ended with a decrease to ~360 µatm and much reduced variability after 13.53 Ma (Fig. 4E). Sea-surface pH values within the entire record range from 7.87 to 8.07± 0.05 (Figs. 4D and S1), while deep-sea pH varies between 7.68 and 7.82 ± 0.07 over the same interval. Interestingly, the pH offset between the surface and deep ocean is nearly identical to the modern gradient and does not change substantially within our Miocene record (Fig. S1).
It is possible that the overall absolute pCO2 values and the glacial rise after 13.82 Ma are overestimated by up to ~80 µatm and ~30 µatm, respectively, since pH might affect planktonic foraminiferal Mg/Ca and hence the sea-surface temperature estimations (Fig. S2) (Gray and Evans, 2019). Since differences in relative pCO2 changes are not substantial and absolute values are within uncertainties between both reconstructions, confirming the overall robustness of our pCO2 record, we report the uncorrected data here.
The pCO2 variations within the studied time interval show a remarkable agreement with variations in δ 13 C that correspond to CM5b and CM6 (Fig. 5). Accordingly, the maxima at ~14.1 and 13.7 Ma as well as the minima at ~13.9 and 13.5 Ma indicate that atmospheric CO2 levels were paced by the 400 k.y. eccentricity cycle, as also demonstrated by evolutive harmonic and power spectral analyses (Fig. 6).

Sensitivity tests
Given that we propose an eccentricity-modulated change in TA input to the ocean, which might have affected our pCO2 reconstruction, we carried out several sensitivity tests. The first tests demonstrate that either constant salinity, a ~1 unit change in salinity, or TA co-varying with salinity using a modern TA-S relationship, has no discernible effect on pCO2 estimates (Fig. 7A). As ODP Site 1092 was apparently influenced by a northward migration of the Southern Ocean fronts, we tested the potential impact on calculated pCO2 by varying TA with temperature. Even at an unrealistically high total change of 400 µmol/kg, calculated relative pCO2 changes are very similar after the EAIS expansion at ~13.9 Ma (Fig. 7B).
Given the modern sea-surface TA gradient across the polar frontal system of ~100 µmol/kg (higher in the south), we conclude that the increasing influence of southern-sourced surface waters on pCO2 estimates at our study site is comparatively small. When in our simulations TA was varied with pH or eccentricity (both on short and long periodicity) by ±100 µmol/kg, pCO2 reconstructions are almost identical among each other (Figs. 7C-7H). On the other hand, at TA variations of ±400 µmol/kg, the pCO2 estimates may differ substantially between the different scenarios, but the general shape of the record remains more or less the same. Given that total variations in TA of 800 µmol/kg are very unlikely, we suggest that periodic changes in monsoon-driven TA do not affect our pCO2 estimates significantly.

Comparison with other pCO2 records
Our pCO2 data broadly agree with reconstructed long-term trends based on planktonic foraminiferal δ 11 B (Foster et al., 2012) and paleosol δ 13 C (Ji et al., 2018), despite a few differing values before CM5b and during CM6 (Fig. 3). The boron isotope 6 160 165 170 175 180 data from the Mediterranean Sea by Badger et al. (2013) also show decreasing pCO2 values after 13.8 Ma, but no rebound after 13.75 Ma and no rise after 13.85 Ma, although these discrepancies may be due to the different sampling resolutions. In contrast, the pCO2 reconstructions from high-resolution alkenone δ 13 C data from Badger et al. (2013) are similar to the estimates of Super et al. (2018), which are considerably lower than the boron-based reconstructions and do not exhibit significant changes CM6. This might be related to the insensitivity of the alkenone proxy at low to moderate pCO2 levels (Badger et al., 2019). It is worth noting that our pCO2 record reveals a striking match with the box model output presented in Ma et al. (2011), suggesting a linear response of weathering-controlled nutrient input to changes in orbital parameters with attendant variations in ocean carbon reservoir during the middle Miocene.

The role of eccentricity in the carbon cycle
The pCO2 increase preceding and following the EAIS expansion closely tracks the global δ 13 C history of CM5b and CM6 ( Fig. 4), emphasizing a fundamental link between changes in the ocean carbon cycle and atmospheric carbon dioxide. The CM peaks correspond to minima in the 400 k.y. eccentricity cycle, suggesting that δ 13 C variations were related to changes in monsoon and weathering intensity (Holbourn et al., 2007;Ma et al., 2011). Monsoon intensity is paced by solar radiation variations caused by precessional cycles but their amplitude variations are modified by eccentricity. More precisely, higher eccentricity results in larger precession amplitudes and hence in larger wet/dry variations in the tropics (e.g., Wang, 2009), stronger physical and chemical weathering, and ultimately in an increased input of particulate organic material, dissolved inorganic carbon, alkalinity and nutrients to the oceans (e.g., Clift and Plumb, 2008;Ma et al., 2011;Wan et al., 2009). The potential influence of alkalinity variations on our pCO2 record is shown in Fig. 7 and evaluated in section 3.2.
The increased nutrient supply enhanced the primary production and organic carbon burial, which in turn lowered pCO2. At the same time, the burial of shallow and total calcium carbonate (CaCO3) increased, while CaCO3 burial in the deep ocean decreased (Holbourn et al., 2007). However, as the net burial of CaCO3 (enriched in 13 C) in relation to Corg (depleted in 13 C) increased at high eccentricity, δ 13 CDIC values decreased, which is in line with the global isotope signature (Holbourn et al., 2007;Ma et al., 2011).
The increase in carbonate burial in shallow seas should have removed alkalinity from seawater, and thus lowered pH and released CO2 to the atmosphere (e.g., Zeebe and Wolf-Gladrow, 2001) during 100 k.y. eccentricity maxima, but the resolution of our record is too low to detect this. On longer timescales (400 k.y. cycles), however, increased alkalinity input from rivers, dissolution of deep-sea carbonates, and the enhanced burial of Corg in tropical regions during eccentricity maxima might have contributed to the long-term decrease in pCO2. This agrees well with the box model output from Ma et al. (2011), suggesting that during the long-eccentricity maxima pCO2 decreased when the monsoon was most intense and the organic carbon burial and river fluxes were high. Conversely, weathering and nutrient supply to the ocean in low latitudes decreased when eccentricity was low, resulting in a net decrease in the CaCO3-to-Corg burial ratio and hence in CO2 release and higher δ 13 C.

Enhanced glacial deep-water ventilation
The MMCT bears similarities to the Eocene-Oligocene transition (EOT) ~34 Ma ago, when Antarctic ice-expansion was also accompanied by a large positive δ 13 C excursion, in terms of both magnitude and duration (Coxall et al., 2005). The δ 13 C excursion following Antarctic glaciation was accompanied by a transient pCO2 rise of more than 300 µatm (Pearson et al., 2009;Heureux and Rickaby, 2015), which bears similarity to the pCO2 rise during the mid-Miocene EAIS expansion, although the amplitude was much larger. A modeling study suggested that a sea-level-fall-induced shelf-basin carbonate burial fractionation and a temporal enhancement of deep-water formation (by 150 %) in the Southern Ocean were sufficient to produce both the positive δ 13 C excursion and a transient pCO2 rise by ~80 µatm in the EOT (McKay et al., 2016).
Similarly, ocean circulation may have intensified during the Miocene EAIS expansion, as the latitudinal temperature gradient increased and atmospheric circulation strengthened. This is supported by benthic foraminifers, Mn/Ca and XRF data from Southeast Pacific sites, which indicate improved deep-water ventilation and carbonate preservation after 13.9 Ma, particularly during colder climate phases (Holbourn et al., 2013). Improved deep-water ventilation may also have led to enhanced advection of silica-rich waters toward low latitudes, culminating in an increased upwelling and diatom productivity between 14.04 and 13.96 Ma and between 13.84 and 13.76 Ma in the equatorial East Pacific, as indicated by massive opal accumulation found at IODP Site U1338 (Holbourn et al., 2014). We speculate that increased upwelling may have resulted in CO2 outgassing, if the supply exceeded consumption by primary productivity.

Conclusions
We conclude that pCO2 variations across the MMCT were paced by 400 k.y. eccentricity cycles, where pCO2 decreased at high eccentricity and rose when eccentricity was low (Fig. 4). At high eccentricity, the global monsoon and hence weathering intensity increased (e.g., Wang, 2009), causing an increased input of dissolved inorganic carbon, alkalinity and nutrients to the oceans. The resulting increase in the CaCO3-to-Corg burial ratio lowered pCO2 and decreased δ 13 CDIC. The pCO2 decrease between 14.08 and 13.87 Ma might have facilitated the initiation of EAIS expansion starting at ~13.9 Ma.
Conversely, the CO2 rise after the onset of EAIS expansion, possibly caused by decreased weathering fluxes at low eccentricity and hence net decrease in the CaCO3-to-Corg burial ratio, could have acted as a negative feedback on the progressing Antarctic glaciation. In this way, the radiative forcing due to the temporary pCO2 rise may have helped to stabilize the climate system on its way to the late Cenozoic "icehouse" world. Our results also highlight the need for more high-resolution pCO2 records across the middle Miocene climate transition.

Data availability
The boron isotope data collected for this study are available from Table S1, the tie points used for the revised age model are listed in Table S2, and the modified R code of 'MgCaRB' (Gray and Evans, 2019) is available from S1 in the supplement.    were re-calculated using a TA of 2000 µmol/kg and a δ 11 Bsw value of 37.80 ‰ for consistency with our data. In addition, the age models from these studies were also revised for direct comparison with our record (Fig. 3). All other records are based on their original data and age models.