Constraining the evolution of Neogene ocean carbonate chemistry using the boron isotope pH proxy

Abstract Over the course of the Neogene, the Earth underwent profound climatic shifts from the sustained warmth of the middle Miocene to the development of Plio-Pleistocene glacial–interglacial cycles. Major perturbations in the global carbon cycle have occurred alongside these shifts, however the lack of long-term carbonate system reconstructions currently limits our understanding of the link between changes in CO2, carbon cycling, and climate over this time interval. Here we reconstruct continuous surface ocean pH, CO2, and surface ocean aragonite saturation state using boron isotopes from the planktonic foraminifer Trilobatus trilobus and we perform a sensitivity analysis of the key variables in our calculations (e.g. δ 11 Bsw, [Ca]sw, CCD). We show that the choice of δ 11 B sw influences both seawater pH and CO2 while [Ca]sw reconstructed dissolved inorganic carbon exerts a significant influence only on CO2. Over the last 22 Myr, the lowest pH levels occurred in the Middle Miocene Climate Optimum (MMCO; 17–14 Myr ago) reaching ∼ 7.6 ± 0.1 units in all our scenarios. The extended warmth of the MMCO corresponds to mean CO2 and aragonite saturation state levels of 470–630 ppm and 2.7–3.5, respectively. Despite a general correspondence between our CO2 record and climate, all CO2 scenarios show a peak at ∼9 Ma not matched by corresponding changes in climate reconstructions. This may suggest decoupling (i.e. significant CO2 change without a discernible climate response) for a limited interval in the Late Miocene (11.6–8.5 Ma), although further refinement of our understanding of the temporal evolution of the boron isotopic composition of seawater is necessary to fully evaluate the nature of the relationship between CO2 and climate. Nonetheless, from our long-term view it is clear that low-latitude open ocean marine ecosystems are unlikely to have experienced sustained surface pH and saturation levels below 7.7 and 1.7, respectively, during the past 14 million years (66% CI).


