Revisiting the Middle Eocene Climatic Optimum “Carbon Cycle Conundrum” With New Estimates of Atmospheric pCO2 From Boron Isotopes

The Middle Eocene Climatic Optimum (MECO) was a gradual warming event and carbon cycle perturbation that occurred between 40.5 and 40.1 Ma. A number of characteristics, including greater‐than‐expected deep‐sea carbonate dissolution, a lack of globally coherent negative δ13C excursion in marine carbonates, a duration longer than the characteristic timescale of carbon cycle recovery, and the absence of a clear trigger mechanism, challenge our current understanding of the Earth system and its regulatory feedbacks. This makes the MECO one of the most enigmatic events in the Cenozoic, dubbed a middle Eocene “carbon cycle conundrum.” Here we use boron isotopes in planktic foraminifera to better constrain pCO2 changes over the event. Over the MECO itself, we find that pCO2 rose by only 0.55–0.75 doublings, thus requiring a much more modest carbon injection than previously indicated by the alkenone δ13C‐pCO2 proxy. In addition, this rise in pCO2 was focused around the peak of the 400 kyr warming trend. Before this, considerable global carbonate δ18O change was asynchronous with any coherent ocean pH (and hence pCO2) excursion. This finding suggests that middle Eocene climate (and perhaps a nascent cryosphere) was highly sensitive to small changes in radiative forcing.