Introduction
Boron isotopes (δ 11 B) in marine carbonates have been exploited as a proxy for the marine carbonate system (i.e. pH, pCO 2 ) on multiple timescales within the Cenozoic (e.g. Sanyal et al., 1995;Palmer et al., 1998;Pearson and Palmer, 2000;Foster et al., 2012;Penman et al., 2013;Anagnostou et al., 2016;Gutjahr et al., 2017;Chalk et al., 2017). Plio-Pleistocene foraminiferal δ 11 B derived reconstructions exhibit glacial-interglacial variability in pH and pCO 2 , in good agreement with independent records (e.g. ice cores) (Seki et al., 2010;Bartoli et al., 2011;Hönisch et al., 2009;Henehan et al., 2013;Martínez-Botí et al., 2015a; have been determined using Earth System modelling (Anagnostou et al., 2016;Gutjahr et al., 2017). According to Anagnostou et al. (2016), the early Eocene CO 2 concentration was 1400 ± 470 ppm decreasing to 550 ± 190 ppm in the early Oligocene. The uncertainty window on these estimates reflects uncertainty in calcite saturation and δ 11 B sw estimates, δ 11 B c measurement error, and possible atmosphere-ocean disequilibrium. Despite these relatively large uncertainties, it is clear that early Cenozoic pCO 2 was much higher than in the Plio-Pleistocene and CO 2 change was an important driver of climate change in the early Cenozoic. However, application of the δ 11 B-pCO 2 proxy during the more recent Neogene (2.6 to 23 Myrs ago) transition to a bipolar icehouse world is complicated by the fact that most reconstructions indicate a relatively dramatic shift in δ 11 B sw in the late Miocene to Pliocene. While there is first-order agreement between different records for the late Miocene δ 11 B sw rise, the exact timing, magnitude and rate of change of this shift remains uncertain (Pearson and Palmer, 2000;Lemarchand et al., 2002;Raitzsch and Hönisch, 2013;Greenop et al., 2017), with implications for the structure of δ 11 B-based reconstructions of pH and pCO 2 during the period. Furthermore, the history of carbonate chemistry (i.e. DIC and alkalinity) in this interval is poorly understood and currently is limited to model estimates that use pCO 2 estimates from GEOCARB carbon cycle models and/or paleo-pCO 2 proxy estimates paired with estimates of the calcite compensation depth (CCD) (Tyrrell and Zeebe, 2004;Ridgwell, 2005).
During the Neogene period, the Earth underwent profound changes from the warmth of the Middle Miocene Climatic Optimum (∼14-17 Myrs ago; MMCO) to the expansion of the Antarctic ice sheet at 13.9 Ma (Lear et al., 2010). Dynamic changes in the carbon cycle have been proposed to be associated with these climatic shifts, involving changes in silicate weathering, volcanic emissions and the organic carbon subcycle (Derry and France-Lanord, 1996;France-Lanord and Derry, 1997;Raymo and Ruddiman, 1992;Raymo, 1994;Foster et al., 2012). The evolution of the carbonate system during this time interval therefore represents a critical gap in our understanding of climate-carbon cycle interactions. Previous CO 2 reconstructions derived from a disparate set of proxies (e.g. plant stomata, marine alkenones, paleosols, marine carbonate δ 11 B) have produced a wide range of estimates, yet with a general consensus that MMCO CO 2 levels were higher than preindustrial, and that CO 2 subsequently declined across the mid-Miocene climate transition (Kürschner et al., 2008;Foster et al., 2012;Zhang et al., 2013). Since CO 2 and pH are highly correlated variables in the ocean carbonate system, the discrepancies between existing CO 2 reconstructions (Pagani et al., 1999a(Pagani et al., , 1999bFoster et al., 2012;Zhang et al., 2013) imply a lack in understanding of the ocean's history of pH change and the evolution of the carbon cycle and also, to some extent, an uncertainty in proxy systematics. Further, the degree of coupling between climate and CO 2 in the late Miocene is poorly constrained, limiting our ability to determine the importance of CO 2 variability in setting the evolution of Miocene climate (Pagani et al., 1999a(Pagani et al., , 1999bLaRiviere et al., 2012). Thus, the lack of available records that constrain the ocean's carbonate chemistry is a critical impediment to establishing the relationship between the longterm Earth System evolution of carbon cycle and climate (Hönisch et al., 2012). Additionally, without detailed knowledge of longterm changes in seawater carbonate chemistry we are blind to its role in driving evolutionary adaptation and resilience of marine organisms to environmental change (Hönisch et al., 2012;Hain et al., 2015).
Here we use a new dataset of 40 planktic boron isotope measurements in conjunction with existing benthic data and new constraints on the evolution of seawater δ 11 B (Lemarchand et al., 2002;Raitzsch and Hönisch, 2013;Greenop et al., 2017) to recon-struct surface pH and ocean carbon content for the past 22 million years. Combining our surface pH record with independent reconstructions of seawater calcium and magnesium concentrations, changes in the calcite compensation depth, deep-sea pH, and temperature, we calculate both atmospheric CO 2 and surface aragonite saturation state (Ω arag ). We also discuss the relative importance of the variable sources of uncertainty in these parameters on our reconstructions. With extensive uncertainty propagation using a Monte Carlo approach and sensitivity analysis we reconstruct longterm pH, aragonite saturation state and atmospheric CO 2 since the MMCO (0-17 Myrs ago). Our new record of ocean carbonate chemistry allows us to examine the role of CO 2 in the climate system over the course of the Neogene and contextualize future carbon emission pathways against this geological perspective.

Materials and methods
The boron isotope proxy relies on the principle that the ratio of the dominant boron species in seawater, borate ion and boric acid, varies as a function of pH (Dickson, 1990), which also controls the isotopic composition of the borate ion that is incorporated into a marine carbonate such as a planktic foraminiferal calcite (Vengosh et al., 1991;Hemming and Hanson, 1992;Foster and Rae, 2016). The foraminifer-derived pH record may then be used to estimate CO 2 levels and Ω arag using assumptions about other carbonate system parameters. Early attempts at reconstructing pH and atmospheric CO 2 from fossil foraminifera over this interval (Spivack et al., 1993;Palmer et al., 1998;Pearson and Palmer, 2000) are qualitatively useful but have been superseded by a revised understanding of boron isotope systematics and analytical developments (see discussion in Seki et al., 2010;Foster et al., 2012;Anagnostou et al., 2016).

Sample locations and site details
In this study we analysed the CaCO 3 shells of a single planktic foraminifera species, Trilobatus trilobus. This species was formerly assigned to the genus Globigerinoides but is no longer considered part of that genus (Spezzaferri et al., 2015). The species is sometimes referred to as T. sacculifer, especially those individuals that bear a final sac-shaped final chamber, but the feature is ecophenotypic (André et al., 2013) and the name trilobus has priority. In this study we avoided individuals with sac-shaped chambers. Unusually for planktic foraminifera, T. trilobus displays no significant cryptic genetic diversity (André et al., 2013). It lives in the upper part of the water column in association with dinoflagellate photosymbionts, although it probably sinks towards the end of its life cycle at gametogenesis, somewhat complicating the shell chemistry (e.g. Birch et al., 2013). The species is widely used in boron isotope studies because it is one of the few to survive almost the entire Neogene and Quaternary periods (Spezzaferri et al., 2015) and a species-specific calibration for δ 11 B-pH exists (discussed further below; Sanyal et al., 2001;Seki et al., 2010;Foster et al., 2012;Greenop et al., 2014).

Age model
Our age model for ODP Site 761 uses biostratigraphic datums (Holbourn et al., 2004) on the CK95 timescale (Cande and Kent, 1995;Supplementary Table S1). The age model for ODP Site 926 is a polynomial fit through nannofossil and planktic foraminifer biostratigraphic datums based on the shipboard report (Site 926 Shipboard Scientific Party, 1995) and post-cruise refinement Pearson and Chaisson, 1997) adapted to the CK95 timescale. Ages for ODP Site 872 was calculated by linear interpolation between reliable biostratigraphic datums (Pearson, 1995) from the early to mid-Miocene. Note that for ODP Site 872, age estimates for sample depths 105.71 and 106.79 mbsf have large uncertainty (±1 Ma) as evidence for reworking is present in this section, however this is not significant at the resolution of this study for the early Miocene. Nannofossil and planktic foraminifer biostratigraphic datums tied to the CK95 timescale were used for the age model at ODP Site 1000A (Site 1000Shipboard Scientific Party, 1997a, 1997b. The age models from published boron isotope records from ODP Sites 761 Greenop et al., 2014), 926 , and the Blue Clay Formation (Badger et al., 2013) were put on the CK95 timescale in a similar manner using shipboard reports, published biostratigraphic datums (Holbourn et al., 2004;Site 926 Shipboard Scientific Party, 1995;Site 999 Shipboard Scientific Party, 1997a, 1997bHilgen et al., 2009). Previously published records from the Plio-Pleistocene (<5 Myr) (Hönisch et al., 2009;Seki et al., 2010;Bartoli et al., 2011;Martínez-Botí et al., 2015a;Stap et al., 2016;Chalk et al., 2017) are on their original orbitally tied age models.

Analytical methodology
Between 100 and 150 shells of T. trilobus were picked from the 300-355 μm size fraction and cleaned following the oxidative cleaning method (Barker et al., 2003). Samples were dissolved in 0.5M HNO 3 and stored for boron isotope analysis in Teflon vials. Following Foster (2008) and as detailed in Foster et al. (2013), boron separation from the calcium matrix was performed using the boron specific resin Amberlite IRA 743. Boron isotope analysis was performed on a Thermo-Fisher Scientific Neptune multi collector inductively coupled plasma mass spectrometer (MC-ICPMS) at the University of Southampton with a sample-standard bracketing routine Henehan et al., 2013).
Calibration studies have shown that the borate species is predominantly incorporated in the calcium carbonate shell and thus the boron isotopic composition of borate can be used to estimate ocean pH (Vengosh et al., 1991;Hemming and Hanson, 1992;Sanyal et al., 2001). To calculate pH from the δ 11 B of T. trilobus we account for the species-specific offset between δ 11 B borate in seawater and δ 11 B of the foraminiferal calcite due to vital effects (equation (3)) (Foster and Rae, 2016;Henehan et al., 2016). We use the calibration slope from Sanyal et al. (2001) This species-specific calibration curve is used to estimate the δ 11 B of borate ion. The measurement uncertainty of δ 11 B trilobus estimated by the external reproducibility is fully propagated with an appropriately scaled normal distribution.
To estimate δ 11 B borate and surface ocean pH, we applied sizefraction and species-specific δ 11 B borate -δ 11 B calibrations (Supplementary Table S3). As previously published δ 11 B analyses in this compilation were performed using a range of techniques, we apply inter-instruments offsets to the published negative thermal ionization mass spectrometry (N-TIMS) data to correct to the MC-ICPMS and in doing so account for the well-documented instrument dependent offset in these NTIMS data Farmer et al., 2016) due to differences in chemical composition of sample (i.e. B/Ca) and other reasons that are currently unclear (Farmer et al., 2016). Some δ 11 B data from Hönisch et al. (2009) were performed on two different TIMS, the NBS design instrument at Stony Brook and Thermo Scientific Triton at LDEO. To account for this and correct the data to the MC-ICPMS, we use an offset of −1.1h, as determined from an interlaboratory comparison of five modern and fossil samples measured on both the NBS design and Triton NTIMS (Hönisch et al., 2009). To align the δ 11 B datasets collected on the Triton (Hönisch et al., 2009;Bartoli et al., 2011;Stap et al., 2016) with the Neptune MC-ICPMS used here, we adopt the −1.12h, offset presented in Farmer et al. (2016) based on average of 5 fossil samples measured on both the Triton NTIMS and Neptune ( Supplementary Fig. S1). Using this approach, the degree of the δ 11 B borate overlap between datasets from both instruments is still relatively poor in the region 2.5-3.0 Ma. To overcome this, we estimate pH and then tune the NTIMS offsets for overlapping pH datasets until there is a good agreement between NTIMS and MC-ICPMS. The best fit is achieved when we apply an additional −1.0h offset to the Triton planktic δ 11 B datasets. This additional offset suggests that further work is needed to provide a mechanistic understanding of the offset between MC-ICPMS and NTIMS techniques. We do not include the dataset of Stap et al. (2016) in the compilation as even with this treatment of the NTIMS offset the δ 11 B borate values exhibit substantial scatter and do not align with any of the other records ( Supplementary Fig. S1). We applied no offsets to the benthic δ 11 B data from Greenop et al. (2017) and Raitzsch and Hönisch (2013), measured on the MC-ICPMS and NTIMS, respectively, as the data agree well in their raw state (see Greenop et al., 2017).
Tandem trace element analyses were performed on a small aliquot (<10%) of each sample measured here on the Thermo-Fisher Scientific Element 2 single collector inductively coupled plasma mass spectrometer at University of Southampton following the method outlined in Henehan et al. (2013). Over the duration of this study, analytical reproducibility for Mg/Ca was ±2.7% (2σ ).
T. trilobus Mg/Ca values were converted to sea surface temperatures (SST) using the calibration of Anand et al. (2003) with a correction for seawater Mg/Ca (Horita et al., 2002;Brennan et al., 2013) following Hasiuk and Lohmann (2010) using a power constant of 0.41 (Delaney et al., 1985;Evans and Mueller, 2012) (Zachos et al., 2008); Labels indicate the timing of the Middle Miocene Climate Optimum (MMCO), Middle Miocene Climate Transition (MMCT), and Northern Hemisphere Glaciation (NHG); (B) Boron isotopic composition (δ 11 B) of planktonic foraminifer T. trilobus (300-355 μm) record from this study (blue) and published data from T. trilobus and G. ruber (Hönisch et al., 2009;Seki et al., 2010;Bartoli et al., 2011;Foster et al., 2012;Badger et al., 2013;Greenop et al., 2014;Martínez-Botí et al., 2015a;Chalk et al., 2017) (grey). Note the inverted axis with low δ 11 B corresponding to low pH and pCO 2 . Error bars denote analytical external reproducibility (2 s.d.); (C) Sea surface temperatures (SSTs) derived from foraminiferal Mg/Ca from this study (blue) and published SST data derived from Mg/Ca (Hönisch et al., 2009;Bartoli et al., 2011;Foster et al., 2012;Badger et al., 2013;Greenop et al., 2014;Martínez-Botí et al., 2015a;Stap et al., 2016;Chalk et al., 2017) and alkenone U k 37 values (Seki et al., 2010) (grey). Sites are coloured as in (B). (D) δ 11 B sw reconstructions from Lemarchand et al. (2002) with ±1h confidence interval, Raitzsch and Hönisch (2013) using the fractionation factor of Klochko et al. (2006) and 2 sigma uncertainty, and Greenop et al. (2017) with 95% confidence intervals. These reconstructions are referred to hereafter as "δ 11 B sw -L02", "δ 11 B sw -RH13", and "δ 11 B sw -G17". (E) Surface ocean pH reconstructions using the three δ 11 B sw reconstructions in (D). Note there is an additional δ 11 B sw 'RH13' reconstruction presented in Raitzsch and Hönisch (2013) not shown here calculated using an apparent fractionation factor. Here we use the Klochko et al. (2006) derived δ 11 B sw value to be consistent with the δ 11 B calculations in this study. 66% confidence intervals are plotted. Dashed lines represent δ 11 B sw -G17 estimates from >17 Ma to highlight uncertainties associated with the δ 11 B sw reconstruction. Pre-industrial and 2017 levels are highlighted by blue and black asterisks respectively on the pH axis. Pleistocene and Eocene pH estimates are shown as vertical bars on the right side of panel E (Hönisch et al., 2009;Foster, 2008;Anagnostou et al., 2016). (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.) mentary Fig. S2). SST estimates were calculated in a similar manner for previously published Mg/Ca records (Hönisch et al., 2009;Seki et al., 2010;Bartoli et al., 2011;Foster et al., 2012;Badger et al., 2013;Greenop et al., 2014;Martínez-Botí et al., 2015a;Chalk et al., 2017). The following equations were used to correct for Mg/Ca seawater and estimate temperature (in units of • C) from foraminiferal Mg/Ca (Mg/Ca test ): where equation (4) is the generic planktic foraminifera calibration of Anand et al. (2003) with a correction applied to the seawater Mg/Ca ratio to account for past Mg/Ca changes. This correction, given in equation (5) (Barker et al., 2003). In all cases Al/Ca was <100 μmol/mol and in the majority of cases Mn/Ca was <1000 μmol/mol indicating cleaning was efficient. Some Mn/Ca values were higher than our 1000 μmol/mol threshold, for example at ODP site 926, but in all cases they do not covary with Mg/Ca or δ 11 B ( Supplementary   Fig. S3). The ODP site 761 Mg/Ca presented in this study is of too low resolution to judge the influence of contamination through correlation of Mg/Ca and Mn/Ca so instead we look to benthic foraminiferal Mg/Ca data from the same site which shows a lack of correlation between Mn/Ca and Mg/Ca (see Lear et al., 2015 for further discussion). ODP site 872 shows no positive correlation between Mn/Ca and Mg/Ca.

Solving for pH and the carbonate system
To convert our reconstruction of surface δ 11 B borate into pH, CO 2 and carbonate saturation state requires three steps: (1) calculate the acid/base chemistry equilibrium constants based on independently known temperature, salinity and seawater major ion composition; (2) calculate surface water pH using the appropriate equilibrium constants, reconstructed δ 11 B borate and independently estimated seawater δ 11 B; and (3) calculate surface water carbonate chemistry based on reconstructed pH and an independently known carbonate system parameter. The uncertainties associated with the various input parameters need to be propagated through the full set of calculations. The input parameters are presented in Fig. 3 and below we describe these calculations in more detail. All required acid/base chemistry equilibrium constants, including the boric acid dissociation constant pK * B , are calculated using the MyAMI model (Hain et al., 2015) based on reconstructed SST, sea surface salinity (SSS), water depth, and seawater major ion composition. We note that the calcium-to-bicarbonate interaction   Table 1) to He and Morse (1993;β (0) Ca-HCO 3 = 0.2 and β (1) Ca-HCO 3 = 0.3, see their Table 5) to bring MyAMI into agreement with experimental data of Pytkowicz and Hawley (1974). SST estimates are derived from Mg/Ca ratios from T. trilobus for ODP Sites 761, 926, 872, 1000A, and Malta from this study and previously published records Badger et al., 2013;Greenop et al., 2014). SST estimates for ODP Site 999 (Seki et al., 2010) are based on a published record derived from the U K 37 proxy. In our CO 2 calculations we assign an uncertainty in SST of ±1 • C (1σ ), which is propagated into  To calculate surface pH from our boron isotope data we follow previous work (e.g., Vengosh et al., 1991;Foster and Rae, 2016) using eqns. (6) and (7): where δ 11 B is the ratio of 10 B and 11 B expressed in delta notation (in per mil, h), 11 B/ 10 B NIST951 is the isotopic ratio of the NIST SRM 951 boric acid standard ( 11 B/ 10 B = 4.04367) (Catanzaro et al., 1970). The boric acid dissociation constant pK * B is calculated as described above. The boron isotope fractionation factor α B is taken to be 1.0272 (Klochko et al., 2006). The isotopic composition of borate, δ 11 B borate , is based on our T. trilobus data and calibration equation (3) and previously published planktic δ 11 B data and specified calibration equations (Supplementary Table S3). As further described in section 3.2 below, we carry out all calculations for three distinct scenarios for secular change in seawater boron isotopic composition (δ 11 B sw ) using the reconstructions by (a) Lemarchand et al. (2002) ("L02"), (b) Raitzsch and Hönisch (2013) ("RH13") and (c) Greenop et al. (2017) ("G17"). Each record is used with its respective reported uncertainty assuming normal distribution. To calculate CO 2 and carbonate saturation state (Ω CaCO 3 ) from our pH reconstruction requires independent knowledge of ambient dissolved inorganic carbon concentration DIC, or any other carbonate system parameter: The different K* terms represent the equilibrium constants governing carbonate chemistry, which we can calculate using the MyAMI model as described above. The [Ca 2+ ] concentration is internally consistent with the MyAMI calculations for both seawater major ion composition scenarios described above (i.e., constant modern or based on Horita et al., 2002 andBrennan et al., 2013). The dissolved inorganic carbon concentration of seawater, DIC, is not generally known through Earth history, and given residence time constraints it may have experienced long-term changes over our study interval (Ridgwell, 2005). As estimating past DIC change in detail is beyond the scope of this study here we devise a datadriven approach to account, to the first order, for long-term drift in DIC in lieu of simply presuming constant DIC across the study interval. To this end we use published benthic foraminiferal δ 11 B data (Raitzsch and Hönisch, 2013;Greenop et al., 2017) to generate a deep ocean pH record using equation (7), using the same three scenarios for δ 11 B sw (i.e. L02, RH13, G17) that we also use to generate our surface ocean pH records (Lemarchand et al., 2002;Raitzsch and Hönisch, 2013;Greenop et al., 2017). We then combine deep ocean pH with estimates of deep ocean carbonate ion concentration ([CO 3 2− ]) to solve for deep ocean DIC. Deep ocean [CO 3 2− ] is calculated based on past changes in the CCD (using Sime et al., 2007 andPälike et al., 2012 as two scenarios) in conjunction with the [Ca 2+ ] defined by the two scenarios for seawater major ion changes (i.e., constant modern seawater composition versus trends based on data by Horita et al., 2002 andBrennan et al., 2013). We define deep ocean carbonate ion concentration as: where K * CaCO 3 is calculated using the MyAMI model for the hydrostatic pressure at the reconstructed depth of the CCD and the calcite saturation index Ω calcite assumes a nominal value of 1.
Hence, the implications of past changes in the CCD for carbonate chemistry are represented via the strong pressure dependence of K * CaCO 3 , which increases by ∼16% per kilometre water depth. We then combine deep ocean carbonate ion concentration and the benthic foraminifera-based pH record to yield estimated DIC: Our data-driven approach makes the assumption that there is no difference in the carbonate ion concentration between the CCD and the core sites of the deep ocean pH reconstruction. This treatment is a significant simplification and certain to introduce error in the DIC estimate. However, we note that observed modern ocean deep ocean gradients in carbonate chemistry are relatively small when compared for example to the depth dependence of K * CaCO 3 that we include to represent past CCD changes. Furthermore, by assuming a nominal Ω calcite of 1 at the pressure of the CCD we effectively introduce systematic bias towards greater estimated DIC because bottom waters at the CCD exhibit moderate undersaturation (Ω < 1) so as to facilitate the near complete dissolution of calcite on the seafloor. With these caveats in mind we note that under all scenarios the DIC estimated from the youngest data is lower than deep ocean DIC in the modern ocean (Supplementary Fig. 4), in particular in the scenarios based on the δ 11 B sw record of Raitzsch and Hönisch (2013). Therefore, we caution that the approach taken here to estimate DIC appears to be biased towards lower absolute values by about 10-20% depending on scenario. However, as is evident from the comparison of the different sensitivity scenarios, the change in estimated DIC over the course of the Neogene is dominantly driven by the counteracting effects of secular [Ca 2+ ] decline and the tendency of deep ocean pH to increase, which act to raise and lower estimated DIC, respectively. While our approach to estimating past DIC has unaccounted bias, the first-order changes in the DIC estimates over time can be attributed to trends in the data-driven input parameters and is reasonable first approach. Therefore, we carry out linear regression on the output of a 10,000 member Monte-Carlo simulation of the deep ocean input data, extracting for each scenario the estimated DIC rate of change and the intercept uncertainty. Based on these regression results, we estimate surface ocean DIC change relative to its known modern concentration for each respective scenario of δ 11 B sw , CCD and seawater major ion composition. That is, we set surface DIC to its modern reference value of 2000 μmol/kg plus the regressed change over time plus a nominal uncertainty based on the regression intercept value ( Supplementary Fig. S5) (i.e., 1σ of 190 to 270 μmol/kg depending on scenario; see Supplementary  Table S4 for regression results). We note that in the early Miocene interval (>17 Myr) of the Greenop et al. (2017) δ 11 B sw reconstruction large uncertainties exist in the δ 11 B sw value for the oldest samples which cause an edge effect in the smoothing algorithm used in that study (Greenop et al., 2017). Thus, for the regressions using the Greenop et al. (2017) δ 11 B sw scenario we only use benthic δ 11 B data <16 Ma to regress DIC and extrapolated further back in time. This is necessary because this possible bias would have been propagated into the regression of DIC change, thus affecting our reconstruction for the entire record.
We note that applying the DIC change regressed from the benthic data to the surface we implicitly assume a constant DIC gradient ( DIC). An alternative would be to regress deep alkalinity and assume a constant alkalinity gradient ( ALK) instead, but we find that these two alternative approaches yield nearly identical results (see Supplementary Fig. S6). We argue that the surface pH that we reconstruct is the dominant driver of reconstructed pCO 2 and surface carbonate saturation because it is a strong constraint on the ALK-to-DIC difference. Whether DIC or ALK are assumed to be constant therefore makes only a minimal difference for our purpose, and the deviation between these two assumptions falls well within the range of our propagated uncertainties.
DIC estimates vary between δ 11 B sw scenarios across the Neogene: both L02 and G17 show a long-term DIC rise (+20 to +37 μmol/kg per million years, depending on scenario, Supplementary Table S5) in contrast to RH13 scenario, which results in essentially constant DIC (+3 to +11 μmol/kg per million years, depending on CCD scenario, Supplementary Table S5) throughout the study period. That is, the scenarios based on the three available δ 11 B sw reconstructions collectively support a modest rate of longterm DIC increase over the study period, which integrated over the duration of the Neogene could amount to a significant expansion of the active carbon pool of the climate system. Along these lines, our assertion of significant Neogene DIC pool expansion agrees with the previously published model results from Tyrrell and Zeebe (2004), Ridgwell (2005), and Caves et al. (2016) (Supplementary  Fig. S7).
The approach taken here is similar to Foster et al. (2012) which used a simple model anchored by paleo-pH depth profile and CCD reconstructions to estimate surface alkalinity and calculate pH and pCO 2 . Using this approach, Foster et al. (2012) estimate surface alkalinity at 1293 ± 200 umol/kg at a single data point (12.72 Ma) and extrapolate across the mid-Miocene. The resultant average mid-Miocene DIC ∼1100 ± 180 umol/kg are much lower than estimates derived in this study ( Supplementary Fig. S7). Building on this, Greenop et al. (2014) adopted the approach of Foster et al. (2012) but increased the uncertainty to ±300 umol/kg to account for likely rapid CCD changes not resolved in available CCD reconstructions. Based on our new data-driven analysis and existing model based estimates (Ridgwell, 2005;Tyrrell and Zeebe, 2004) we argue that the approach by Foster et al. (2012) yielded DIC that was too low. Consequently, due to the higher DIC estimates applied here, the Mid-Miocene pCO 2 estimates calculated from planktic δ 11 B data of Foster et al. (2012) and Greenop et al. (2014), as well as from our new data from ODP site 761 and 926, are significantly higher than previously published. While highlighting again the caveats in our estimates of surface DIC we argue that this approach is an improvement in deriving long-term DIC change for the boron isotope-based reconstruction of pCO 2 and carbonate saturation across the Neogene.