Introduction
The Middle Eocene Climatic Optimum (MECO;~40.1-40.5 Ma) was a global warming event during which marine bulk and benthic carbonate δ 18 O values steadily declined by roughly 1‰ in over 400 kyr, usually interpreted as 3-6°C of global temperature rise (Bohaty et al., 2009;Bohaty & Zachos, 2003). While sometimes referred to as a hyperthermal event (e.g., Arreguín-Rodríguez et al., 2016;Pomar et al., 2017), the MECO differs fundamentally from the true hyperthermal events earlier in the Paleogene, such as the Paleocene-Eocene Thermal Maximum (PETM,~56 Ma) or Eocene Thermal Maximum 2 (ETM2;~53.7 Ma). First, these hyperthermal events were rapid in their onset (<10 kyr; Kirtland Turner et al., 2017) and were followed by a gradual return to roughly pre-event temperatures as the silicate weathering feedback drew down atmospheric CO 2 (e.g., Kelly et al., 2005;Penman, 2016;Penman et al., 2016;Zachos et al., 2005), whereas the MECO saw a steady decline in carbonate δ 18 O values over several 100 kyr (interpreted as gradual warming), followed by a rapid return to pre-event conditions (Bohaty et al., 2009). A second important distinction is that the earlier hyperthermal events are clearly marked by sharp drops in δ 13 C that are expressed globally in both marine CaCO 3 and terrestrial records (e.g., Abels et al., 2012;Kennett & Stott, 1991;Koch et al., 1992;Westerhold, Röhl, Donner, Frederichs, et al., 2018;Westerhold, Röhl, Donner, & Zachos, 2018;Zachos et al., 2005). These carbon isotope excursions (or CIEs) are thought to be manifestations of rapid injections of isotopically light C into the Earth system, perhaps from volcanic sills intruding into organic rich sediments (Gutjahr et al., 2017;Storey et al., 2007;Svensen et al., 2004) and/or the release of seafloor methane clathrates (Dickens et al., 1995). A small negative δ 13 C excursion is observed at the peak of the 400 kyr MECO warming trend in most marine records but is not uniformly expressed (e.g., at ODP Site 1263; Bohaty et al., 2009) and appears in some cases to lag behind minimum δ 18 O values (e.g., at ODP Site 702; Bohaty et al., 2009). Furthermore, during the first few 100 kyr of warming, patterns of carbon isotope change are inconsistent between sites: global bulk carbonate δ 13 C values display inverse trends in each hemisphere, with progressively lighter δ 13 C during MECO warming in northern latitudes, minimal δ 13 C change in the tropics, and progressively heavier δ 13 C during warming toward the southern high latitudes (Bohaty et al., 2009). Thus, unlike the hyperthermal events, the MECO is difficult to attribute to a simple, rapid injection of exogenous isotopically light carbon.
One characteristic that the MECO perhaps unexpectedly shares with the hyperthermal events is extensive dissolution of deep-sea sediments, representing a shoaling of the carbonate compensation depth (CCD; Bohaty et al., 2009;Bohaty & Zachos, 2003). During the hyperthermal events, this carbonate dissolution occurred in response to CO 2 being injected into the ocean-atmosphere system faster than the capacity of silicate weathering on land to draw it down and input alkalinity to the ocean (see, e.g., Berner et al., 1983;Colbourn et al., 2015;Kirtland Turner et al., 2017). MECO warming, however, was nearly an order of magnitude slower than the PETM, and as such, chemical weathering feedbacks should have kept pace with warming (Colbourn et al., 2015;Sluijs et al., 2013). This would be expected to have driven a deepening of the CCD instead of the~1 km shoaling observed. Moreover, at the MECO, there is little signal of a "carbonate overshoot" following the event (Bohaty et al., 2009;Sluijs et al., 2013), as is seen after the PETM . Since these carbonate overshoots are thought to reflect the Earth's silicate weathering feedback drawing down injected CO 2 , the absence of any such event after the MECO might suggest either (a) a relatively small addition of CO 2 and/or (b) a weakened silicate weathering response to CO 2 rise during this interval. The latter scenario is favored by a recent study by van der Ploeg et al. (2018), who observed a  Gradstein et al., 2012). In each case, excursions toward higher benthic foraminiferal δ 13 C coincide with instances where δ 18 O changes diverge from the underlying trend of cooling and ice growth in the middle-late Eocene. "E-OT" as marked in the gray band denotes the Eocene-Oligocene transition. reduction in the 187 Os/ 188 Os composition of seawater over the MECO. While this observation could be consistent with (a) a rise in hydrothermal activity, (b) an increase in weathering of mafic silicate rocks, or (c) a decrease in weathering of felsic radiogenic continental rocks, a reduction in continental weathering is more easily reconciled with deep-sea carbonate dissolution. Consequently, van der Ploeg et al. hypothesize that after tens of millions of years of intense weathering under greenhouse Eocene climates, Earth's silicate weathering feedback was no longer strong enough to buffer a carbon cycle imbalance at the MECO (van der Ploeg et al., 2018).
While a weakened silicate weathering feedback is an intriguing hypothesis, it does pose its own questions. First, the idea that the baseline middle Eocene earth surface was no longer as readily weatherable is based on Li isotopes measured in foraminifera (Misra & Froelich, 2012). However, even leaving potential non-seawater controls on foraminiferal δ 7 Li aside (Roberts et al., 2018;, there is more than one possible solution that could have resulted in light seawater δ 7 Li, including sustained high rates of regolith chemical weathering and soil production (Vigier & Goddéris, 2015). Second, one of the most pronounced features of the MECO is a rapid (10s of kyr) drop in temperatures after the peak of the warming trend, which would indicate some unexplained reactivation of the silicate weathering feedback after million years of hypothesized gradual decline.
Another outstanding question surrounding the MECO is that of causality. The long, drawn-out nature of the MECO warming might suggest a sustained CO 2 release, and hence perhaps a volcanic source, by analogy to the several 100 kyr-long late Maastrichtian warming event (LMWE) thought to have been caused by Deccan outgassing (e.g., Barnet et al., 2018;Henehan, Hull et al., 2016). A number of potential sources of volcanic CO 2 at the MECO exist (summarized in van der Ploeg et al., 2018), but the synchronicity of any one volcanic degassing event with MECO warming is yet to be demonstrated. Besides this, the coincidence of the MECO with a 2.4 Myr very long eccentricity cycle minimum (Westerhold & Röhl, 2013) might implicate an orbital trigger instead. These long eccentricity 2.4 Myr cycles have recently been shown to coincide with numerous excursions in δ 13 C and δ 18 O during the early Cenozoic (Barnet et al., 2019;Kocken et al., 2019), and carbon cycle models have been able to reproduce significant changes in the Earth's carbon cycle due to oscillating carbon reservoirs on these timescales, inducing sizeable changes in pCO 2 and the CCD (Kocken et al., 2019). Observed cyclical changes in the Pacific CCD in the Eocene (Lyle et al., 2005;Pälike et al., 2012) appear to have a roughly 2.4 Myr periodicity (Kocken et al., 2019), and the MECO follows the most pronounced and extreme of these periodic carbonate accumulation events (CAEs) (Lyle et al., 2005). It is possible that a combination of volcanic outgassing and a favorable orbital configuration may have been key to the amplification of the MECO as a major carbon cycle perturbation, as has been suggested for the PETM and LMWE (Barnet et al., 2019). Alternatively, it is also conceivable that orbital configuration alone could have independently induced carbon cycle changes (as with the Early-Late Paleocene Event; Barnet et al., 2019), perhaps through changing rainfall patterns and hence silicate weathering intensity (Westerhold & Röhl, 2013). Indeed, repeating patterns in carbonate δ 13 C and δ 18 O values ( Figure 1) might favor some such internally modulated, pseudo-cyclic carbon cycle imbalance rather than a stochastic carbon injection from a volcanic source.
Constraining atmospheric pCO 2 across the MECO is central to understanding both the causal drivers of the MECO and the Earth system's enigmatic response to the event. The only published pCO 2 reconstruction over the MECO to date is that of Bijl et al. (2010), derived from alkenone δ 13 C. While undoubtedly pioneering, the record poses even more questions of its own, in that taken at face value it indicates several high-amplitude fluctuations on the order of several thousand μatm of CO 2 over the MECO, asynchronous with a near-linear decline in δ 18 O. This would appear incongruous in the context of background pCO 2 levels in the middle Eocene reconstructed from boron isotopes (Anagnostou et al., 2016) and would require extremely large masses of carbon to have been injected into-and subsequently removed from-the Earth system on geologically short timescales with no discernible influence on carbonate δ 13 C values. Here, we revisit this issue using the boron isotope (δ 11 B) pH proxy (see, e.g., Foster & Rae, 2016;Hemming & Hanson, 1992) to Figure 2. Locations of our study sites (colored markers) on a continental reconstruction from~40 Ma from www.odsn.de (Hay et al., 1999). The location of ODP Site 1172 (where alkenone-based pCO 2 reconstructions from Bijl et al., 2010, were derived) is also marked as a white star.
generate new, independent estimates of ocean pH and pCO 2 change at this time. We then explore the significance of our record in constraining causal mechanisms, global weathering feedbacks, and climate sensitivity over this enigmatic event.

Site and Species Selection
We analyzed single-species mixed-layer foraminifera primarily at the equatorial Pacific ODP Site 865 (1,300-1,500 m paleo-depth; Bralower et al., 1995), which has been the subject of previous boron isotope investigations Pearson & Palmer, 1999, 2000. To ascertain to what extent our resultant pCO 2 record was globally representative, we compared these data with lower-resolution measurements from three other deep-sea drill sites, spanning a wide range of oceanographic settings ( Figure 2): ODP Sites 1260 (equatorial Atlantic;~2,500-3,200 m paleo-depth; Sexton et al., 2006;Shipboard Scientific Party, 2004a), 1263 (midlatitude South Atlantic;~2,000 m paleo-depth; Shipboard Scientific Party, 2004b), and 702 (Southern Ocean;~2,250 m paleo-depth; Katz & Miller, 1991;Shipboard Scientific Party, 1988). High-resolution bulk carbonate or benthic foraminiferal δ 13 C and δ 18 O stratigraphies are available for each site (Bohaty et al., 2009;Boscolo Galazzo et al., 2014;Edgar et al. (subm.); Edgar et al., 2010;Sexton et al., 2006b), allowing cross-correlation. All sites are open-ocean in setting and remain above the local CCD across the event, although at ODP Site 1260 multiple clay layers are evident within the MECO (Edgar et al., 2010), indicating that this site was close to the local lysocline. It is likely, therefore, that our sampling at ODP Site 1260 misses the peak of the MECO event. Calcareous nannofossil and foraminiferal communities from ODP Site 865 indicate that sediments were deposited under oligotrophic open-ocean conditions (Bralower et al., 1995), suggesting that [CO 2 ] aq reconstructed at these sites was likely in (near) equilibrium with the atmosphere with respect to CO 2 . Abundant siliceous microfossils suggest more productive surface ocean conditions at ODP Sites 1260 and 702, while productivity at ODP Site 1263 appears to vary temporally over the event (Boscolo Galazzo et al., 2015).
At each site, monospecific separates (136-343 tests, 1.8-4.3 mg CaCO 3 ) of planktic foraminifera were picked for δ 11 B analyses from the 250-300 μm sieve size range to minimize any ontogenetic or metabolic effects. We selected the shallowest-dwelling foraminiferal species possible at each site, given the sample size requirements for analysis. Specifically, these were Globigerinatheka index (at ODP Site 702); Acarinina praetopilensis (ODP Site 865); Morozovelloides crassatus (ODP Site 1260); and Globigerinatheka kugleri (ODP Site 1263). Foraminifera at all sites are recrystallized, but the degree of diagenetic recrystallization varies between sites (see Sr/Ca data, supporting information). The most pervasive alteration is seen at ODP Site 865, but even extensive recrystallization in Eocene sediments from this site appears to have had little influence on fossil foraminiferal δ 11 B values .

Age Models
All ages are reported on the Geological Timescale 2012 (Gradstein et al., 2012). Sites were aligned using a combination of stable isotopic, magneto-stratigraphic, and biostratigraphic tie points. An internally consistent age model for all sites was constructed based on ODP Site 702, which has the most complete magnetostratigraphy (Clement & Hailwood, 1991) and most highly resolved bulk stable isotope stratigraphy in this study (Bohaty et al., 2009). All stable isotope tie points were calibrated at this site based on linear interpolation between magnetochron boundaries. Although magnetostratigraphy at all sites was developed with biostratigraphic controls in mind, we used bioevent tie points for age control only outside of the focal interval, at depths where no stable isotope stratigraphy or reliable magnetic data were available for correlation. This is because bioevents are difficult to employ for correlation across such a large latitudinal range used in this study, given likely faunal differences between sites at any one time. Immediately surrounding and within the MECO itself, stable isotope tie points were most useful for correlation. Our age models therefore represent a compromise between number of robust (easy-to-pick) tie points and alignment of major stable isotope features at all sites. A full list of the stable isotope tie points used in this study and the age assigned to each are shown in Table S1 and are plotted in the depth domain in Figure S1.

Analytical Methods
Boron isotope analyses were carried out on a Thermo Neptune Multi-Collector ICP-MS at the University of Southampton following Foster (2008). External reproducibility is estimated from 11 B signal intensity (primarily dictated by sample size) using a relationship between intensity and long-term reproducibility at the University of Southampton from Greenop et al. (2017). Sample cleaning and preparation follow established protocols (see Anagnostou et al., 2016;Henehan et al., 2016, and references therein). For each sample, a subsample of cleaned, crushed foraminiferal carbonate was taken for δ 13 C and δ 18 O analyses, while an aliquot of dissolved material was taken prior to boron column chemistry for trace and minor element analysis (B, Mg, Al, and Sr). Trace element analysis was carried out using a Thermo Element ICP-MS at the University of Southampton. Long-term reproducibility (2 SD) was ±<5% for Al/Ca ratios measurements, and ±<3% for Mg/Ca and Sr/Ca. Al/Ca ratios for samples analyzed here are all <105 μmol/mol, below an operational threshold of~140 μmol/mol commonly used to indicate clay contamination (Rae et al., 2011). Foraminiferal δ 13 C and δ 18 O analyses from ODP Sites 702, 865, 1260, and 1263 were made on a Thermo Scientific Delta V Advantage mass spectrometer coupled to a GasBench II at Cardiff University and are reported relative to the Vienna Pee Dee Belemnite (VPDB) standard with an external analytical precision (1σ) of 0.06‰ for δ 13 C and 0.07‰ for δ 18 O. Further samples from IODP Site U1408 were measured for δ 13 C and δ 18 O on a Thermo Kiel IV coupled to a Thermo MAT 253 mass spectrometer at Yale University, with typical external precision (1σ) of <0.03‰ for δ 13 C and <0.04‰ for δ 18 O based on replicate analyses of in-house (TS, MERC, CM, and PX) and international (NBS19) standards.
Vital effects in extinct species can be gauged by comparison of δ 11 B, δ 13 C, δ 18 O, and B/Ca ratios across a foraminiferal assemblage ( . Planktic foraminiferal stable oxygen (a) and carbon (b) isotope gradients with test size for species used in this study, measured at IODP Site 1408. Preliminary age assigned to this time-slice is between 40.0 and 40.5 Ma or planktic foraminiferal Zone E12. Thermocline dwelling, symbiont-barren species Subbotina linaperta and benthic foraminifera Nuttallides truempyi are shown for comparison.
using intra-specific body size δ 13 C and δ 18 O gradients (Birch et al., 2012;D'Hondt et al., 1994;Edgar et al., 2013;Ezard et al., 2015;Pearson et al., 1993;Sexton, Wilson, & Pearson, 2006b). To constrain vital effects in the species measured here, we made new δ 13 C and δ 18 O measurements of well-preserved "glassy" foraminifera from IODP Site U1408 where all of the species measured here coexisted in a single time-slice (Norris et al., 2014) (Figure 3). This new δ 13 C and δ 18 O data support a shallow mixed-layer habitat for A. praetopilensis and M. crassatus. However, our data support a somewhat deeper mixed-layer habitat for G. index and G. kugleri, as suggested elsewhere (Anagnostou et al., 2016;Edgar et al., 2013Edgar et al., , 2015Pearson et al., 1993Pearson et al., , 2001Pearson et al., , 2006Sexton, Wilson, & Pearson, 2006b). Globigerinathekids record cooler δ 18 O temperatures than M. crassatus and A. praetopilensis (Figure 3a), with G. index in particular indicating progressively deeper habitat in size fractions larger than those analyzed here, perhaps reflecting the addition of gametogenic crusts at depth (Pearson et al., 1993;Premoli Silva et al., 2006;Sexton, Wilson, & Pearson, 2006b). With regards δ 13 C, globigerinathekids at Site U1408 record lower values than M. crassatus or A. praetopilensis, indicating weaker photosynthetic activity in their microenvironment and/or a deeper mean habitat. While G. kugleri displays a strong positive δ 13 C gradient with increasing size (Figure 3b) consistent with symbiont-bearing foraminifera (Ezard et al., 2015;Pearson et al., 1993;Spero et al., 1991), G. index does not, suggesting that high latitude G. index had few or no photosynthetic symbionts. Thus, our focal species likely have varying degrees of influence from photosymbionts. By analogy with extant taxa, these species are expected to display variable offsets from aqueous δ 11 B borate (Henehan, Foster, et al., 2016). We therefore assign calibrations by analogy with modern taxa based on our δ 13 C and δ 18 O data and published observations (Anagnostou et al., 2016;Edgar et al., 2013;Pearson et al., 1993Pearson et al., , 2001Sexton, Wilson, & Pearson, 2006b). For A. praetopilensis and M. crassatus, we apply a calibration based on the low-latitude, mixed-layer dwelling, dinoflagellate-bearing Trilobatus sacculifer (Foster et al., 2012). For deeper mixed-layer G. kugleri, we apply a calibration based on the similarly deeper-dwelling symbiont-bearing species Orbulina universa (Henehan, Foster, et al., 2016). For G. index, which lacks any gradient in δ 13 C with size ( Figure 3b), we apply the symbiont-barren calibration for Globigerina bulloides (Martínez-Botí et al., 2015). In each case, calibration lines were adjusted to account for a middle Eocene seawater δ 11 B (δ 11 B sw ) of 38.7‰ using the approach of Greenop et al. (2019). This is necessary, as otherwise modern calibrations conflate the effects of pH and δ 11 B sw on deep time records (see Greenop et al., 2019, for more details).
We note also that reduced body size δ 13 C gradients in Acarinina during peak MECO at ODP Sites 1051 and 748 have been interpreted as symbiont "bleaching" (Edgar et al., 2013; see also Wade et al., 2008). Since this could modify the pH of the foraminiferal microenvironment and artificially produce or accentuate a surface ocean pH drop, we tested whether bleaching occurred at the MECO by examining body size δ 13 C gradients at tropical ODP Sites 865 and 1260, where temperature-driven "bleaching" might be most evident. However, we found no evidence for a reduction in either the δ 13 C gradient between symbiont-bearing and symbiont-barren taxa or the body size δ 13 C gradient in the species analyzed here ( Figure S2). This suggests that a loss of symbionts did not occur at these low-latitude sites and that any δ 11 B changes should therefore reflect changes in ambient seawater pH rather than changing vital effects.

pH and Temperature Calculations and Uncertainty Propagation
To calculate pH accurately from δ 11 B, a number of auxiliary parameters are required, including δ 11 B sw , salinity, and temperature. Using the approaches of Anagnostou et al. (2016) (viz., assuming reasonable bounds on surface ocean saturation and vertical DIC gradients) and a δ 11 B sw -corrected T. sacculifer calibration for their mixed-layer dwelling, symbiont-bearing early and middle Eocene foraminifera, we calculate absolute bounds of δ 11 B sw of 38.5-38.9‰ at this time. We assign a salinity of 34 ± 1 psu to all samples (note that this parameter has only a small effect on pH and pCO 2 estimates; e.g., Henehan et al., 2013). To estimate sea surface temperature, we used planktic foraminiferal Mg/Ca ratios, as they appear more robust in the face of diagenetic recrystallization than δ 18 O (Sexton et al., 2006a). However, modern Mg/Ca temperature calibrations cannot be directly applied to middle Eocene foraminifera, because of the different Mg/Ca sw at this time and the (sizeable) effect of pH on foraminiferal Mg/Ca (e.g., Gray & Evans, 2019). All our temperature calculations are therefore corrected for the nonlinear response of foraminiferal Mg/Ca to changing Mg/Ca sw following . Since pH calculations require temperature estimates, and vice versa, we iteratively correct both parameters, similar to methods given by Gray and Evans (2019).
Low Eocene seawater pH means that pH corrections on Mg/Ca-derived temperatures are sizeable ( Figure S3). To our knowledge, these are the first Mg/Ca sw and pH-corrected temperatures for this period and thus provide new constraints on sea surface temperatures across a latitudinal gradient in the middle Eocene.
For uncertainty propagation, we used Monte Carlo simulations in R (R Core Team, 2015). For pH estimates, 500 δ 11 B calcite values were randomly generated from within a normal distribution described by 2σ analytical reproducibility for each data point (see Table S1). Each value was then converted to δ 11 B borate according to each species' δ 11 B sw -specific calibration (see section 2.4). Major ion chemistry-specific equilibrium constants (from Hain et al., 2015) for each value were calculated according to preliminary Mg/Ca-derived temperatures (from 500 simulations of Mg/Ca within 3% analytical uncertainty, corrected for Mg/Ca sw ; Evans, Brierley, et al., 2016), but not pH (Evans, Wade, et al., 2016), and 500 randomly generated values of salinity (within 34 ± 1). pH was then calculated for each simulated value and used to derive the "excess Mg/Ca" due to low pH (according to the linear fit of Evans, Wade, et al., 2016). pH-corrected temperatures were then calculated for each Monte Carlo replicate and subsequently used to recalculate equilibrium constants and pH. This was iterated until simulated replicates converged on a unique solution to within 0.002 pH. Final temperature and pH estimates were calculated as the mean of the 500 Monte Carlo simulations, with 95% interquantile range of these replicates taken as uncertainty. In these calculations, we include uniformly distributed uncertainty on δ 11 B sw within the range of 38.5-38.9‰ (Anagnostou et al., 2016), when presenting true uncertainty on absolute values of pH. However, the long residence time of boron in the ocean (~10 to 20 Myr) means that δ 11 B sw would have remained constant across the time interval studied here, and so uncertainty on relative changes in pH is considerably smaller. To isolate the effect of δ 11 B sw on relative pH changes and provide the basis for LOSCAR simulations discussed later in section 2.7, we also calculated pH across a range of discrete values of δ 11 B sw . We do not attempt to propagate uncertainty in Eocene major ion chemistry or (small) salinity effects on foraminiferal Mg/Ca (Hönisch et al., 2013) since quantitative constraints over this interval are limited.

pCO 2 Calculations
To calculate pCO 2 from surface ocean pH, one requires one other carbonate system parameter (Zeebe & Wolf-Gladrow, 2001). Because we have no suitable proxy record to obtain this second parameter, we explore a number of possible scenarios using simplified assumptions to represent end-member pCO 2 possibilities. First, we calculate pCO 2 assuming a constant calcite saturation state (Ω calcite ) in the surface ocean at each site. This assumption represents a scenario in which alkalinity increases driven by elevated silicate weathering and carbonate dissolution could effectively keep pace with the pH decline calculated from δ 11 B. Specifically, in our Monte Carlo simulations, we allow starting Ω calcite to vary within a uniform distribution between 5.75 and 7.5 at Sites 865 and 1260; 5.25 and 7.25 at Site 1263; and 3.75 and 5.25 at Site 702, with differences reflecting each site's paleolatitude (following Anagnostou et al., 2016). For each Monte Carlo simulation, Ω calcite was then held constant through the MECO, and the carbonate system solved to calculate pCO 2 , according to the carbonate equilibrium constants and calcite solubility product calculated for each data point.
As a second end-member scenario, we performed pCO 2 calculations using the assumption that surface ocean total alkalinity (TAlk) remained constant throughout the event. By contrast with the constant Ω calcite scenario, this assumption implies that silicate weathering and carbonate compensation were too weak to significantly affect surface ocean carbonate chemistry during the MECO and that changes in pH reflect the addition of DIC alone. This assumption of constant TAlk would be more consistent with a postulated weakened silicate weathering feedback strength across the MECO (van der Ploeg et al., 2018). To obtain a pre-event surface ocean TAlk estimate, we took pre-event equilibrium values from simulations of the geochemical box model LOSCAR, matched to pre-event mean LOESS pH fit (see section 2.7). We then held this alkalinity constant at within ±150 μMol of LOSCAR's pre-MECO surface ocean TAlk estimate throughout the MECO.

Carbon Cycle Modeling MECO Acidification
Regrettably, the complex, spatially variable and as-yet poorly constrained array of potential influences on carbonate δ 13 C values during the MECO mean that reverse modeling the event (e.g., as done for the PETM by Gutjahr et al., 2017) is not yet feasible. However, because LOSCAR (Long-term Ocean Sediment Carbon Reservoir, v2.0.4; Zeebe, 2012) approximates some critical processes for interpreting pH recordsnamely, the changes in chemical weathering and carbonate burial one would normally associate with pCO 2 rise-we could use it to further explore feasible changes in pCO 2 over the MECO across a broad range of possible δ 11 B sw . More specifically, we leveraged LOSCAR to obtain estimates for pCO 2 from MECO carbon cycle model scenarios that were constrained only by δ 11 B-derived pH. To do this, we used LOSCAR's "PALEO" configuration (with a Tethys ocean and Paleogene circulation patterns and bathymetry) but with some modifications. First, we replaced the default corrections for [Mg] and [Ca] (Ben-Yaakov & Goldhaber, 1973) with updated equilibrium constants that better quantify ion pairing effects (Hain et al., 2015) by accounting for the activities of all ions in solution using the larger and more up-to-date MIAMI database (Millero & Pierrot, 1998), using a seawater [Mg] and [Ca] of 38 and 17 mM, respectively. Second, to more finely resolve deep water carbonate dissolution and compensation depth, we modeled one sediment level every 100 m, rather than the default 500 m (following Henehan, Hull, et al., 2016).
To obtain pre-MECO alkalinity estimates for any given δ 11 B sw scenario, the model was equilibrated by adjusting equilibrium pCO 2 in 10 Myr-long spin-up runs with a default strength silicate weathering feedback (N Si = 0.2) until we attained a fit between LOSCAR's low-latitude surface Pacific box and boron-derived pre-MECO surface ocean pH values (to within 0.002 pH units). This target pH for a given δ 11 B sw was calculated as the mean of the first 17 values of the mean LOESS fit prior to the pH drop at the peak MECO. For each value of δ 11 B sw /starting pH, we then iteratively added different masses of CO 2 to the atmosphere over the timescale of our peak MECO pH excursion (~100 kyr) until we produced a surface Pacific pH decrease that matched the minimum value in the δ 11 B-derived LOESS pH curve to within 0.002 pH units. This exercise was repeated across a broad range of δ 11 B sw (37.4-39.5‰), and in each case, ΔpCO 2 was calculated as the difference in pCO 2 from spin-up to peak MECO. CO 2 forcing change (ΔF) was then calculated from ΔpCO 2 , according to the relationship of Myhre et al. (1998). Our primary simulations used default parameterized climate sensitivity in LOSCAR of 3°C warming per doubling of atmospheric pCO 2 . Since calculated warming influences CO 2 solubility, DIC speciation, and carbonate solubility calculations in LOSCAR, we also repeated our experiments with a prescribed 2°C warming over the 100 kyr carbon release event to test for the influence of this parameter, with negligible effect on our CO 2 estimates (at δ 11 B sw = 38.7‰, ΔpCO 2 changes by <0.002 doublings).

Boron Isotope Records
Boron isotope measurements are shown in Figure 4. The main feature of our record at our main site, ODP Site 865, is a transient~1‰ δ 11 B decrease focused around minimum MECO δ 18 O values (~40.22-40.10 Ma, hereafter "peak MECO"), followed by return to pre-event values (Figure 4). Despite some discrepancies (Figures 4b-4d), our lower-resolution sites generally corroborate the observed trends at Site 865 well, especially when considering the stratigraphic challenges of correlating each record (see, e.g., Figure S1). With no further assumptions, this drop in δ 11 B values constitutes the first empirical evidence for surface ocean acidification over the MECO and at face value supports a rise of atmospheric pCO 2 associated with the event. Notably, though, the δ 11 B excursion at peak MECO where our samples are best resolved is brief (~120 kyr from onset to recovery), with limited, if any, change during the preceding 280 kyr of global δ 18 O decline (~40.5-40.22 Ma, hereafter "MECO onset").

Calculated pH Change at the MECO
To examine whether our choices of vital effect calibrations are reasonable, we compare pH calculated from δ 11 B with and without consideration of vital effects in Figure 5. We observe very large latitudinal offsets in δ 11 B values (~3.3‰; Figure 4). Without applying calibrations for vital effects at the MECO, even accounting for temperature (and hence pK * B ) variations with latitude, gives a latitudinal pH gradient of up to 1.2 (Figure 5a), that cannot reasonably be explained even considering a possibly deeper depth habitat for

10.1029/2019PA003713
Paleoceanography and Paleoclimatology G. index. Such a large pH offset would translate to considerable [CO 2 ] aq variation between sites and CO 2 disequilibrium between ocean and atmosphere of >4,500 μatm CO 2 at Site 702 (for alkalinity >1,000 μmol/kg), which is not plausible. Besides this, at the pH values inferred at Site 702, for an ocean alkalinity ranging from 1,000-3,000 μmol/kg, Ω calcite would reach as low as 0.2-0.4 at peak MECO. This would require G. index and other calcifiers at Site 702 to have persisted in waters highly undersaturated with respect to calcite, with few if any analogues for such behavior today. In contrast, our pH estimates agree closely once vital effects are accounted for as described above (Figure 5b), suggesting our assignment of calibrations from modern analogue species according to δ 13 C and δ 18 O is reasonable. These data therefore indicate that vital effects in foraminiferal δ 11 B were present at least in middle Eocene globigerinathekids and thus may well have been present in other extinct species too. We note, however, that the magnitude of vital effects in A. praetopilensis and M. crassatus predicted by the Foster et al. (2012) T. sacculifer calibration is very small, and so the agreement between calibrated globigerinathekids and uncalibrated A. praetopilensis and M. crassatus is not significantly worse than when modern analogue calibrations are applied to all species (see Figure S4). Further research, considering among other things changes in δ 11 B with size fraction (see, e.g., Henehan, Foster, et al., 2016), is therefore required to determine whether modern analogue calibrations must be applied to Morozovella and Acarinina species.
As with any geological observation, corroboration of a signal at multiple sites globally allows more confidence in results. Additionally, multisite analyses allow us to test the feasibility of vital effect corrections in the Eocene (Figure 5), reducing an otherwise large structural uncertainty in our pH and pCO 2 estimates. However, combining data from four geographically disparate sites as we do here does present complications, with the obvious introduction of aliasing (Pisias & Mix, 1988) from our low-resolution sites (e.g., ODP Site 1260, where the peak MECO is not sampled), noise from uncertainty in stratigraphic correlation (e.g., see Figure S1), and possible regional hydrographic differences (e.g., variable productivity, and hence perhaps air-sea disequilibrium, at ODP Site 1263 over the MECO;Boscolo Galazzo et al., 2015). In all likelihood, all of these processes are at work within our records. The intersite pH variability, taken at face value, at any given time within our record ( Figure 5) is sometimes larger than could easily be simulated with earth system models, suggesting that inaccuracy in stratigraphic correlation likely plays a role in adding noise to our record. Consequently, so as not to risk overinterpreting scatter from these various sources of error at any one site and to facilitate modeling work, we calculate a smoothed, globally averaged surface ocean ΔpH LOESS curve for the MECO (Figure 6).
To calculate the LOESS curve, we combined all sites' δ 11 B-pH measurements using a best-fit LOESS where the degree of smoothing was optimized using a generalized cross-validation algorithm (Golub et al., 1979) which objectively gauges the most statistically supported value for the LOESS span (using code from M. Friendly, http://www.datavis.ca).
Because the rates of pH change are important, we incorporated some uncertainty in the age domain, based on the sedimentation rate estimates between tie points, along with uncertainty in calculated pH. Specifically, 500 independent LOESS lines were fitted to each of the 500 replicate Monte Carlo pH data sets outlined above, with 500 simulated ages for each data point generated from within the age range of the sampling interval. For a δ 11 B sw of 38.5-38.9‰, our data indicate a pH decrease of 0.12 ± 0.04 over <100 kyr coincident with the peak MECO ( Figure 6). The objective smoothing algorithm used in our global LOESS fit does however dampen the transient reduction of up to~0.18 pH units implied by a single point at our most highly resolved site, ODP 865 (Figure 4), and so it is possible that further higher resolution studies may find our estimates of pH change to be conservative.
A somewhat surprising aspect of our LOESS global mean pH record is the lack of a pronounced decline in mean ocean pH during the onset phase of the MECO, when global temperature records suggest gradual warming (e.g., Bohaty et al., 2009). There are some apparent discrepancies between sample sites during this interval: At ODP Site 1263, δ 11 B drops more in line with global δ 18 O decline than observed at Sites 865 or 1260. However, since higher resolution sites show no such rise, and statistically none of the four sites show a significant trend in pH before the peak of the MECO ( Figure S5), we avoid overinterpreting intersite differences here.
3.3. pCO 2 Change Over the MECO pCO 2 calculated from surface ocean pH is shown for a constant saturation state ( Figure 7a) and a constant LOSCAR-derived TAlk of 1,750 ± 150 μmol/kg (Figure 7b). In the case of a constant Ω, pre-MECO pCO 2 (as defined in this instance as the mean of the LOESS values during the earliest 17 pre-MECO data points) is~669 ± 145 μatm, with the maximum of the pCO 2 LOESS rising to~1,062 ± 127 μatm at the peak MECO (equating to 0:67 þ0:17 − 0:12 doublings). Assuming a constant (±150 μmol/kg) alkalinity instead, our boron isotope data would suggest a pCO 2 rise from~563 ± 67 μatm to~770 ± 62 μatm at the peak MECO, a rise that is~47% smaller than if calcite saturation remained constant. As these two end-member scenarios represent extreme limits to the possible real-world response, we suggest that the pCO 2 increase during the MECO was likely somewhere between these two scenarios.
From our LOSCAR simulations (Figure 8), it may be seen that while absolute values of pCO 2 calculated are sensitive to δ 11 B sw (Figure 8b), relative change in greenhouse forcing in terms of CO 2 doublings is relatively insensitive to changes in this parameter (i.e., 0.55 to 0.75 CO 2 doublings for δ 11 B sw between 37.5‰ and 39.5‰; Figure 8c). This has been similarly demonstrated at the Eocene-Oligocene transition (Pearson et al., 2009) and at the PETM (Penman & Zachos, 2018). Using our LOSCAR modeling approach, instead of assuming constant Ω calcite or alkalinity, our best-estimate middle Eocene δ 11 B sw range of 38.5-38.9‰ produces ΔpCO 2 scenarios from pre-MECO (i.e., background) to the peak MECO of 500 μatm rising to 765 μatm for a δ 11 B sw of 38.5‰ and 630 μatm rising to 990 μatm for a δ 11 B sw of 38.9‰ (Figure 8b). Expressed as doublings of CO 2 , these end-member scenarios for ΔpCO 2 equate to between 0.62 and 0.67, which translates to a

10.1029/2019PA003713
Paleoceanography and Paleoclimatology CO 2 forcing change (ΔF) of between 2.3 and 2.5 W m −2 (Figure 8c). To quantify the absolute uncertainty bounds on pCO 2 within our best-estimate range of δ 11 B sw , we repeated these LOSCAR modeling exercises for 1 and 2 SE of the 500 Monte Carlo LOESS values for pre-event pH and peak MECO pH minimum (i.e., gray shaded regions in Figure 6), to calculate a pre-event pCO 2 of 565 þ110 − 80 , and a peak MECO pCO 2 of 880 þ130 − 100 , as plotted in Figure S6. These pre-MECO pCO 2 levels are somewhat lower in an absolute sense compared to those values calculated assuming constant, pre-ascribed Ω calcite . This is because our LOSCAR simulations produce values of Ω calcite in the surface Pacific of~5 throughout the MECO, which is on the lowest end of our assumed feasible Eocene range of Ω calcite and thus implies lower equilibrium TAlk values at this time. However, despite this difference in absolute values, a rise of 0:67 þ0:17 − 0:12 doublings for the constant Ω calcite scenario falls squarely in the range of estimates produced by LOSCAR. By contrast, assuming constant alkalinity produces a smaller rise in pCO 2 upon a similar pre-MECO baseline (0:45 þ0:06 − 0:05 ), because without any concurrent rise in TAlk, a smaller rise in DIC is required to affect the same change in pH.

MECO CO 2 Change and Carbon Cycling
Our new boron isotope record helps to resolve some enigmatic features of the MECO. For example, an abrupt (<10 kyr) rise in atmospheric pCO 2 centered around the peak of the MECO (Figure 7) could help to explain why the CCD shoaled even at the peak of a~400 kyr interval of apparent gradual warming-a key element of the so-called MECO "carbon cycle conundrum" (Sluijs et al., 2013). A shorter CO 2 pulse could have surpassed the timescales on which chemical weathering responds to maintain oceanic CaCO 3 saturation state (~10 4 -10 5 year) (Colbourn et al., 2015), thereby prompting deep-sea carbonate dissolution similar to that seen during hyperthermals (e.g., Zachos et al., 2005). The source of such a carbon injection remains enigmatic however, particularly given the lack of any global excursion toward negative δ 13 C at this time (Bohaty et al., 2009;Boscolo Galazzo et al., 2014;Edgar et al., 2010). Our lower estimates for background pCO 2 and change in pCO 2 over the MECO relative to previous reconstructions (Bijl et al., 2010) alleviate this problem somewhat, in necessitating lower masses of CO 2 injection into, and removal from, the surface ocean-atmosphere system. This also helps to reconcile a rapid cooling and drawdown of CO 2 after the MECO with the absence of any carbonate overshoot due to silicate weathering, such as is seen at the PETM . We note that enhanced burial of organic carbon, such as is observed in the Peri-Tethys of Italy (Spofforth et al., 2010) and the Crimea-Caucasus (Benyamovskiy, 2012), and consistent with globally increased δ 13 C in carbonates following the MECO (Bohaty et al., 2009), may have also played an important role in carbon drawdown.
Our finding that there was limited change in atmospheric pCO 2 during the onset phase of the MECO, however, even while carbonate δ 18 O (Bohaty et al., 2009) and biomarker-based proxies (Bijl et al., 2010;Cramwinckel et al., 2018) suggest global temperatures were rising steadily over >200 kyr is a notable feature. To place a probabilistic bound on the maximum pCO 2 rise permitted by our data during the onset interval, we used a Monte Carlo approach, simulating 1,000 data sets of pH across the "MECO onset" (40.5-40.21 Ma) at our preferred estimate of δ 11 B sw by randomly sampling within the range of uncertainty of each data point. We couple these to a constant approximate pre-MECO surface ocean alkalinity estimate from LOSCAR In panel b, the atmospheric pCO 2 required for LOSCAR to attain "pre-event" pH values (solid red line and triangles) is shown. We then iteratively simulated injection of different carbon masses over the timescales of our pH drop (~100 kyr) until we attained a surface Pacific pH that matched the pH minimum in our LOESS curve calculated for that δ 11 B sw and noted the atmospheric pCO 2 attained by LOSCAR (dashed red line and inverted triangles). Panel c shows the change in pCO 2 in panel b, expressed in terms of CO 2 doublings, and as CO 2 forcing change, ΔF (calculated according to Myhre et al., 1998). Our best estimate of δ 11 B sw (38.5-38.9‰) is shown by the vertical blue bar through (a), (b) and (c) and implies 0.6-0.7 doublings of pCO 2 over the MECO event.
(1,750 ± 150 μmol/kg) to calculate the trend in pCO 2 predicted over the MECO onset interval. Using this approach, the mean simulated ΔpCO 2 over the MECO onset was +29 μatm. In the parlance of the IPCC (IPCC Core Writing Team, 2014), we find that an increase of >230 μatm (~0.4 doublings) would be very unlikely (i.e., <5% chance, see Figure S7). Reconciling this small a magnitude of pCO 2 change with warming of 3°C, we suggest, would require one (or a combination of) (a) very high climate sensitivity to CO 2 in the middle Eocene and/or some non-CO 2 climate forcing, (b) dynamic changes in ocean alkalinity that we do not properly account for during the onset of the MECO, or (c) a nonthermal (i.e., ice volume) component to carbonate δ 18 O change.
While calculating a precise estimate of climate sensitivity is difficult given the limited spatial coverage of available surface temperature records, the observed~0.5‰ decrease in benthic foraminiferal δ 18 O values could equate to~2°C global mean surface temperature change using the relationships for the pre-Pleistocene described by Hansen, Sato, Russell, and Kharecha (2013), assuming no ice volume component. Using our upper limit of 0.3 CO 2 doublings during the MECO onset, we would estimate an Earth system sensitivity of >6.7°K per CO 2 doubling. Such sensitivity would be far higher than model estimates of equilibrium climate sensitivity (IPCC Core Writing Team, 2014) and well outside of the range of climate sensitivities previously calculated over the Cenozoic (Rohling & Members, 2012). Moreover, any such heightened climate sensitivity would have to have lowered again at the peak MECO; otherwise, a 0.65 doubling of pCO 2 would have produced >4°C of warming on top of the gradual MECO warming trend, something which is not observed in marine temperature records.
Large-scale changes in ocean alkalinity are another hypothetical candidate for raising atmospheric CO 2 without detectable change in surface pH (and thus foraminiferal δ 11 B) during the MECO onset. A reduction in surface ocean TAlk might be possible over the MECO interval if there were a reduction in silicate weathering, as suggested by Os isotopes and the observed shallowing of the CCD (van der Ploeg et al., 2018). However, without a concurrent reduction in DIC, such a scenario would be associated also with a reduction in global surface ocean pH, for which we see little evidence. For surface ocean aqueous pCO 2 , and hence atmospheric pCO 2 , to have risen without significant pH change, TAlk and dissolved inorganic carbon (DIC) would have had to rise in a ratio of approximately 1:0.9. A process (or processes) by which such stoichiometry would arise is unclear and might necessitate a complex combination of weathering fluxes and organic carbon burial fluxes. However, all else being equal, such a rise in ocean TAlk and DIC would have increased deep ocean carbonate saturation state, counter to observed MECO CCD change indicative of decreased saturation. Moreover, alkalinity would have to increase by~1,000 μmol/kg and [DIC] bỹ 900 μmol/kg to increase pCO 2 by 300 μatm: magnitudes which seem unfeasibly large to reconcile with a lack of pronounced shift in δ 13 C. Thus, while possible factors controlling alkalinity such as dynamic changes in calcification fluxes (Boudreau et al., 2018;Henehan, Hull, et al., 2016) or sulfide oxidation (Calmels et al., 2014;Torres et al., 2017) merit further investigation with more complex earth system models, at this stage, it seems difficult to decouple pH and pCO 2 during the MECO onset given the constraints of the marine carbonate cycle.
Alternatively, it is possible that there was some contribution from ice melt to the trend toward more negative δ 18 O during the MECO onset, as has been hypothesized elsewhere (e.g., Lyle et al., 2005;Tripati et al., 2005). Besides amplifying changes in carbonate δ 18 O, ice growth and subsequent melt can also provide a mechanism to explain CCD change, by changing the balance of carbonate burial onto shelves versus the deep sea through sea level fluctuation (Lyle et al., 2005;Sluijs et al., 2013). As yet, however, evidence for earlier, pre-MECO glaciation is only suggestive and not conclusive. Certainly, some high-latitude influence on climate and carbon cycle at this time is evidenced by a strengthened influence of obliquity in climate records (Bosboom et al., 2014;Westerhold et al., 2014). Furthermore, the second "MECO-like" event shown in Figure 1 does roughly coincide with a known ephemeral Antarctic glaciation event at~37.3 Ma (Scher et al., 2014, and references within), perhaps by analogy lending some suggestive support for dynamic ice also around the MECO. A steeper slope in the relationship between high-latitude TEX 86 and δ 18 O in the middle Eocene compared to the early Eocene (Bijl et al., 2009) and comparison of Mg/Ca and δ 18 O (Billups & Schrag, 2003;Lear, 2000) could also support an additive nonthermal component to seawater δ 18 O fluctuations beginning prior to the MECO. Furthermore, although glaciation thresholds are highly model-specific (Gasson et al., 2014), most general circulation models (GCMs) simulate some continental ice on Antarctica (e.g., DeConto & Pollard, 2003;Gasson et al., 2014) at pre-MECO pCO 2 (e.g., Anagnostou et al., 2016). Temperature proxy records from the distal Antarctic peninsula (paleolatitudẽ 70°;van Hinsbergen et al., 2015) indicate coastal mean annual temperatures of~10-15°C prior to the MECO (Douglas et al., 2014;Judd et al., 2019), which when scaled to the continental interior should have been sufficiently cold to prompt glacial inception in the Gamburtsev mountains (Rose et al., 2013). Additionally, various lines of sedimentological, micropaleontological, and seismographic evidence (of varying degrees of certitude) exist for alpine or marine-terminating glaciers on Antarctica prior to or around the MECO (e.g., Birkenmajer, 1991;Ehrmann & Mackensen, 1992;Eittreim et al., 1995;Gulick et al., 2017;Margolis & Kennett, 1970), although we note that recently some of these interpretations of seismic features have been contested (Sauermilch et al., 2019). Further, while we stress that extensive northern hemisphere ice sheets were not present (Edgar et al., 2007), indicators for middle Eocene Arctic sea ice (Darby, 2014;Eldrett et al., 2007;Stickley et al., 2009;Tripati et al., 2008) would seem incongruous with an ice-free Antarctic interior. Finally, the onset of synchronous changes in inferred sea level and δ 18 O (Browning et al., 1996;Pekar et al., 2005) that follow 1.2 Myr obliquity cycles (Boulila et al., 2011) coincides with an inflection in the relationship between estimated sea level and global deep-sea temperatures around 42-44 Ma (Gasson et al., 2012;Kominz et al., 2008;Lear, 2000) that would also support ice volume change. We note though that this is strongly reliant on Kominz et al.'s New Jersey margin sea level reconstruction being representative of global eustasy (Gasson et al., 2012), and in such marginal marine settings, it is difficult to confidently distinguish global eustatic change from steric effects or sediment supply, especially given the flatter shelves likely present at this time (Sømme et al., 2009). In If middle Eocene ice reservoirs were present and were sensitive to orbital configuration, changes in ocean circulation, or small rises in pCO 2 , ice melt during the MECO onset would have lowered δ 18 O sw and raised sea level, and partitioned more carbonate deposition onto flooded shelves and away from the deep ocean (Lyle et al., 2005;Sluijs et al., 2013;Tripati et al., 2005). Such a scenario could therefore reconcile our observed stable surface ocean pH, shoaling CCD and a decline in carbonate δ 18 O during the MECO onset. On one hand our new estimates of pre-MECO pCO 2 fall well below most modeled thresholds required for significant Antarctic glaciation (although model dependent; Gasson et al., 2014), and evidence for sizeable glacio-eustatic fluctuations around this time (e.g., Browning et al., 1996) suggests this is not wholly unreasonable. On the other hand, as previously discussed, conclusive evidence for substantial middle Eocene ice reservoirs is still lacking. Published organic biomarker temperature proxies (TEX 86 and U k 37 ) should allow another means to test this hypothesis. In support of an elevated Earth system sensitivity scenario, bulk carbonate δ 18 O and organic temperature proxies at ODP Site 1172 move more or less in unison (Bijl et al., 2010) during the period in which we see little change in global ocean pH. Lower-resolution TEX 86 data over the MECO at nearby ODP Site 1170 also show MECO warmth, albeit more muted . However, a possible incursion of the warm East-Australian current into this area around the MECO is indicated by fossil assemblages (Cramwinckel, Woelders, et al., 2019), which may complicate these signals. At equatorial ODP Site 959 a clear MECO warming event is expressed in TEX 86 , but differing sample resolution makes it difficult to discern how clear the correlation between bulk carbonate δ 18 O and TEX 86 is at this site . Additionally, a lack of age model tie points between 40.02 and 42.84 Ma at this site makes it difficult to definitively discern whether the warming expressed in TEX 86 was gradual and began during the global MECO onset or was rapid and focused around the peak MECO, when our boron isotope data indicate CO 2 rise. Finally, at ODP Site 1263, there is no apparent correlation between TEX 86 and concurrent carbonate δ 18 O values (Boscolo Galazzo et al., 2014). As yet, then, evidence is not sufficiently conclusive to definitively accept or reject the hypothesis of ice melt contribution to δ 18 O change during the onset of the MECO. However, the extremely high climate sensitivity to pCO 2 otherwise implied by our boron data suggests that a component of ice melt is a hypothesis worthy of further interrogation.

Caveats, Conclusions, and Prospects
By combining with LOSCAR simulations, we derive boron-based estimates of pCO 2 between 40.21-41.0 Ma of~550 μatm (2σ þ110 − 80 μatm; see Figure S6) that provide a precise estimate of the boundary conditions on which the MECO was superimposed. Additionally, our estimate of peak MECO pCO 2 of~870 μatm ( 2 σ þ130 − 100 μatm; see Figure S6) tightly constrains CO 2 forcing change (ΔF) to between 2.3 and 2.5 W m −2 (Figure 8c). These data confirm the finding of Bijl et al. (2010) that there was a rise in pCO 2 over the MECO, albeit to a lower degree and superimposed upon a lower baseline pCO 2 than indicated by the alkenone δ 13 C-pCO 2 proxy. A number of outstanding puzzles surrounding the MECO clearly remain, however. The question of reduced silicate weathering at the MECO (van der Ploeg et al., 2018) remains open, for example, and requires further investigation. A reduced silicate weathering feedback would have implications for our LOSCAR-derived pCO 2 estimates, in that we use default configuration weathering in LOSCAR to predict the evolution of TAlk expected over the course of a 100 kyr CO 2 injection, allowing us to convert pH to pCO 2 . As demonstrated in Figure 7, the effect of a reduced silicate weathering feedback over the MECO would be to necessitate a smaller mass of CO 2 to be input to effect the same change in ocean pH and δ 11 B. However, there are a number of reasons why we do not attempt to explore this in greater depth with LOSCAR here. First, while a temporarily weakened and subsequently reactivated silicate weathering feedback is an elegant hypothesis, it is not the only solution that could explain changes in Os isotopes and changes in CCD. For example, rising sea level could conceivably partition more carbonate deposition onto shelves and shoal the CCD (as already demonstrated in LOSCAR; Sluijs et al., 2013) while simultaneously changing the relative contributions of radiogenic versus nonradiogenic rocks to global weathering fluxes. Effects on marine Si cycling, an often-overlooked implication of mechanisms that invoke a dynamic silicate weathering feedback, could provide an independent means to interrogate the hypothesis of reduced silicate weathering at the MECO. During earlier Paleocene and Eocene warming events, the silicate weathering feedback on land drove enhanced burial and preservation of biogenic silica in the oceans as Si input fluxes were elevated (Penman, 2016;Penman et al., 2019). By analogy, a weakened silicate weathering feedback prior to the MECO might predict reduced burial of silica in the middle Eocene ocean. However, opal mass accumulation rates in the Pacific were among their highest Eocene levels in radiolarian zones RP14 and RP15, preceding the MECO (Moore et al., 2008). Moreover, the response to the MECO in the North Atlantic (Witkowski et al., 2014), the Southern Ocean (Witkowski et al., 2012), and the equatorial Indian Ocean (Savian et al., 2016) was a pronounced rise in the deposition/preservation of siliceous plankton, supportive perhaps of a global rise in marine [Si] concurrent with warming. These records do not therefore immediately appear consistent with a weakened silicate weathering feedback at the MECO. However, we would encourage further investigation into the magnitude and geographical extent of this enhanced silica burial and to look for chert deposits analogous to those seen after earlier warming events (Penman, 2016;Penman et al., 2019) to independently test the weakened silicate weathering hypothesis. Furthermore, interrogating other weathering-sensitive isotope systems such as Li, Si, or Sr may help in further constraining changes in silicate weathering over this event.
Besides questions over the dynamics of the silicate weathering feedback, we suggest that even implementing a weakened and then restored silicate weathering feedback in LOSCAR (as in van der Ploeg et al., 2018) would not fully capture the complexity of the MECO event. In such simulations (e.g., Sluijs et al., 2013;van der Ploeg et al., 2018) CO 2 drawdown is not as rapid as the recovery of pH and global temperatures observed after the MECO and is usually accompanied by an overshoot in simulated CCD, for which (at least in the Pacific and Indian Oceans) there is scant sedimentary evidence. This may suggest a role for organic matter burial (Moebius et al., 2015;Spofforth et al., 2010) and/or ocean carbon storage (Sexton et al., 2011). To simulate these forcings realistically, and to accurately reproduce MECO surface ocean δ 13 C changes (which vary widely both in magnitude and direction of change globally), would require a more complex model with dynamic ocean circulation. In short it is likely that during the MECO and its recovery more complex biogeochemical feedbacks were at work than are represented in the LOSCAR simulations of Sluijs et al. (2013), Van der Ploeg (2018), or our current contribution. Therefore, we present our scenarios for MECO pCO 2 change as current best estimates, mindful that lower TAlk change due to weakened silicate weathering might reduce the required mass of C injection but that higher sampling resolution around the short pH minimum (statistically underweighted in our global LOESS as it is defined by a single data point) might increase it.
In summary, these data offer robust independent confirmation of pCO 2 rise during the MECO, the bulk of which occurs at the peak of the 400 kyr event. Although our new pCO 2 estimates do not allow us to ascertain the causal driver of the MECO, a number of informative observations can be made. First, the mass of C that was transferred to and from the surface ocean-atmosphere system over the MECO was considerably smaller than that required to explain previous MECO ΔpCO 2 estimates based on the alkenone δ 13 C proxy (Bijl et al., 2010), and if the silicate weathering feedback were diminished, this mass could be smaller still.
This smaller mass could bring internal reorganization of carbon reservoirs into play in explaining the MECO, rather than invoking a large exogenous carbon injection. Second, the limited evidence we see globally for any substantial rise in pCO 2 during the onset stage of the MECO suggests a climate, and perhaps nascent cryosphere, that was highly sensitive to perturbation. While we suggest that a contribution to global δ 18 O change from ice melt during the initial phase of the MECO may be the most parsimonious means to reconcile CCD shoaling and limited pCO 2 rise, we recognize that further work is required to ascertain the feasibility of such a scenario. For instance, more precise dating of glacial deposits and more extensive circum-Antarctic temperature proxy records are of vital importance. Furthermore, we encourage further comparison of δ 18 O with concurrent independent temperature proxies (e.g., TEX 86 and clumped isotopes) to further isolate δ 18 O sw change over the MECO. Modeling work to explore the potential δ 18 O range of nascent middle Eocene ice would also be beneficial, to constrain the volumes of ice melt required. Finally, while a trigger for ice melt without significant pCO 2 rise is uncertain, we suggest that ocean circulation is a prime target for investigation, given (a) spatially variable bulk carbonate δ 13 C changes (Bohaty et al., 2009), (b) spatial variability in circum-Antarctic temperature (Douglas et al., 2014;Judd et al., 2019), (c) geochemical and micropaleontological evidence for changing high-latitude ocean circulation (Cramwinckel, Woelders, et al., 2019;Scher & Delaney, 2010), and (d) enhanced sensitivity of global thermohaline circulation to orbital configuration (Vahlenkamp et al., 2018) around the MECO. In sum, while our new estimates of pCO 2 across the MECO move our understanding of this event forward, there remains considerable work to be done to fully elucidate the drivers of this enigmatic event.