δ 11 B record of T. trilobus from the Neogene
Our δ 11 B record displays a long-term trend in which low δ 11 B is associated with high Mg/Ca in the same planktic foraminifera, and with low benthic oxygen isotope ratios (δ 18 O). All else being equal, this suggests that low pH (and thus high CO 2 ) is associated with high sea surface temperatures and some combination of relatively warm bottom waters and/or low global ice volume (Fig. 2). Reasonable agreement between δ 11 B from geographically diverse sites suggests the record mainly reflects global carbon cycle variations rather than complications relating to diagenesis and variations in the extent of surface water (dis)equilibrium with respect to atmospheric CO 2 (Fig. 1). The broad correlations evident in Fig. 2 alone suggest a tight coupling between the carbon cycle (that predominantly drives δ 11 B c change) and climate (that predominantly drives δ 18 O b change) through the Neogene. Our new δ 11 B data presented in this study shows similar variability to previously published δ 11 B records from the Miocene Badger et al., 2013;Greenop et al., 2014) and Plio-Pleistocene (Hönisch et al., 2009;Seki et al., 2010;Bartoli et al., 2011;Martínez-Botí et al., 2015a;Chalk et al., 2017) and fills in data gaps in the early Miocene and late Miocene ( Fig. 2; Supplementary Fig. S1). From here on we compile δ 11 B estimates from this study alongside previously published δ 11 B records (with appropriate analytical offsets as required, see above) and calculate carbonate system variables to examine Neogene changes in the carbon cycle.

Scenarios and sensitivity testing
The calculations carried out in order to reconstruct seawater carbonate chemistry and atmospheric CO 2 require independent knowledge of the histories of seawater boron isotope composition (δ 11 B sw ) and major ion composition (especially [Ca 2+ ] and [Mg 2+ ]) as well as the CCD. For each of these parameters we make at least two different choices and carry out the full set of calculations, including Monte-Carlo uncertainty propagation, for all twelve combinations. In this section we describe all scenarios and explain the agreements/disagreements between them. The output data for all scenarios is provided as a tabulated dataset in the supplemental material (Supplementary Table S6). In the section below, we show that the uncertainties in our reconstructions are dominated by the δ 11 B sw reconstructions, not the CCD, and hence only consider these scenarios in our main discussion (see section 3.3).  Fig. 2D). The divergence in the early Miocene (>17 Myr) and late Miocene amongst the three δ 11 B sw scenarios yields divergent pH and pCO 2 reconstructions for these time intervals but has only minor impact on our reconstruction of surface Ω arag (Fig. 4). The G17 reconstruction has likely too low δ 11 B sw before 17 Myr due an edge effect (Greenop et al., 2017) and thereby yields anomalously high pH and low pCO 2 for this interval. Conversely, the RH13 reconstruction agrees well with the G17 record from the middle Miocene to present but is likely too high in the early Miocene as these authors assume a linear increase in surface pH through time (Fig. 2D). Further, because Raitzsch and Hönisch (2013) presuppose knowledge of pH in constructions ("δ 11 B sw -L02"; "δ 11 B sw -RH13"; "δ 11 B sw -G17"; Greenop et al., 2017;Raitzsch andLemarchand et al., 2002) and 'Pälike' and 'fluid inclusions' scenarios. 66% confidence intervals are plotted. Dashed lines represent δ 11 B sw -G17 estimates from >17 Ma to highlight uncertainties associated with the δ 11 B sw reconstruction. DIC estimates from the model output of Ridgwell (2005) and Caves et al. (2016) are plotted in panel C. Note CO 2 , Ω arag and ALK were calculated based on surface pH reconstruction and long-term DIC trends. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.) their approach, caution should be used when interpreting reconstructions of pH and pCO 2 that are based on their δ 11 B sw record (Greenop et al., 2017). Both measurement-based records, RH13 and G17, fall within the uncertainty envelope of L02 modelled reconstruction. However, L02 does not show the same multi-million year variability as RH13 and G17 across the Late Miocene, suggesting that the ocean's input and output of boron that the Lemarchand et al. (2002) model is based on requires further refinement. Overall, we are encouraged by the good agreement between the three independent δ 11 B sw records from the Middle Miocene to present, which builds confidence in our reconstructions, but we caution that further refinement is needed in δ 11 B sw of the early and late Miocene. Carbonate system data for all δ 11 B sw scenarios are provided in the output table (Supplementary Table S6) and discussed further below (see Section 3.3). As the 'δ 11 B sw − G17' > 17 Myr reconstruction is biased by the edge effect we use a dashed line to represent this portion of the pH, CO 2 and Ω arag reconstructions.
As for the history of CCD change, which factors into our estimation of Neogene DIC change, comparing equivalent scenarios with contrasting CCD reconstruction (either Pälike et al., 2012or Sime et al., 2007 we find only small differences in estimated DIC change ( Supplementary Fig. S8), with the scenarios based on the Sime et al. (2007) record consistently resulting in slightly greater rates of long-term DIC increase because of the tendency of longterm CCD deepening in that record (see Supplementary Table S4; Supplementary Fig. S5). This difference between CCD scenarios has only very modest implications for ocean carbon chemistry, consistent with earlier findings (Tyrrell and Zeebe, 2004;Ridgwell, 2005). Consequently, regardless of which CCD reconstruction is used, we arrive at essentially identical reconstructions of CO 2 and Ω arag ( Supplementary Fig. S8). This is because the long-term trend in CCD over our study interval (i.e., <500 m) implies only a <8% change in the deep ocean calcium-to-carbonate ion product, DIC, CO 2 and surface saturation state, whereas long-term fractional changes in [Ca 2+ ], benthic [H + ] and surface [H + ] are all much greater. Due to the minor impact the CCD reconstruction has on our results we focus our discussion below on only the scenarios using the CCD reconstruction of Pälike et al. (2012), referred hereafter simply as 'Pälike', because it incorporates data from an array of Pacific DSDP and ODP sediment cores, has better age control, and is of higher resolution.
For our consideration of seawater major ion composition, we carry out our calculations using [Ca 2+ ] and [Mg 2+ ] reconstructions change of Horita et al. (2002) and Brennan et al. (2013) as well as the hypothetical scenario that both concentrations remained unchanged at their modern seawater value throughout the entire study interval (Fig. 3). The reconstructed [Ca] sw derived from fluid inclusions are consistent with calcium isotope model-data comparisons suggesting that [Ca] sw likely decreased through the Neogene (Horita et al., 2002;Heuser et al., 2005;Brennan et al., 2013). Consistent with Hain et al. (2015) we find that the leading order impact of seawater major ion composition operates through [Ca 2+ ] being a direct factor in the CaCO 3 solubility product, not through ion-ion interactions and ion-complex formation affecting the equilibrium constants. As described above, using reconstructed deep ocean pH and CCD (the latter as a constraint on the deep ocean calcium-to-carbonate ion product) the reconstructed ∼40% secular decline in seawater [Ca 2+ ] over the course of the Neogene translates to a reciprocal secular increase in DIC (Equation (10), Supplementary Fig. S4). Conversely, the general tendency of reconstructed deep ocean pH increase imparts a secular decline in estimated DIC because at higher pH less DIC is needed to maintain a given carbonate ion concentration. That is, in the scenarios that include seawater major ion change the opposing effects of [Ca 2+ ] and benthic pH on estimated DIC partially cancel out, leaving a modest residual long-term DIC rise (Supplementary Fig. S9 Fig. S7; Tyrrell and Zeebe, 2004;Ridgwell, 2005;Caves et al., 2016). As a point of comparison, the hypothetical scenario of time-invariant seawater major ion composition eliminates the calcium-driven tendency of secular increase in DIC, yielding instead long-term estimated DIC decline driven by reconstructed long-term benthic pH rise. With regard to our reconstruction of surface carbon chemistry, which is the product of both our extensive surface pH reconstruction and the long-term trends in DIC that we infer from deep ocean data, the opposing trends in DIC change between the "fluid inclusions" and "constant" scenarios impart opposing tendencies upon reconstructed atmospheric CO 2 , with the more realistic "fluid inclusions" scenario exhibiting a reduced secular CO 2 decline ( Supplementary Fig. S9). In contrast, our reconstruction of surface carbonate saturation state is virtually insensitive to these scenarios of changing seawater major ion composition (including calcium concentration) (Supplementary Fig. S9). This is largely because our estimate of DIC change is anchored by deep ocean carbonate saturation state via the CCD, which already accounts for seawater major ion composition changes and the strong effect of [Ca 2+ ] change on the carbonate saturation product. In the discussion below we focus exclusively on the "fluid inclusions" scenario because the "constant" scenario is intended only as a sensitivity test and to illustrate our approach to the calculations.

Neogene pH, CO 2 , and Ω arag reconstructions
The calculation of pH and pCO 2 requires records of seawater chemistry (e.g. δ 11 B sw , [Ca] sw ) and a second carbonate system parameter (here derived via the CCD) which currently have considerable uncertainties. Based on our sensitivity analysis (Section 3.2), we will discuss three scenarios which uses the L02, and RH13, G17 δ 11 B sw reconstructions along with 'Pälike' CCD, and fluid inclusion [Ca 2+ ] and [Mg 2+ ]. Late Miocene pH (∼8.1) are in agreement with earlier model-based estimates (Ridgwell, 2005). The δ 11 B-derived CO 2 record from each scenario are in relatively good agreement with published records from a range of proxies (e.g., alkenones, stomata, paleosols) (Zhang et al., 2014;Foster et al., 2017; Supplementary Fig. S10). However, our record shows more structure, with peaks in CO 2 during the MMCO, late Miocene, and late Pliocene ( Fig. 4; Supplementary Fig. S10). A particular strength of the δ 11 B proxy and our approach is that we also capture simultaneous secular changes in seawater carbonate chemistry.
The most acidic surface conditions (pH∼7.6-7.7) occurred in the MMCO where inter-sample variability likely reflects orbital timescale (∼100-kyr) cycles associated with large (∼200 ppm) CO 2 variations reaching >780-1100 ppm (i.e. max CO 2 value: 891 ±  Greenop et al., 2014). Our MMCO surfacewater pH estimates are similar to the low values reported for the late Eocene (7.6-7.8; Anagnostou et al., 2016), a time of global warmth, high CO 2 and low ice volume (Lear et al., 2000). Although the rest of the record is at lower resolution and contains sample aliasing of orbital-scale variability, some long-term trends are clearly discernible. Following the MMCO, CO 2 decreased and surface pH increased to ∼8.1 at 12.5 Ma, similar to pre-industrial (Figs. 2, 4). Surface pH was ∼8.0-8.1 between 9 and 3 Myr and subsequently increased towards the late Pleistocene glacial/interglacial range of ∼8.3 to ∼8.1 (Foster, 2008;Hönisch et al., 2009;Henehan et al., 2013;Chalk et al., 2017). Unlike the oscillations in pH and CO 2 , our record of low latitude aragonite saturation (Ω arag ) is broadly constant being in the range 3 to 5 for much of the last 24 million years. In detail however mean low latitude Ω arag increased from the MMCO  (Zachos et al., 2008) and (B) CO 2 records from this study using δ 11 B sw -L02, δ 11 B sw -RH13, and δ 11 B sw -G17 reconstructions with 'Pälike' and 'fluid inclusions scenarios; 66% confidence intervals are plotted.
It is worth noting that the rather large uncertainty intervals we report for each of a number of scenarios are the result of conservative propagation of uncertainty in the various input parameters.
Only a small proportion of that uncertainty is due to the δ 11 B measurements and the other proxy data required for the boron isotope pH proxy. Instead, structural uncertainty in seawater composition -major ion, δ 11 B sw , DIC -are compounded in our approach to reconstruct seawater carbonate chemistry over geologic time-scales much longer than the residence time of each of these constituents. To illustrate that point we consider an alternative scenario for the Plio-Pleistocene ( Supplementary Fig. S11), a geologically recent interval where seawater major ion, δ 11 B sw , DIC can be anchored at modern values and their structural uncertainty is usefully assessed through the spread in the respective scenarios. For this alternative scenario we propagate the uncertainties due to the spread of δ 11 B sw reconstructions and the spread in our scenarios of longterm DIC drift under the 'fluid inclusion' major ion scenario, while maintaining the propagation of uncertainties in temperature, salinity and foraminiferal δ 11 B the same. In that scenario we find uncertainty intervals for pH, CO 2 and carbonate saturation states small enough to confidently resolve glacial/interglacial cycles, something not possible with our treatment of the full Neogene dataset (Supplementary Fig. S11; Supplementary Table S7). That is, while our approach to reconstruct seawater carbonate chemistry of the Neogene is able to resolve long-term trends and large magnitude episodes such as the MMCO there remains considerable room for improving the reconstructions by reducing uncertainties related to seawater composition.

Evaluating Miocene CO 2 -climate decoupling
Long-term climate trends across the Neogene, as inferred from benthic foraminiferal δ 18 O b records, show a transition from an extended period of global warmth of the MMCO to the modern icehouse world marked by the Middle Miocene Climate Transition (MMCT) at 13.9 Ma where Antarctic ice sheets expanded and climate cooled (Fig. 5;Zachos et al., 2001Zachos et al., , 2008Cramer et al., 2009;Lear et al., 2010). Following this global cooling step, the lack of trend in the δ 18 O b record suggests that Earth's climate was relatively stable during the late Miocene (11.6-5.3 Ma) (Zachos et al., 2001(Zachos et al., , 2008. However, a recent global SST compilation indicates a long-term surface cooling across this interval reaching near modern ocean temperatures between 5 and 7.4 Myrs ago (Herbert et al., 2016). In contrast to the late Miocene, the early Pliocene was a period of prolonged warmth terminated at 2.7 Ma with the initiation of northern hemisphere glaciation and subsequent intensification of glacial interglacial climate variability in the Pleistocene (Lisiecki and Raymo, 2005). Previous CO 2 records suggest that the transition from the warm MMCO to the cooler post-MMCT was associated with a decline in CO 2 (Kürschner et al., 2008;Foster et al., 2012). The warmth of the early Pliocene has also been previously attributed to higher than pre-industrial atmospheric CO 2 levels (Seki et al., 2010;Bartoli et al., 2011;Martínez-Botí et al., 2015a). However, the role of CO 2 in evolution of global climate across the Late Miocene is still under debate.
Alkenone-derived CO 2 records from Pagani et al. (1999aPagani et al. ( , 1999b) exhibited low CO 2 levels during the early to mid Miocene with no significant shifts in the late Miocene. This led the authors to identify no long-term correlation of climate and CO 2 and suggest a late Miocene decoupling between CO 2 and climate (Pagani et al., 1999a(Pagani et al., , 1999b(Pagani et al., , 2005. However, subsequent work has shown that foraminifera-δ 18 O based SST estimates used in these late Miocene alkenone-CO 2 calculations were too low due to inaccuracies in their estimates of ice volume, calcification depth, and diagenetic overprint (Pagani et al., 2010). Subsequent TEX 86 derived SST estimates indicate relatively higher temperatures that would suggest higher CO 2 levels than previously published, but still supporting a decoupling between climate and CO 2 in the late Miocene (Zhang et al., 2014).
Late Miocene climate records suggesting warmer than modern continental and oceanic temperatures and the relatively low alkenone-based CO 2 levels are therefore difficult to reconcile with a dominant role of CO 2 in setting the evolution of late Miocene climate (Pound et al., 2011;LaRiviere et al., 2012). Specifically, LaRiviere et al. (2012) shows that North Pacific SSTs were significantly higher than today and suggest that in the absence of significantly elevated CO 2 , a complex mechanism involving a deepened thermocline, expansion of tropical warm pool and increased water vapour may have maintained the late Miocene warmth. Here we use our new Neogene CO 2 records to examine the degree of coherency between CO 2 and δ 18 O b and revisit the idea of climate-CO 2 decoupling in the Miocene.
Cross-plots between CO 2 climate forcing (i.e., a linear function of the logarithm of CO 2 , i.e., ln(CO 2 /278); Myhre et al., 1998) and δ 18 O b , for all three δ 11 B sw scenarios, show a correlation across the  and Pliocene (0-5 Ma) (Fig. 6). The good correspondence between CO 2 and δ 18 O b across these time intervals suggests a dominant role for CO 2 in forcing climate. Despite this good coupling between global climate (as recorded by benthic δ 18 O) and CO 2 for the Plio-Pleistocene and Mid Miocene, across the late Miocene (5.3-11.6 Ma) we record CO 2 changes that are not strongly matched by corresponding changes in the δ 18 O b record (i.e. no variance in δ 18 O b that corresponds to the variance in our reconstructed CO 2 ; R 2 = 0.02-0.09; Figs. 5; 6). Regardless of the chosen δ 11 B sw scenario, all our CO 2 reconstructions shows CO 2  (Zachos et al., 2008). CO 2 reconstruction used in this comparison is from all δ 11 B sw scenarios (G17, RH13, L02) and 'Pälike' and 'fluid inclusions' reconstructions. peak at ∼9 Ma that is not seen in the climate records (Fig. 5). This peak in CO 2 at ∼9 Ma is also seen in recent cell-sized corrected alkenones and pennate diatoms (diffusive and ACTI-CO models) CO 2 records (Bolton et al., 2016;Mejía et al., 2017;Supplementary Fig. S12). However, quantification of the CO 2 variations using the diatom method awaits determination of a diatom active carbon uptake model. Nonetheless, the convergence in CO 2 trends from 11.6 to 8.5 Ma and the mismatch between CO 2 and climate trends support the notion of a late Miocene "climate-CO 2 decoupling" (Pagani et al., 1999a(Pagani et al., , 1999b. This could reflect the significance of oceanographic or palaeogeographic changes and related feedbacks in setting the late Miocene climate, as suggested by others (Pagani et al., 1999a(Pagani et al., , 1999bFoster et al., 2010;LaRiviere et al., 2012), which may have masked or counteracted changes in the CO 2 forcing at this time.
However, there remain two principle caveats that prevent definite verification of the apparent late Miocene decoupling of CO 2 and climate with our current CO 2 record: (1) while we present the highest resolution late Miocene CO 2 reconstruction to date, our record still does not resolve CO 2 change at an orbital timescale and may thus contain aliasing of orbital scale CO 2 variability, and (2) during the late Miocene there is growing evidence for rapid and large magnitude δ 11 B sw change, evident in the two measurementbased δ 11 B sw records (RH13 and G17) (Fig. 2). Boron isotope based reconstructions of carbon cycle and CO 2 change during this period critically depends on reconstructions that correctly capture both the pace and timing of this δ 11 B sw shift. Interestingly, there is a correlation of our planktic δ 11 B to δ 18 O b (R 2 = 0.22) that is not evident in CO 2 calculated from the boron isotope data (Supplementary Fig. S13). This observation could be used to argue that some intervals of late Miocene carbon cycle change may in fact have changed in tune with climate but this relationship is currently obscured by the uncertainties in the parameters needed to convert δ 11 B to CO 2 (and pH). Therefore, a detailed and robust record of δ 11 B sw change across the late Miocene is needed to fully elucidate the nature of CO 2 -climate coupling at this time.

Geological context of future ocean acidification
Over the last 250 yrs, anthropogenic emissions have increased atmospheric CO 2 levels by >40% to ∼400 ppm (Keeling et al., 2001;Betts et al., 2016; https://www.co2 .earth) more than 120 ppm above the range of glacial-interglacial cycles of at least the last 800,000 yrs (Petit et al., 1999;Siegenthaler et al., 2005;Luethi et al., 2008;Bereiter et al., 2015). The oceans have absorbed about a third of the CO 2 emitted since the 1750s (Doney et al., 2009;Le Quere et al., 2016), with the CO 2 titrating ocean carbonate ion so as to simultaneously decrease pH and calcium carbonate saturation state. Predictions of future pH and saturation state are widely regarded as sufficient to disrupt physiological processes and reduce calcification, development, and fecundity in a range of marine organisms (Gattuso et al., 2015). Furthermore, reduced pH will impact the availability of essential nutrients (i.e. N) and trace metals (i.e. Fe) critical to marine primary producers for biological functions, altering marine biogeochemical cycles (Hutchins et al., 2007;Hoffmann et al., 2012).
In Fig. 7 we compare our geological record of pH and CO 2 to future predictions based on IPCC Representative Concentration Pathway (RCP) scenarios from the stringent (RCP2.6) to the highcarbon business as usual trajectory (BAU; RCP8.5) which extend to 2100 (Bopp et al., 2013). Note we limit our discussion here to MC-ICPMS based pH, CO 2 , and Ω arag dataset (Seki et al., 2010;Foster et al., 2012;Badger et al., 2013;Greenop et al., 2014;Martínez-Botí et al., 2015a;Chalk et al., 2017; this study), to maintain internal consistency free of poorly defined analytical offsets (see previous discussion). Because ocean acidification will continue beyond 2100, we also make the comparison with estimates up to the year 2500 using emission pathways releasing 250 PgC and 2000 PgC (the near equivalent to RCP2.6 and BAU RCP8.5, respectively), 5000 PgC (a possible upper limit including unconventional fossil fuels such as shale gas), and 10,000 PgC (a possible upper limit on all available carbon resources) (Winkelmann et al., 2015;Bruckner et al., 2014). To compare our low latitude specific Ω arag to the global average as presented in Winkelmann et al. (2015) model, we adjust model Ω arag by a constant offset (i.e. 0.9) based on the pre-industrial difference between low latitude and global model estimates. Current CO 2 reconstructions of 400 ppm (Keeling et al., 2001;Betts et al., 2016; https://www.co2 .earth) are similar to the highest values we reconstruct during the interglacials of the Pliocene and in late Miocene. By 2100, pH levels on the RCP8.5 pathway will be <7.8 and CO 2 near 930 ppm, unprecedented for the last 14 million years (66% CI). If continued unabated (following the 2000 PgC emission pathway) conditions will reach a pH as low as 7.7 by 2150 and CO 2 as high 960 ppm, equivalent to some of the lowest MMCO levels (66% CI). Exhaustion of conventional and unconventional (10,000 PgC) fossil fuel reserves could reduce pH to 7.2 and increase CO 2 to ∼4200 ppm by 2300 (Bruckner et al., 2014;Winkelmann et al., 2015), unprecedented in at least 22 Myr (at 66% confidence). By 2100 Ω arag levels will be <2.75 in the BAU scenarios, lower than any point in the last 8 million years (66% CI) in two scenarios. Further CO 2 emissions predict Ω arag to reach globally undersaturated conditions. The only parallels of such extreme and geologically rapid acidification in the geological past is associated with major mass extinction (Clarkson et al., 2015).
On geologic timescales (>10,000 yrs), the balance between the delivery and removal of alkalinity to the oceans via carbonate weathering and burial acts to partially decouple pH and carbonate saturation state ( Fig. 4; Hönisch et al., 2012). Previous studies suggest that surface ocean saturation state was relatively constant across the last 65 Myrs excluding short-term events such as Eocene-Oligocene Transition and Paleocene-Eocene Thermal Maximum (Tyrrell and Zeebe, 2004;Ridgwell and Zeebe, 2005;Hönisch et al., 2012;Pälike et al., 2012;Zeebe, 2012;Penman et al., 2014). However, these Cenozoic records of Ω rely on carbon cycle models that make significant assumptions about CO 2 across the time interval. Our record of the ocean carbonate system shows broadly constant values for the Neogene in the range 3-5, but with clear minimum Ω arag values, lower than pre-industrial levels, during the warmth of the MMCO when pH values are the lowest in our record (Figs. 2, 4). Pacific % carbonate records (Holbourn et al., 2013) indicate that there were large changes in the CCD on glacial-interglacial timescales during the MMCO indicating a dynamic carbon cycle. The fossil record shows significant global and regional extinction of corals and molluscs in the early and middle Miocene, although arguably not an extraordinary occurrence, while calcifying plankton thrived during the MMCO (Budd, 2000;Warren and Bottier, 2001;Bown et al., 2004;Ezard et al., 2011). Considering the possible consequences of ongoing acidification for marine organisms, however, one must balance these observations with the fact that geological rates of change were much slower leading into the MMCO than current anthropogenic change, allowing for more time for species to adapt (e.g. scleractinian reefs in the middle Miocene; Perrin, 2002;Wilson, 2008;Hönisch et al., 2012).
The existence of past warm climates, combined with a history of glacial-interglacial variability suggests that, if given time, many calcifiers may be able to adapt to slightly more acidic and less saturated waters than occur at present. However there is no clear analogue for the rate and magnitude of change over the coming centuries under BAU emissions scenarios. End-of-century BAU pH predictions are comparable to the long-term pH regime (pH < 7.8) experienced in natural CO 2 seep sites, suggesting future, acclimatized tropical ecosystems may have a lower species diversity with higher number of acidification tolerant species (Hofmann et al., 2011;Fabricius et al., 2011). Reduced pH will alter the biogeochemical cycling and the ocean inventory on nutrients and trace metals, which can significantly impact marine biota function under future CO 2 scenarios (Hutchins et al., 2007;Hoffman et al., 2012). Moreover, ocean acidification coincides with warming and other environmental stresses, potentially exacerbating these negative effects (Gattuso et al., 2015). The current pH of the surface ocean is already probably lower than at any time in the last 2 Myr and the rate of acidification is increasing (Hönisch et al., 2009;IPCC, 2014). Some organisms are already negatively impacted (Bednarsek et al., 2012). The predicted changes in pH and saturation state constitute a substantial threat to marine life, especially for the higher emissions scenarios in which the chemistry of the open ocean is forced, with unprecedented rapidity, into a state that is clearly well be-  (Foster, 2008;Hönisch et al., 2009). Highlighted on the right y-axis is percent difference in [H + ] from preindustrial levels. B) CO 2 , plotted on the log scale reconstruction with Monte Carlo uncertainty estimates (66% CI). Ice core CO 2 data for the last 0.8 Myr are plotted in blue (Petit et al., 1999;Siegenthaler et al., 2005;Luethi et al., 2008;Bereiter et al., 2015).

Conclusions
Here we present a Neogene record of planktic foraminiferal δ 11 B from which we estimate surface water pH and Ω arag as well as atmospheric CO 2 . These calculations also incorporate previous reconstructions of seawater δ 11 B (δ 11 B sw ), major ion composition and deep sea carbonate burial (i.e., CCD). Our sensitivity analysis demonstrates that different δ 11 B sw reconstructions influence our estimates of both pH and CO 2 , most significantly during the early and late Miocene when δ 11 B sw reconstructions diverge. Additionally, the reconstructed secular decline in seawater calcium concentration implies a secular DIC increase that counteracts a portion of reconstructed pH-rise-driven secular CO 2 decline but has only a weak effect on reconstructed surface aragonite saturation state. The choice of CCD record used in the reconstruction of DIC has mostly insignificant effects on our reconstructions. Our pH records suggests that the lowest pH of the past 22 Myr occurred in the Middle Miocene Climatic Optimum, when pH levels were similar to those of the early Cenozoic. This interval of global warmth is associated with average CO 2 in all our scenarios, reaching >470-630 ppm. Following the MMCO there is a general long-term trend towards lower CO 2 as climate cooled towards the Pleistocene. Superimposed on this trend is an apparent CO 2 peak at ∼9 Ma, reproduced in all scenarios, that is not matched by corresponding changes in SST and δ 18 O b records. However, to definitely resolve the degree climate-carbon cycle coupling of the late Miocene will require more detailed records of the pace and magnitude of δ 11 B sw change during that interval. Finally, we show that it is unlikely that marine ecosystems experienced low latitude surface pH and Ω arag below 7.7 and 1.7 respectively during the past 14 million years (66% CI). Under business as usual scenarios MMCO CO 2 and pH levels will be reached by 2100 and Ω arag will be lower than any point in the last 14 million years (66% CI). Understanding the full implication of this conclusion requires long-term laboratory and field studies examining physiological and ecosystem responses to the full range of biogeochemical changes associated with ocean acidification as well as additional observations from the fossil record.