A record of Neogene seawater δ 11 B reconstructed from paired δ 11 B analyses on benthic and planktic foraminifera

. The boron isotope composition ( δ 11 B) of foraminiferal calcite reﬂects the pH and the boron isotope composition of the seawater the foraminifer grew in. For pH reconstructions, the δ 11 B of seawater must therefore be known, but information on this parameter is limited. Here we reconstruct Neogene seawater δ 11 B based on the δ 11 B difference between paired measurements of planktic and benthic foraminifera and an estimate of the coeval water column pH gradient from their δ 13 C values. Carbon cycle model simulations underscore that the 1 pH– 1δ 13 C relationship is relatively insensitive to ocean and carbon cycle changes, validating our approach. Our reconstructions suggest that δ 11 B sw was ∼ 37.5 ‰ during the early and middle Miocene (

G r e e n o p, Ro s a n n a , H ai n, M a t hi s P., S o s di a n, Si n di a M .ORCID: h t t p s://o r ci d.o r g/ 0 0 0 0-0 0 0 2-4 5 9 9-5 5 2 9, Oliver, Kevin I. C., Goo d wi n, P hilip, C h al k, Tho m a s B., L e ar, C a r oli n e H . ORCID: h t t p s://o r ci d.o r g/ 0 0 0 0-0 0 0 2-7 5 3 3-4 4 3 0, Wilso n, P a ul A. a n d Fo s t er, G avi n L. 2 0 1 7.A r e c o r d of N e o g e n e s e a w a t e r 1 1 B r e c o n s t r u c t e d fro m p ai r e d 1 1 B a n alys e s o n b e n t hi c a n d δ δ pl a n k tic fo r a mi nife r a. Cli m a t e of t h e P a s t 1 3 , p p. 1 4 9-1 7 0. 1 0. 5 1 9 4/c p-1 3-1 4 9-2 0 1 7 file P u blis h e r s p a g e : h t t p s:// doi.o r g/ 1 0. 5 1 9 4/c p-1 3-1 4 9-2 0 1 7 < h t t p s :// doi.o r g/ 1 0. 5 1 9 4/ c p-1 3-1 4 9-2 0 1 7 > Pl e a s e n o t e: C h a n g e s m a d e a s a r e s ul t of p u blis hi n g p r o c e s s e s s u c h a s c o py-e di ti n g, fo r m a t ti n g a n d p a g e n u m b e r s m a y n o t b e r efl e c t e d in t his ve r sio n.Fo r t h e d efi nitiv e ve r sio n of t hi s p u blic a tio n, pl e a s e r ef e r t o t h e p u blis h e d s o u r c e.You a r e a d vis e d t o c o n s ul t t h e p u blis h e r's v e r sio n if yo u wi s h t o cit e t hi s p a p er.
Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s. S e e h t t p://o r c a .cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s.Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

Introduction
Key to determining the relationship between CO 2 and climate in the geological past is the calculation of reliable estimates of absolute CO 2 through time.In recent years the boron isotope composition (δ 11 B) of foraminiferal calcite has become a high-profile tool for reconstructing CO 2 beyond the last 800 ky and throughout the Cenozoic Era (Foster, 2008;Hönisch et al., 2009;Pearson et al., 2009;Bartoli et al., 2009;Foster et al., 2012;Badger et al., 2013;Henehan et al., 2013;Greenop et al., 2014;Martínez-Botí et al., 2015a).Yet long-term change in the boron isotope composition of seawater (δ 11 B sw ) is currently poorly constrained and represents a major source of the uncertainty associated with δ 11 B-determined CO 2 estimates (e.g.Pearson et al., 2009).
In the modern ocean, boron is a conservative element with a spatially invariant isotope ratio (39.61 ‰; Foster et al., 2010), but this value is subject to change through geological time.
The residence time of boron in the ocean is estimated to lie between 11 and 17 My (Lemarchand et al., 2000).Therefore, we can expect the uncertainty associated with δ 11 B sw to be an important factor in CO 2 estimates beyond the late Pliocene (∼ 4-5 Ma; Palmer et al., 1998;Lemarchand et al., 2000;Pearson et al., 2009;Foster et al., 2012;Anagnostou et al. 2016).
The ocean boron budget and its isotopic composition are controlled by a number of inputs and outputs (Fig. 1).How-  Lemarchand et al. (2000) and Park and Schlesinger (2002).Isotopic compositions are from Lemarchand et al. (2000), Foster et al. (2010) and references therein.ever, because the magnitude of the boron fluxes between land, the ocean and the atmosphere in modern times are still poorly understood, the residence time and changes in both concentration ([B] sw ) and isotopic composition (δ 11 B sw ) through time remain uncertain.The main inputs of B into the ocean are silicate weathering, and to a lesser extent evaporite and carbonate weathering, delivered to the ocean by rivers (Lemarchand et al., 2000;Rose et al., 2000;Lemarchand and Gaillardet, 2006), hydrothermal vents (You et al., 1993) and fluid expelled from accretionary prisms (Smith et al., 1995).The major loss terms are low-temperature oceanic crust alteration (Smith et al., 1995), adsorption onto sediments (Spivack and Edmond, 1987) and co-precipitation into carbonates (Hemming and Hanson, 1992).In the case of all three outputs, the light 10 B isotope is preferentially removed relative to 11 B, such that the seawater 11 B / 10 B ratio (δ 11 B sw , 39.61 ‰) is significantly greater than that of the cumulative inputs (δ 11 B of ∼ 10.4 ‰; Lemarchand et al., 2000).Our understanding of the modern boron fluxes outlined above, and illustrated in Fig. 1, implies a significant imbalance between inputs and outputs.Consequently, the poorly constrained ocean-atmosphere boron fluxes may also be an important part of the ocean's modern boron mass balance (Park and Schlesinger, 2002).Here, however, we follow Lemarchand et al. (2000) in assuming that atmospheric fluxes are unlikely to have varied significantly on geological timescales and therefore will not be discussed further in reference to the Neogene δ 11 B sw record we present.
Unlike many other isotope systems (e.g.δ 7 Li sw , δ 26 Mg sw , δ 44/40 Ca sw , 87 Sr / 86 Sr) to date, no direct archive has been documented for δ 11 B sw .This is a result of the pH-dependent boron speciation in seawater upon which the δ 11 B-pH proxy is based (Hemming and Hanson, 1992), which imparts a pH dependency on the δ 11 B of all marine precipitates so far ex-amined.Empirical reconstructions of δ 11 B sw must therefore use "indirect" approaches.So far four approaches have been applied to the problem (Fig. 2): (1) geochemical modelling (Lemarchand et al., 2000), (2) δ 11 B analysis of halites (Paris et al., 2010), (3) measurements of benthic foraminiferal δ 11 B coupled to various assumptions about past changes in ocean pH (Raitzsch and Hönisch, 2013), and (4) measurements of δ 11 B in surface-and thermocline-dwelling foraminifera coupled with additional information on the pH gradient of the surface ocean (Palmer et al., 1998;Pearson andPalmer 1999, 2000;Anagnostou et al., 2016).
Geochemical modelling of the changes in the flux of boron into and out of the ocean through time has been used to suggest that δ 11 B sw increased from 37 ‰ at 60 Ma to 40 ‰ ± 1 ‰ today, driven by a combination of processes including changing boron continental discharge (Lemarchand et al., 2000).In the case of approach 2, while modern natural halites reflect δ 11 B sw (39.7 ‰) with no apparent fractionation, measurement of δ 11 B in ancient halites yields isotopic ratios that are significantly lower than all other approaches (Fig. 2; Paris et al., 2010), with implausible variability among samples of the same age (7 ‰ range), thereby casting doubt over the reliability of this approach (Raitzsch and Hönisch, 2013).In the case of approach 3, δ 11 B sw is calculated from globally distributed benthic δ 11 B data, with an imposed degree of deep-ocean pH change (Fig. 2; Raitzsch and Hönisch, 2013).This method hinges on two key assumptions: (a) a near-linear surface water pH increase of 0.39 over the past 50 My taken from the average pH output from a number of modelling studies (Berner and Kothavala, 2001;Tyrrell and Zeebe, 2004;Ridgwell, 2005) and (b) a prescribed constant surface-to-deep ocean pH gradient of 0.3 (Tyrrell and Zeebe, 2004, and modern observations).The modelled surface pH and estimated fixed pH gradient are then used to estimate Paris et al., 2010Pearson and Palmer, 2000Lemarchand et al., 2002Foster et al., 2012Raitzsch and Honisch, 2013 (Klochko) Raitzsch and Honisch, 2013 (Empirical)

Figure 2.
A compilation of published δ 11 B sw records.Seawater composition reconstructed from foraminifera depth profiles (light blue squares and dark blue crosses) from Pearson and Palmer (2000) and Foster et al. (2012) respectively; numerical modelling (green line), with additional green lines showing ±1 ‰ confidence interval (Lemarchand et al., 2000); benthic δ 11 B (purple diamonds and dark purple line showing 5-point moving average uses the fractionation factor of Klochko et al., 2006, light purple line showing 5-point moving average uses an empirical calibration) from Raitzsch and Hönisch (2013); and halites (orange crosses) from Paris et al. (2010).The orange crosses in brackets were discarded from the original study.
deep ocean pH, and then convert benthic foraminiferal δ 11 B measurements into δ 11 B sw .This approach yields broadly similar results to geochemical modelling (Fig. 2).Approach 4 exploits the non-linear relationship between δ 11 B and pH alongside estimated pH gradients in the ocean to constrain δ 11 B sw (Palmer et al., 1998;Pearson andPalmer 1999, 2000), and this is the basis of the approach used in this study.The advantage of this method is that δ 11 B sw can be reconstructed empirically without relying on a priori absolute-pH constraints.The non-linear relationship between δ 11 B and pH means that the pH difference between two δ 11 B data points varies as a function of δ 11 B sw (Fig. 3).Consequently, if the size of the pH gradient can be estimated, then there is only one δ 11 B sw value that is consistent with the foraminiferal δ 11 B measurements and the specified pH gradient, irrespective of the absolute pH (Fig. 3c).Previously this approach was applied to pH variations in the surface ocean and used in studies of Cenozoic pCO 2 to account for changes in δ 11 B sw (determined using δ 11 B in surface-and thermocline-dwelling foraminifera) (Fig. 2) (Palmer et al., 1998;Pearson andPalmer 1999, 2000;Anagnostou et al., 2016).This approach uses a constant pH gradient between the surface and some depth proximal to the oxygen minimum zone and the boron isotope values of a mixed-layerdwelling species and thermocline dweller to calculate a value for δ 11 B sw (Pearson and Palmer, 1999).The resulting record suggests that δ 11 B sw varies between 37.7 ‰ and 39.4 ‰ throughout the Neogene (Fig. 2) (Pearson and Palmer, 2000).
The same method, but using planktic-benthic instead of surface planktic-thermocline planktic δ 11 B gradients to calculate δ 11 B sw , was recently applied to the middle Miocene where it yielded a δ 11 B sw of 37.6 +0.4  −0.5 ‰ (Foster et al., 2012).A further modification to the method of Pearson and Palmer (1999) was also proposed in that study wherein δ 13 C in foraminiferal calcite was used to estimate the surface-todeep pH gradient (Foster et al., 2012).Here, we reconstruct δ 11 B sw for the last 23 My, the Neogene, based on this modified approach.We undertake extensive sensitivity tests using both the CYCLOPS carbon cycle box model and the GE-NIE Earth system model to define the plausible range in the relationship between surface-deep pH difference and δ 13 C difference, which is an essential parameter for this approach.Finally, we employ a Monte Carlo approach for comprehensive propagation of uncertainty in all input parameters, and we focus on reconstructing δ 11 B sw .The implications of our work for understanding the evolution of Neogene ocean pH and atmospheric pCO 2 will be documented elsewhere.

Site locations and age models
Foraminifera from four sites are used to construct the planktic-benthic δ 11 B pairs; Ocean Drilling Program, ODP, Site 758 and ODP Site 999 for the Pleistocene and Pliocene samples and ODP Site 926 and Site 761 for the Miocene samples (Fig. 4) (this study; Foster et al., 2012;Martìnez-Botì et al., 2015a, and a follow up study by Sosdian et al., 2017).We also incorporate the middle-Miocene planktic-benthic pair from Site 761 in Foster et al. (2012).To place all data from all sites on a single age model we use the nanno and planktic foraminifera stratigraphy from sites 999, 926 and 761 (Shipboard Scientific Party, 1995, 1997;Zeeden et al., 2013;Holbourn et al., 2004) updated to GTS2012 (Gradstein et al., 2012).At Site 758 the magnetostratigraphy (Shipboard Scientific Party, 1989) is used and updated to GTS2012 (Gradstein et al., 2012).

Boron isotope analysis and pH calculation
The boron isotope measurements (expressed in delta notation as δ 11 B -per mil variation) were made relative to the boric acid standard standard reference material (SRM) 951; (Catanzaro et al., 1970).Boron was first separated from the Ca matrix prior to analysis using the boron-specific resin Amberlite IRA743 following Foster et al. (2013).The boron isotope composition was then determined using a samplestandard bracketing routine on a Thermo Fisher Scientific Neptune multicollector inductively coupled plasma mass spectrometer (MC-ICPMS) at the University of Southampton (following Foster et al., 2013).The relationship between δ 11 B of CaCO 3 and pH is very closely approximated by the following equation: where pK * B is the equilibrium constant, dependent on salinity, temperature, pressure and seawater major ion composition (i.e.[Ca] and [Mg]); ∝ B is the fractionation factor be-Clim.Past, 13, 149-170, 2017 www.clim-past.net/13/149/2017/tween the two boron species; and δ 11 B sw is the boron isotope composition of seawater.Here we use the fractionation factor of 1.0272, calculated from spectrophotometric measurements (Klochko et al., 2006).No temperature correction was applied because a number of recent studies suggest that it is not significant over our investigated temperature range (Rae et al. 2011;Henehan et al., 2013;Martínez-Botí et al., 2015b;Kaczmarek et al. 2016).
Although the δ 11 B of foraminifera correlates well with pH and hence [CO 2 ] aq , the δ 11 B calcite is often not exactly equal to δ 11 B borate (Sanyal et al., 2001;Foster, 2008;Henehan et al., 2013).The planktic species used to construct the benthicplanktic pairs changes through time since a single species is not available for the entire Neogene (this study; Foster et al., 2012;Martìnez-Botì et al., 2015a, and a follow up study by Sosdian et al., 2017).Here Globigerinoides ruber is used for 0 to 3 Ma, Trilobatus sacculifer (formally Globigerinoides sacculifer and including Trilobatus trilobus; Hembleden et al., 1987;Spezzaferri et al., 2015) for 0 to 20 Ma and Globigerina praebulloides for 22 to 23 Ma.The calibration for G. ruber (300-355 µm) is derived from culturing data supported by core top data (Henehan et al., 2013).The T. sacculifer calibration (300-355 µm) is from a follow up study by Sosdian et al. (2017) where the T. sacculifer calibration of Sanyal et al. (2001) is used with a modified intercept so that it passes through the core top value for T. sacculifer (300-355 µm) from ODP Hole 999A (Seki et al., 2010).Unlike the nonsymbiotic modern G. bulloides, G. praebulloides appears to be symbiotic at least in the latest Oligocene (Pearson and Wade, 2009).Therefore, we apply the T. sacculifer (300-355 µm) calibration to this species.For T. sacculifer (500-600 µm) between 0 and 1 Ma, we use the calibration from Martìnez-Botì et al. (2015b), where the calibration of Sanyal et al. (2001) measured using negative thermal ionization mass spectrometry (NTIMS) is corrected for the offset between MC-ICPMS and NTIMS using a comparison of core top T. sacculifer measured by the two different methods from adjacent sites (Foster, 2008;Sanyal et al., 1995).In order to constrain deep-water pH, analysis was conducted on benthic foraminifera Cibicidoides wuellerstorfi or Cibicidoides mundulus depending on which species was most abundant in each sample.The δ 11 B of both Cibicidoides species shows no offset from the theoretical δ 11 B of the borate ion and therefore no calibration is needed to adjust for species-specific offsets (Rae et al., 2011).
As mentioned above, in addition to δ 11 B calcite , temperature, salinity, water depth (pressure) and seawater major ion composition are also needed to calculate pH from δ 11 B. We use the MyAMI Specific Ion Interaction Model (Hain et al., 2015) to calculate the appropriate equilibrium constants based on existing [Ca] and [Mg] reconstructions (Horita et al., 2002;Brennan et al., 2013).Sea surface temperature (SST) is calculated from tandem Mg / Ca analyses on an aliquot of the δ 11 B sample (with a conservative 2σ uncertainty of 2 • C).Adjustments were made for changes in Mg / Ca sw using the records of Horita et al. (2002) and Brennan et al. (2013).We corrected for changes in dependence on Mg / Ca sw following Evans and Müller (2012), using H = 0.41 calculated from T. sacculifer (where H describes the power relationship between test Mg / Ca incorporation and Mg / Ca sw ; Delaney et al., 1985;Hasiuk and Lohmann, 2010;Evans and Müller, 2012) 2011), using the calibration of Lear et al. (2010), and applying this change to the modern bottom-water temperature at each site taken from the nearest Global Ocean Data Analysis Project (GLODAP) site (with a conservative 2σ uncertainty of 2 • C).Salinity is held constant at modern values determined from the nearest GLODAP site (2σ uncertainty of 2 ‰ uncertainty) for the entire record.Note that temperature and salinity have little influence on the calculated pH, and the uncertainty in δ 11 B sw is dominated by the uncertainty in the δ 11 B measurement and the estimate of the pH gradient.
The majority of the δ 13 C data were measured at Cardiff University on a Thermo Finnigan MAT 252 coupled with a Kiel III carbonate device for automated sample preparation.Additional samples were measured on a gas source mass spectrometer Europa GEO 20-20, University of Southampton, equipped with an automated carbonate preparation device and a Finnigan MAT 253 gas isotope ratio mass spectrometer connected to a Kiel IV automated carbonate preparation device at the Zentrum für Marine Tropenökologie (ZMT), Bremen.The Pliocene benthic δ 13 C from Site 999 were taken from the nearest sample in Haug and Tiedemann (1998).In almost all cases δ 13 C was analysed on the same foraminiferal species as δ 11 B and Mg / Ca (38/44 samples).Where this was not possible another surface dweller or benthic foraminifera was used from the same depth habitat.C. wuellerstorfi or C. mundulus were measured in all cases for benthic δ 13 C. Stable isotope results are reported relative to the Vienna Pee Dee Belemnite (VPDB) standard.We use a carbon isotope vital effect for G. ruber (+0.94 ‰; Spero et al., 2003), T. sacculifer and G. praebulloides (+0.46 ‰;Spero et al., 2003;Al-Rousan et al., 2004;), C. mundulus (+0.47 ‰;McCorkle et al., 1997), and C. wuellerstorfi (+0.1 ‰;McCorkle et al., 1997) to calculate the δ 13 C of dissolved inorganic carbon (DIC).

Carbon isotopes as a proxy for vertical ocean pH gradient
The use of δ 13 C in foraminiferal calcite to estimate the surface-to-deep pH gradient requires knowledge of the slope of the pH-δ 13 C relationship in the past.In this section we briefly outline the main factors that contribute to the pHδ 13 C relationship in order to underpin our analysis of extensive carbon cycle model simulations.
The production, sinking and sequestration into the ocean interior of low-δ 13 C organic carbon via the soft-tissue component of the biological pump leads to a broad correlation between δ 13 C, [CO 2− 3 ] and macronutrients in the ocean (e.g.Hain et al., 2014a).The remineralization of this organic matter decreases δ 13 C and titrates [CO 2− 3 ], thereby reducing pH, while nutrient concentrations are increased.In waters that have experienced more soft tissue remineralization, both pH and δ 13 C will be lower (Fig. 5a, b).This is the dominant reason for the positive slope between δ 13 C and pH in the modern ocean (e.g.Foster et al., 2012;Fig. 5c).
Another significant factor affecting the spatial distribution of both δ 13 C and pH is seawater temperature, which affects both the equilibrium solubility of DIC and the equilibrium isotopic composition of DIC.Warmer ocean waters have decreased equilibrium solubility of DIC and so increased local [CO 2− 3 ] and pH (Goodwin and Lauderdale, 2013), while warmer waters have relatively low equilibrium δ 13 C values (Lynch-Stieglitz et al., 1995).This means that a spatial gradient in temperature acts to drive δ 13 C and pH in opposite directions: warmer waters tend to have higher pH but lower δ 13 C.These opposing temperature effects act to reduce the pH difference between two points with greatly different temperature to below the value expected based on δ 13 C alone.In other words, when using δ 13 C differences to estimate the pH gradient between the warm low-latitude surface waters and cold deep waters, the appropriate pHδ 13 C gradient will be less than expected when only considering the effect of organic carbon production, sinking and sequestration.For this reason, in our modelling analysis we focus on the warm surface to cold bottom pHδ 13 C rather than the slope of the overall pH-δ 13 C relationship, with the latter expected to be greater than the former.
In the modern ocean, and for the preceding tens of millions of years, the two dynamics described above are likely dominant in setting spatial variation in δ 13 C and pH (and [CO 2− 3 ]).However, other processes will have had a minor effect on either pH or δ 13 C.For instance, the dissolution of CaCO 3 shells increases [CO 2− 3 ] and pH (Broecker and Peng, 1982), but does not significantly affect δ 13 C (Zeebe and Wolf-Gladrow, 2001).Moreover, the long timescale of airsea isotopic equilibration of CO 2 combined with kinetic isotope fractionation during net carbon transfer is an important factor in setting the distribution of δ 13 C on a global ocean scale (Galbraith et al., 2015;Lynch-Stieglitz et al., 1995).Meanwhile, the effect of CO 2 disequilibrium on [CO 2− 3 ] and pH is modest (Goodwin and Lauderdale, 2013).

Modelling the pH-δ 13 C relationship
After correcting for the shift in δ 13 C due to anthropogenic activity, or Suess effect (Keeling 1979), modern global ocean observations demonstrate a near-linear relationship between global ocean data of in situ seawater pH and δ 13 C DIC , with a slope of 0.201 ± 0.005 (2σ ) (Foster et al., 2012;Fig. 5c).This empirically determined slope might well have been different in past oceans with very different nutrient cycling, carbon chemistry and circulation compared to today, and it does not appropriately represent the temperature effect de-  (Schlitzer 2016).pH data are from the CARINA data set (CARINA group, 2009) and the δ 13 C data are from the GLODAP data compilation (Key et al., 2004).(c) pH and δ 13 C DIC relationships in the modern ocean are adapted from Foster et al. (2012).Data are from all the ocean basins spanning approximately 40 • N to 40 • S. Because of anthropogenic acidification and the Suess effect, only data from > 1500 m are plotted.Also included in the plot are the data from a transect in the North Atlantic (from 0 to 5000 m) where the effects of anthropogenic perturbation on both parameters have been corrected (Olsen and Ninneman, 2010).scribed above (i.e.warm surface to cold bottom-water pHδ 13 C).Here we use an ensemble approach with two independent carbon cycle models to investigate changes in the pHδ 13 C regression.Below we provide pertinent information on the GENIE and CYCLOPS model experiments.
We use the Earth system model GENIE-1 (Edwards and Marsh, 2005;Ridgwell et al. 2007) to assess the robustness of the pHδ 13 C relationship and its sensitivity to physical and biogeochemical ocean forcing.The configuration used here is closely related to that of Holden et al. (2013), in which the controls on oceanic δ 13 C distribution were assessed, with an energy and moisture balance in the atmosphere, simple representations of land vegetation and sea ice, and frictional geostrophic ocean physics.In each of the 16 vertical levels in the ocean, increasing in thickness with depth, there are 36 × 36 grid cells (10 • in longitude and nominally 5 • in latitude, with higher resolution at low latitudes).Modern ocean bathymetry and land topography is applied in all simulations.The ocean biogeochemical scheme (Ridgwell et al., 2007) is based on conversion of DIC to organic carbon associated with phosphate uptake with fixed P : C : O stoichiometry.Organic carbon and nutrients are remineralized according to a remineralization profile with a predefined e-folding depth scale.This depth scale, as well as the rain ratio of inorganic to organic carbon in sinking particulate matter, is among the parameters examined in the sensitivity study.In these sim-ulations, there is no interaction with sediments.As a result of this, the steady state solutions reported here are reached within the 5000-year simulations, but they are not consistent with being in secular steady state with regard to the balance of continental weathering and ocean CaCO 3 burial.
The sensitivity study consists of seven sets of experiments, each varying a single model parameter relative to the control simulation with preindustrial atmospheric pCO 2 .This enables us to assess which processes, if any, are capable of altering the oceanic relationship between pH and δ 13 C and the uncertainty in the predictive skill of this relationship due to spatial variability.These experiments are therefore exploratory in nature and intended to study plausible range rather than determine magnitude of past changes.The seven parameters varied are (1) the ocean alkalinity reservoir; (2) the ocean's carbon reservoir; (3) the parameter "S.Lim gas exchange", which blocks air-sea gas exchange south of the stated latitude, significant here because of the dependence of δ 13 C on surface disequilibrium (Galbraith et al., 2015); (4) inorganic-to-organic carbon rain ratio, controlling the relationship between DIC and alkalinity distributions; ( 5 ment equivalent to freshwater hosing, leading to a shut down of the Atlantic meridional overturning circulation at low values; and (7) remineralization depth scale of sinking organic matter, which affects the vertical gradient of both pH and δ 13 C.A wide range of parameter values is chosen for each parameter in order to exceed any plausible changes within the Cenozoic.
For the second exploration of the controls on the slope of the pHδ 13 C relationship we use the CYCLOPS biogeochemical 18-box model that includes a dynamical lysocline, a subantarctic zone surface box and a polar Antarctic zone box (Sigman et al., 1998;Hain et al., 2010Hain et al., , 2014b)).The very large model ensemble with 13 500 individual model scenarios is designed to capture the full plausible range of (a) glacial-interglacial carbon cycle states by sampling the full solution space of Hain et al. (2010)  (5) range in seawater [Ca] concentration from 1× to 1.5× modern as per reconstructions (Horita et al., 2002); (6) Pacific CCD is set to the range of 4.4-4.9km via changes in the weathering flux, as per sedimentological evidence (Pälike et al., 2012); and (7) atmospheric CO 2 is set from 200 to 1000 ppm by changes in the "weatherability" parameter of the silicate weathering mechanism.
The ensemble spans predicted bulk ocean DIC between 1500 and 4500 µmol/kg, a wide range of ocean pH and CaCO 3 saturation states consistent with the open system weathering cycle, and widely different states of the oceanic biological pump.All 13 500 model scenarios are run for 2 million years after every single weatherability adjustment, part of the CCD inversion algorithm, guaranteeing the specified CCD depth and steady state with regard to the balance of continental weathering and ocean CaCO 3 burial for the final solution (unlike the GENIE simulations, CaCO 3 burial was entirely neglected due to computational cost of the long model integrations it would require).The inverse algorithm typically takes at least 10 steps to conversion, resulting in ∼ 300 billion simulated years for this ensemble.This range of modelling parameters was chosen to exceed the range of carbonate system and ocean circulation changes that can be expected for the Neogene based on records of [Ca] and [Mg] (Horita et al., 2002), CCD changes (Pälike et al., 2012), atmospheric CO 2 (Beerling and Royer, 2011) and records of glacial-interglacial circulation change (Curry and Oppo, 2005).

Assessing uncertainty
δ 11 B sw uncertainty was calculated using a Monte Carlo approach where pH was calculated for deep and surface waters at each time slice using a random sampling (n = 100 00) of the various input parameters within their respective uncertainties as represented by normal distributions.These uncertainties (2σ uncertainty in parentheses) are temperature (±2 • C), salinity (±2 units on the practical salinity scale), [Ca] (±4.5 mmol kg −1 ), [Mg] (±4.5 mmol kg −1 ), δ 11 B planktic (±0.15-0.42‰) and δ 11 B benthic (±0.21-0.61‰).For the estimate of the surfaceto-sea floor pH gradient, we use the central value of the pHδ 13 C relationship diagnosed from our CYCLOPS and GENIE sensitivity experiments (i.e.0.175 ‰ −1 ; see Sect.3.2 below) and then we assign a ±0.05 uncertainty range with a uniform probability (rather than a normal distribution) to the resulting surface-to-sea floor pH estimate (see also Table 2).Thus, the magnitude of this nominal uncertainty is equivalent to a 0.14-0.21‰ −1 pH / δ 13 C uncertainty range that spans the vast majority of our CYCLOPS and GENIE simulations and the prediction error (RMSE) of fitting a linear relationship to the GENIE pH and δ 13 C output (see Sect. (4) where [ 11 B] is the intensity of 11 B signal in volts and Eqs. ( 4) and ( 5), used with 10 11 and 10 12 resistors respectively.From the 10 000 Monte Carlo ensemble solutions of our 22 benthic-planktic pairs, we construct 10 000 randomized records of δ 11 B sw as a function of time.Each of these randomized δ 11 B sw records are subjected to smoothing using the locally weighted scatterplot smoothing (LOWESS) algorithm with a smoothing parameter (span) of 0.7.The purpose of the smoothing is to put some controls on the rate at which the resulting individual Monte Carlo δ 11 B sw records are allowed to change, which in reality is limited by the seawater boron mass balance (∼ 0.1 ‰ per million years; boron residence time is 11-17 million years; Lemarchand et al., 2000).Our choice of smoothing parameter allows for some of the individual Monte Carlo records to change as fast as ∼ 1 ‰ per million years, although in reality the average rate of change is much smaller than this (see Sect. 3.3).Consequently, this method removes a significant amount of uncorrelated stochastic noise (resulting from the uncertainty in our Set via silicate weatherability 200, 300, 400, 500, 600, 700, 800, 900, 1000 ppm * The six parameters assume 5, 5, 5, 2, 9 and 6 values, yielding 13 500 distinct parameter combinations.* * These parameters are intended to span the full range of ocean carbon cycling over late Pleistocene glacial-interglacial cycles, as described in more detail in Hain et al. (2010).* * * We selected representative time slices based on local extrema in the CCD reconstruction of Pälike et al. (2012) and we combine these with appropriate reconstructed calcium concentrations based on Horita et al. (2002).These choices are intended to capture the range of long-term steady state conditions of the open system CaCO 3 cycle relevant to our study interval.* * * * These atmospheric CO 2 levels are chosen to span a range wider than expected for the study interval.Following the silicate-weathering-feedback paradigm, long-term CO 2 is fully determined by the balance of geologic CO 2 sources and silicate weathering, whereby faster acting processes of the open system CaCO 3 cycle compensate relative to that CO 2 level.All else equal, high CO 2 levels, low calcium concentrations and deep CCD correspond to high bulk ocean carbon concentrations (Hain et al., 2015), with many of the individual simulations of this ensemble exceeding 4000 µM DIC.±4.5 mmol kg −1 Following Horita et al. (2002) input parameters), while not smoothing away the underlying signal.As a result of anomalously low δ 11 B differences (< 1 ‰) between benthic and planktic pairs, two pairs at 8.68 and 19 Ma were discarded from the smoothing.It may be possible that preservation is not so good within these intervals and the planktic foraminifera are affected by partial dissolution (Seki et al., 2010).The spread of the ensemble of smoothed δ 11 B sw curves represents the combination of the compounded, propagated uncertainties of the various inputs (i.e.Monte Carlo sampling) with the additional constraint of gradual δ 11 B sw change over geological time imposed by the inputs and outputs of boron to the ocean and the total boron inventory (i.e. the smoothing of individual Monte Carlo members).Various statistical properties (i.e.mean, median, standard deviation (σ ), various quantiles) of this δ 11 B sw reconstruction were evaluated from the ensemble of smoothed δ 11 B sw records.Generally, for any given benthic-planktic pair the resulting δ 11 B sw estimates are not perfectly normally distributed and thus we use the median as the metric for the central tendency (i.e.placement of marker in Fig. 10).

δ 11 B benthic and planktic data
Surface and deep-ocean δ 11 B broadly show a similar but inverse pattern to δ 13 C and temperature throughout the Neogene (Fig. 6).The δ 11 B benthic record decreases from ∼ 15 ‰ at 24 Ma to a minimum of 13.28 ‰ at 14 Ma before increasing to ∼ 17 ‰ at present day (Fig. 6).This pattern and the range of values in benthic foraminiferal δ 11 B are in keeping with previously published Neogene δ 11 B benthic records measured using NTIMS (Raitzsch and Hönisch, 2013).This suggests that our deep-water δ 11 B record is representative of large-scale pH changes in the global ocean.
While the surface δ 11 B planktic remained relatively constant between 24 and 11 Ma at ∼ 16 ‰, there has been a significant increase in δ 11 B between the middle Miocene and present (values increase to ∼ 20 ‰) (Fig. 6b).The reconstructed surface water temperatures show a long-term decrease through the Neogene from ∼ 28 to 24 • C, aside from during the Miocene Climatic Optimum (MCO) where maximum Neogene temperatures were reached (Fig. 6c).Following Cramer et al. (2011), deep-water temperatures decrease from ∼ 12 to 4 • C at present day and similarly show maximum temperatures in the MCO.Surface and deep-water δ 13 C DIC both broadly decrease through the Neogene and appear to co-vary on shorter timescales (Fig. 6e, f).

The relationship between δ 13 C and pH gradients
In the global modern ocean data, after accounting for the anthropogenic carbon, the empirical relationship between in situ pH and DIC δ 13 C is well described by a linear function with a slope of 0.201 ± 0.005 (2σ ) (Fig. 5; Foster et al., 2012).However, this slope is only defined by surface waters in the North Atlantic due to a current lack of modern data where the impact of the Suess effect has been corrected (Olsen and Ninneman, 2010).Consequently, we are not currently able to determine the slope between the warm surface and cold deep ocean in the modern ocean at our sites.Instead, here we use the two modelling experiments to define this slope.In the control GENIE experiment (green star, Fig. 7), the central value for the slope of the pH-δ 13 C relationship is slightly greater than 0.2 ‰ −1 for the full three-dimensional data regression (not shown) and about 0.175 ‰ −1 for the warm surface-cold deep pHδ 13 C relationship (Fig. 7); this is consistent with the theory for the effect of temperature gradients (see Sect. 2.3).For both ways of analysing the GENIE output, the prediction uncertainty of the regressions, the RMSE, is ∼ 0.05 ‰ −1 under most conditions (open red circles in Fig. 7), with the exception of cases where large changes in either DIC or alkalinity (ALK) yield somewhat larger changes in the relationship between pH and δ 13 C (see below).
In our CYCLOPS model ensemble, the central values of the slopes of the full three-dimensional pH-δ 13 C regressions and of the warm surface-to-cold deep pHδ 13 C are 0.2047 ‰ −1 (1σ of 0.0196 ‰ −1 ; Fig. 8a) and 0.1797 ‰ −1 (1σ of 0.0213 ‰ −1 ; Fig. 8b) respectively.If we restrict our analysis of the CYCLOPS ensemble to only the Atlantic basin warm surface-to-cold deep pHδ 13 C, where most of our samples come from, we find a relationship of only 0.1655 ‰ −1 (1σ of 0.0192 ‰ −1 ; Fig. 8c).That is, overall, we find near-perfect agreement between modern empirical data and our GENIE and CYCLOPS experiments.Encouraged by this agreement we select the warm surface-to-cold deep pHδ 13 C central value of 0.175 ‰ −1 to estimate the surface-sea floor pH difference from the planktic-benthic foraminifera δ 13 C difference.To account for our ignorance as to the accurate value of pHδ 13 C in the modern ocean, its temporal changes over the course of the study interval and the inherent prediction error from using a linear pHδ 13 C relationship, we assign a nominal uniform uncertainty range of ±0.05 around the central pH estimate for the purpose of Monte Carlo uncertainty propagation.Our analysis also suggests that where surface-to-thermocline plankticplanktic gradients are employed, the plausible pHδ 13 C range should be significantly higher than applied here to account for the relatively lower temperature difference.Based on the appropriate pHδ 13 C relationship, we reconstruct a time varying surface-to-deep pH gradient, which ranges between 0.14 and 0.35 pH units over our study interval (Fig. 9); we also apply a flat uncertainty of ±0.05.The reconstructed pH gradient remains broadly within the range of the modern values (0.19 to 0.3), although there is some evidence of multimillion-year scale variability (Fig. 9).
As a caveat to our usage of the pHδ 13 C relationship, we point to changes in that relationship that arise in our GE-NIE sensitivity experiments where carbon and alkalinity inventories are manipulated, which can yield values outside of what is plausible.We note that our CYCLOPS ensemble samples a very much wider range of carbon and alkalinity inventories, with pHδ 13 C remaining inside that range.While CYCLOPS simulates the balance between weathering and CaCO 3 burial, which is known to neutralize sudden carbon or alkalinity perturbations on timescales much less than 1 million years, the configuration used for our GENIE simulations does not and is therefore subject to states of ocean carbon chemistry that can safely be ruled out for our study inter- val and likely for most of the Phanerozoic.The differing outputs from CYCLOPS and GENIE in the DIC and ALK experiments shows that pHδ 13 C depends on background seawater acid-base chemistry in ways that are not yet fully understood.That said, the generally coherent nature of our results confirms that we likely constrain the plausible range of pHδ 13 C for at least the Neogene, if not the entire Cenozoic, outside of extreme events such as the Palaeocene-Eocene Thermal Maximum.

δ 11 B sw record through the Neogene
Using input parameter uncertainties as described in section 2.5 yields individual Monte Carlo member δ 11 B sw estimates between 30 and 43.5 ‰ at the overall extreme points and typically ranging by ∼ 10 ‰ (dashed in Fig. 10a) for each time point.This suggests that the uncertainties we assign to the various input parameters are generous enough not to predetermine the quantitative outcomes.Fig. 10a).The δ 11 B sw for Plio-Pleistocene time points cluster around ∼ 40 ‰, while middle-late Miocene values cluster around ∼ 36.5 ‰.The estimates at individual time points are completely independent from each other, such that the observed clustering is strong evidence for an underlying longterm signal in our data, albeit one that is obscured by the uncertainties involved in our individual δ 11 B sw estimates.The same long-term signal is also evident when pooling the individual Monte Carlo member δ 11 B sw estimates into 8-millionyear bins and evaluating the mean and spread (2σ ) in each bin (Fig. 10b).This simple treatment highlights that there is a significant difference between our Plio-Pleistocene and middle Miocene data bins at the 95 % confidence level and that δ 11 B sw appears to also have been significantly lower than present day during the early Miocene.

Data smoothing
The ∼ 1 to 4 ‰ likely ranges for δ 11 B sw would seem to be rather disappointing given the goal to constrain δ 11 B sw for pH reconstructions.However, most of that uncertainty is stochastic random error that is uncorrelated from time point to time point.Furthermore, we know from mass balance considerations that δ 11 B sw of seawater should not change by more than ∼ 0.1 ‰ per million years (Lemarchand et  Surface-sea floor
Figure 9.The pH gradient between surface and deep water through time calculated from the δ 13 C gradient and using a flat probability derived from the low-latitude ensemble regressions from the CY-CLOPS model.The modern pH gradients at each site are also plotted.
2000) because of the size of the oceanic boron reservoir compared the inputs and outputs (see Fig. 1).We use this as an additional constraint via the LOWESS smoothing we apply to each Monte Carlo time series.One consideration is that every individual Monte Carlo δ 11 B sw estimate is equally likely and the smoothing should therefore target randomly selected individual Monte Carlo δ 11 B sw estimates, as we do here, rather than smoothing over the likely ranges identified for each time point.In this way the smoothing becomes an integral part of our Monte Carlo uncertainty propagation and the spread among the 10 000 individual smoothed δ 11 B sw curves carries the full representation of propagated input uncertainty conditional on the boron cycle mass balance constraint.A second consideration is that the smoothing should only remove noise, not underlying signal.As detailed above, for this reason the smoothing parameter we chose has enough freedom to allow the δ 11 B sw change to be dictated by the data, with only the most extreme shifts in δ 11 B sw removed.We also tested the robustness of the smoothing procedure itself (not shown) and found only marginal changes when changing algorithm (LOESS versus LOWESS, with and without robust option) or when reducing the amount of smoothing (i.e.increasing the allowed rate δ 11 B sw change).The robustness of our smoothing is further underscored by the good correspondence with the results of simple data binning (Fig. 10b).

Comparison to other δ 11 B sw records
The comparison of our new δ 11 B sw record to those previously published reveals that despite the differences in methodology, the general trends in the records show excel- Figure 10.The calculated δ 11 B sw from the benthic-planktic δ 11 B pairs using a pH gradient derived from δ 13 C .The uncertainty on each data point is determined using a Monte Carlo approach, including uncertainties in temperature, salinity, δ 11 B and pH gradient (see text for details).Data are plotted as box and whisker diagrams where the median and interquartile ranges as plotted in the box and whiskers show the maximum and minimum output from the Monte Carlo simulations.The line of best fit is the probability maximum of a LOWESS fit given the uncertainty in the calculated δ 11 B sw .The darker shaded area highlights the 68 % confidence interval and the lighter interval highlights the 95 % confidence interval.The bottom panel shows box plots of the mean and 2 standard error (SE) of "binning" the individual δ 11 B sw measurements into 8 My intervals.The middle line is the mean and the box shows the 2 SE of the data points in that bin.The smoothed record is also plotted for comparison where the line of best fit is the probability maximum of a LOWESS fit given the uncertainty in the calculated δ 11 B sw .The darker shaded area highlights the 68 % confidence interval and the lighter interval highlights the 95 % confidence interval.The black dot is the modern value of 39.61 ‰ (Foster et al., 2010).lent agreement.The most dominant common feature of all the existing estimates of Neogene δ 11 B sw evolution is an increase through time from the middle Miocene to the Plio-Pleistocene (Fig. 11).While the model-based δ 11 B sw record of Lemarchand et al. (2000) is defined by a monotonous and very steady rise over the entire study interval, all three measurement-based records, including our own, are characterized by a single dominant phase of increase between roughly 12 and 5 Ma.Strikingly, the Pearson and Palmer (2000) record falls almost entirely within our 95 % likelihood envelope, overall displaying very similar patterns of long-term change but with a relatively muted ampli-This study Pearson and Palmer, 2000Lemarchand et al., 2002Foster et al., 2012Raitzsch and Honisch, 2013(Klochko) Raitzsch and Honisch, 2013 (Empirical)
Figure 11.The δ 11 B sw curve calculated using the variable pH gradient derived from δ 13 C.The median (red line), 68 % (dark red band) and 95 % (light red band) confidence intervals are plotted.Plotted with a compilation of published δ 11 B sw records.Seawater composition reconstructed from foraminifera depth profiles (light blue squares and dark blue crosses) from Pearson and Palmer (2000) and Foster et al. (2012) respectively; numerical modelling (green line), with additional green lines showing ±1 ‰ confidence interval (Lemarchand et al., 2000); and benthic δ 11 B (purple diamonds and dark purple line showing 5-point moving average is using the fractionation factor of Klochko et al., 2006, light purple line showing 5-point moving average using an empirical calibration) from Raitzsch and Hönisch (2013).All the published δ 11 B sw curves are adjusted so that at t = 0 and the isotopic composition is equal to the modern (39.61 ‰).
tude and overall rate of change relative to our reconstruction.Conversely, some of the second-order variations in the reconstruction by Raitzsch and Hönisch (2013) are not well matched by our reconstruction.However, the dominant episode of rapid δ 11 B sw rise following the middle Miocene is in almost perfect agreement.We are encouraged by these agreements resulting from approaches based on very different underlying assumptions and techniques, which we take as indication for an emerging consensus view of δ 11 B sw evolution over the last 25 My and as a pathway towards reconstructing δ 11 B sw further back in time.
The record by Pearson and Palmer (2000) is well correlated to our reconstruction, but especially during the early Miocene there is a notable ∼ 0.5 ‰ offset (Fig. 11).This discrepancy could be due to a number of factors.Firstly, the applicability of this δ 11 B sw record (derived from δ 11 B data measured using NTIMS) to δ 11 B records generated using the MC-ICPMS is uncertain (Foster et al., 2013).In addition, this δ 11 B sw record is determined using a fractionation factor of 1.0194 (Kakihana et al., 1977), whereas recent experimental data have shown the value to be higher (1.0272 ± 0.0006, Klochko et al., 2006), although foraminiferal vital effects are likely to mute this discrepancy.Thirdly, given our understanding of the δ 11 B difference between species-size fractions (Foster, 2008;Henehan et al., 2013), the mixed species and size fractions used to make the δ 11 B measurements in that study may have introduced some additional uncertainty in the reconstructed δ 11 B sw .Conversely, there is a substantial spread between our three time points during the earliest Miocene, which combined with the edge effect of the smoothing, gives rise to a widening uncertainty envelope during the time of greatest disagreement with Pearson and Palmer (2000).This could be taken as an indication that our reconstruction, rather than that of Pearson and Palmer, is biased during the early Miocene.
The δ 11 B sw record calculated using benthic δ 11 B and assumed deep-ocean pH changes (Raitzsch and Hönisch, 2013) is also rather similar to our δ 11 B sw reconstruction.The discrepancy between the two records in the early Miocene could plausibly be explained by bias in our record (see above) or may in part be a result of the treatment of surface water pH in the study of Raitzsch and Hönisch (2013) and their assumption of a constant surface-deep pH gradient (see Fig. 9).The combined output from two carbon cycle box models is used to make the assumption that surface ocean pH near-linearly increased by 0.39 over the last 50 My.The first source of surface water pH estimates is from the study of Ridgwell et al. (2005), where CO 2 proxy data (including some derived using the boron isotope-pH proxy) is used, leading to some circularity in the methodology.The second source of surface water pH estimates is from Tyrrell and Zeebe (2004)   lem does not apply.While this linear pH increase broadly matches the CO 2 decline from proxy records between the middle Miocene and present, it is at odds with the CO 2 proxy data during the early Miocene that show that CO 2 was lower than the middle Miocene during this interval (Beerling and Royer, 2011).Consequently, the proxy CO 2 and surface water pH estimates may not be well described by the linear change in pH applied by Raitzsch and Hönisch (2013) across this interval, potentially contributing to the discrepancy between our respective δ 11 B sw reconstructions.
Our new δ 11 B sw record falls within the broad uncertainty envelope of boron mass balance calculations of Lemarchand et al. (2000), but those modelled values do not show the same level of multimillion-year variability of either Raitzsch and Hönisch (2013) or our new record.This therefore suggests that the model does not fully account for aspects of the changes in the ocean inputs and outputs of boron through time on timescales less than ∼ 10 million years.
In line with the conclusions of previous studies (e.g.Raitzsch and Hönisch, 2013), our data show that the δ 11 B sw signal in the fluid inclusions (Paris et al., 2010) is most likely a combination of the δ 11 B sw and some other factor such as a poorly constrained fractionation factor between the seawater and the halite.Brine-halite fractionation offsets of −20 to −30 ‰ and −5 ‰ are reported from laboratory and natural environments (Vengosh et al., 1992;Liu et al., 2000).These fractionations and riverine input during basin isolation will drive the evaporite-hosted boron to low-δ 11 B isotope values such that the fluid inclusion record likely provides a lower limit for the δ 11 B sw through time (i.e.δ 11 B sw is heavier than the halite fluid inclusions of Paris et al., 2010).For this halite record to be interpreted directly as δ 11 B sw , a better understanding of the factor(s) controlling the fractionation during halite formation and any appropriate correction need to be better constrained.

Common controls on the seawater isotopic ratios of B, Mg, Ca and Li
Our new record of δ 11 B sw has some substantial similarities to secular change seen in other marine stable isotope records (Fig. 12).The lithium isotopic composition of seawater (δ 7 Li sw ; Misra and Froelich, 2012) and the calcium isotopic composition of seawater as recorded in marine barites (δ 44/40 Ca sw ; Griffith et al., 2008) both increase through the Neogene, whereas the magnesium isotopic composition of seawater (δ 26 Mg sw ) decreases (Pogge von Strandmann et al., 2014), suggesting a similar control on the isotopic composition of all four elements across this time interval (Fig. 12).To further evaluate the correlation between these other marine isotope records and δ 11 B sw , we interpolate and cross plot δ 11 B sw and the δ 7 Li sw , δ 44/40 Ca sw and δ 26 Mg sw records.This analysis suggests that the isotopic composition of δ 11 B sw , δ 7 Li sw , δ 26 Mg sw and δ 44/40 Ca sw are well correlated through the Neogene, although there is some scatter in these relationships (Fig. 13).Although the Sr isotope record shows a similar increase during the Neogene (Hodell et al., 1991), we focus our discussion on δ 11 B sw , δ 7 Li sw , δ 26 Mg sw and δ 44/40 Ca sw , given that the factors fractionating these stable isotopic systems are similar (see below).
To better constrain the controls on δ 11 B sw , δ 7 Li sw , δ 26 Mg sw and δ 44/40 Ca sw , it is instructive to compare the size and isotopic composition of the fluxes of boron, lithium, calcium and magnesium to the ocean (Table 3).The major flux of boron into the ocean is via riverine input (Lemarchand et al., 2000), although some studies suggest that atmospheric input may also play an important role (Park and Schlesinger, 2002).The loss terms are dominated by adsorption onto clays and the alteration of oceanic crust (Spivack and Edmond, 1987;Smith et al., 1995).Similarly, the primary inputs of lithium into the ocean come from hydrothermal sources and riverine input and the main outputs are ocean crust alteration and adsorption onto sediments (Misra and Froelich, 2012).The three dominant controls on magnesium concentration and isotope ratio in the oceans are the riverine input, ocean crust alteration and dolomitization (Table 3) (Tipper et al., 2006b).The main controls on the amount of calcium in the modern ocean and its isotopic composition are the balance between riverine and hydrothermal inputs and removal through CaCO 3 deposition and alteration of oceanic crust (Fantle and Tipper, 2014;Griffith et al., 2008).Dolomitization has also been cited as playing a potential role in controlling δ 44/40 Ca sw , although the contribution of this process through time is poorly constrained (Griffith et al., 2008).Analysis of the oceanic fluxes of all four ions suggests that riverine input may be an important factor influencing the changing isotopic composition of B, Li, Ca and Mg over the late Neogene (Table 3).In the case of all four elements, a combination of the isotopic ratio of the source rock and isotopic fractionation during weathering processes are typically invoked to explain the isotopic composition of a particular river system.However, in most cases the isotopic composition of the source rock is found to be of secondary importance (Rose et al., 2000;Kisakũrek et al., 2005;Tipper et al., 2006b;Millot et al., 2010).For instance, the δ 11 B composition of rivers is primarily dependent on isotopic fractionation during the reaction of water with silicate rocks and to a lesser extent the isotopic composition of the source rock (i.e. the proportion of evaporites and silicate rocks; Rose et al., 2000).While some studies have suggested that the isotope composition of rainfall within the catchment area may be an important factor controlling the δ 11 B in rivers (Rose-Koga et al., 2006), other studies have shown atmospheric boron to be a secondary control on riverine boron isotope composition (Lemarchand and Gaillardet, 2006).The source rock also appears to have limited influence on the δ 7 Li composition of rivers, and riverine δ 7 Li varies primarily with weathering intensity (Kisakũrek et al., 2005;Millot et al., 2010).The riverine input of calcium to the oceans is controlled by the composition of the primary continental crust (dominated by carbonate weathering) and a recycled component, although the relative influence of these two processes is not well understood (Tipper et al., 2006a).In addition, vegetation may also play a significant role in the δ 44/40 Ca of rivers (Fantle and Tipper, 2014).For Mg, the isotopic composition of the source rock is important for small rivers; however, lithology is of limited significance at a global scale in comparison to fractionation in the weathering environment (Tipper et al., 2006b).Given the lack of evidence of source rock as a dominant control on the isotopic composition of rivers, here we focus on some of the possible causes for changes in the isotopic composition and/or flux of riverine input over the Neogene.
In this regard, of the four elements discussed here, the Li isotopic system is the most extensively studied.Indeed, the change in δ 7 Li sw has already been attributed to an increase in the δ 7 Li sw composition of the riverine input (Hathorne and James, 2006;Misra and Froelich, 2012).The causes of the shift in δ 7 Li riverine have been variably attributed to (1) an increase in incongruent weathering of silicate rocks and sec-ondary clay formation as a consequence of Himalayan uplift (Misra and Froelich, 2012;Li and West, 2014), (2) a reduction in weathering intensity (Hathorne and James, 2006;Froelich and Misra, 2014;Wanner et al., 2014), (3) an increase in silicate weathering rate (Liu et al., 2015), (4) an increase in the formation of floodplains and the increased formation of secondary minerals (Pogge von Strandmann and Henderson, 2014), and (5) a climatic control on soil production rates (Vigier and Godderis, 2015).In all five cases the lighter isotope of Li is retained on land in clay and secondary minerals.A mechanism associated with either an increase in secondary mineral formation or the retention of these minerals on land is also consistent across Mg, Ca and B isotope systems.For instance, clay minerals are preferentially enriched in the light isotope of B (Spivack and Edmond, 1987;Deyhle and Kopf, 2004;Lemarchand and Gaillardet, 2006) and Li (Pistiner and Henderson, 2003), and soil carbonates and clays are preferentially enriched in the light isotope of Ca (Tipper et al., 2006a;Hindshaw et al., 2013;Ockert et al., 2013).The formation of secondary silicate minerals, such as clays, is assumed to preferentially take up the heavy Mg isotope into the solid phase (Tipper et al., 2006a, b;Pogge von Strandmann et al., 2008;Wimpenny et al., 2014), adequately explaining the inverse relationship between δ 11 B sw and δ 26 Mg sw .Consequently, the increased formation or retention on land of secondary minerals would alter the isotopic composition of the riverine input to the ocean in the correct direction to explain the trends in all four isotope systems through the late Neogene (Fig. 13).While the relationships between the different isotope systems discussed here suggest a common control, the influence of carbonate and dolomite formation on Ca and Mg isotopes are also likely to have played a significant role in the evolution of these isotope systems (Tipper et al., 2006b;Fantle and Tipper, 2014).Consequently, a future model of seawater chemistry evolution through the Neogene must also include these additional factors.Further exploration is also needed to determine the influence of residence time on the evolution of ocean chemistry.Nonetheless, given the similarities between the geochemical cycles of B and Li, and despite the large difference in residence time (Li = 1 million years, B = 11-17 million years), the correlation between these two records is compelling and would no doubt benefit from additional study.

Conclusions
Here we present a new δ 11 B sw record for the Neogene (0-23 Ma) based on paired planktic-benthic δ 11 B measurements.Our new record suggests that δ 11 B sw (i) was ∼ 37.5 ‰ at the Oligocene-Miocene boundary (23 Ma), (ii) remained low through the middle Miocene (16-12 Ma), (iii) rapidly increased to the modern value during the late Miocene (between 12 and 5 Ma), and (iv) plateaued at modern values over the Plio-Pleistocene (5 Ma to present).Despite some dis-agreements, and different uncertainties associated with each approach, the fact that our new record and both of the published data-based reconstructions capture the first-order late Miocene δ 11 B sw rise suggests that consensus is building for the δ 11 B sw evolution through the Neogene.This emerging view on δ 11 B sw change provides a vital constraint required to quantitatively reconstruct Neogene ocean pH, ocean carbon chemistry and atmospheric CO 2 using the δ 11 B-pH proxy.When our new δ 11 B sw record is compared to changes in the seawater isotopic composition of Li, Ca and Mg, the shape of the records across the Neogene is remarkably similar.For all four systems, riverine input is cited as a common and key control of the isotopic composition of the respective elements in seawater.When we compare the isotopic fractionation of the elements associated with secondary mineral formation, the trends in the δ 26 Mg sw , δ 44/40 Ca sw , δ 11 B sw and δ 7 Li sw records are all consistent with an increase in secondary mineral formation through time.While a more quantitative treatment of these multiple stable isotope systems is required, the δ 11 B sw record presented here provides additional constraints on the processes responsible for the evolution of ocean chemistry through time.

Figure 3 .
Figure 3. Schematic diagram showing the change in pH gradient with a 3 ‰ change in δ 11 B for δ 11 B sw of (a) 39.6 and (b) 37.5 ‰.Arrows highlight the different pH gradients.Note how a δ 11 B difference of 3 ‰ is translated into different pH gradients depending on the δ 11 B sw .Calculated using B T = 432.6 µmol kg −1(Lee et al., 2010) and α B = 1.0272(Klochko et al., 2006).(c) The pH change for a δ 11 B change of 3 ‰ at a range of different δ 11 B sw .

Figure 4 .
Figure 4. Map of study sites and mean annual air-sea disequilibrium with respect to pCO 2 .The black dots indicate the location of the sites used in this study.ODP sites 758, 999, 926 and 761 used in this study are highlighted with water depth.Data are from Takahashi et al., 2009, plotted using Ocean Data View(Schlitzer, 2016).

Figure 5 .
Figure 5. Latitudinal cross section through the Atlantic showing (a) pH variations and (b) the δ 13 C composition.Data are plotted using Ocean Data View(Schlitzer 2016).pH data are from the CARINA data set(CARINA group, 2009)  and the δ 13 C data are from the GLODAP data compilation(Key et al., 2004).(c) pH and δ 13 C DIC relationships in the modern ocean are adapted fromFoster et al. (2012).Data are from all the ocean basins spanning approximately 40 • N to 40 • S. Because of anthropogenic acidification and the Suess effect, only data from > 1500 m are plotted.Also included in the plot are the data from a transect in the North Atlantic (from 0 to 5000 m) where the effects of anthropogenic perturbation on both parameters have been corrected(Olsen and Ninneman, 2010).
) "Antarctic shelf FWF", a freshwater flux adjustment (always switched off in control experiments with GE-NIE) facilitating the formation of brine-rich waters, which produces a high-salinity poorly ventilated deep ocean at high values; (6) "Atlantic-Pacific FWF", a freshwater flux adjust-R.Greenop et al.: A record of Neogene seawater δ 11 B 3.2 below).The uncertainty in the δ 11 B measurements is calculated from the long-term reproducibility of Japanese Geological Survey Porites coral standard (JCP; δ 11 B = 24.3‰) at the University of Southampton using the equations 2σ = 2.25exp −23.01[ 11 B] + 0.28 exp −0.64[ 11 B]

Table 3 .
The average δ 11 B, δ 26 Mg, δ 44/40 Ca and δ 7 Li composition of major fluxes into and out of the ocean.Colour coding reflects the relative importance of each of the processes (darker shading reflects greater importance).The colour coding for boron is based onLemarchand et al. (2000) and references therein, lithium onMisra and Froelich (2012) and references therein, magnesium on Tipper et al. (2006b) and calcium on Fantle and Tipper (2014) and Griffin et al. (2008) and references therein.Modern δ 26 Mg sw and δ 11 B sw are from Foster et al. (2010) and δ 7 Li sw is from Tomascak (2004).The δ 44/40 Ca presented here was measured relative to seawater and hence seawater has a δ 44/40 Ca sw of zero per mill by definition.The isotopic ratio of each process is (a) Lemarchand et al. (2000) and references therein; (b) Misra and Froelich (2012) and references therein; (c) Burton and Vigier (2012); (d) Tipper et al. (2006b); (e) Wombacher et al. (2011); (f) includes dolomitization; (g) removal through hydrothermal activity; (h) Griffith et al. (2008); (i) Fantle and Tipper (2014) and references therein; (j) dolomitization may be an important component of the carbonate flux.n/a = not applicable.

Figure 6 .
Figure 6.δ 11 B planktic , temperature and δ 13 C DIC estimates for the surface and deep ocean through the last 23 million years.(a) δ 11 B planktic surface; (b) δ 11 B borate deep from benthic foraminifera (blue) from this study and (green) Raitzsch and Hönisch (2013).The error bars show the analytical external reproducibility at 95 % confidence for this study.For the Raitzch and Hönisch (2013) data the error bars represent propagated uncertainties of external reproducibilities of time equivalent benthic foraminifer samples from different core sites in different ocean basins; (c) Mg / Ca based temperature reconstructions of surface dwelling planktic foraminifera; (d) deep-water temperature estimates from Cramer et al. (2011); (e) δ 13 C DIC surface record; and (f) δ 13 C DIC benthic record.Squares depict ODP Site 999, triangles are ODP Site 758, diamonds are ODP Site 926 and circles are ODP Site 761.Species are highlighted by colour: Orange is T. trilobus, purple G. ruber, pink G. praebulloides, dark blue Cibicidoides wuellerstorfi and light blue Cibicidoides mundulus.The two benthic-planktic pairs that were removed prior to smoothing are highlighted with arrows.

Figure 7 .
Figure7.The output from GENIE sensitivity analysis showing the warm surface-to-cold deep pHδ 13 C relationship.A pre-industrial model set-up was taken and perturbations were made to alkalinity inventory, carbon inventory, Antarctic shelf freshwater flux (Sv), Atlantic-Pacific freshwater flux, S. Lim gas exchange (blocks air-sea gas exchange south of the stated latitude), remineralization depth scale (m) and rain ratio -as described in the methods section.Blue circles depict the pHδ 13 C relationship (where the colours reflect the CO 2 level of each experiment) and red open circles show the RMSE.The green stars are the pHδ 13 C relationship for the control experiment conducted at 292.67 ppm CO 2 .The green (open) points show the RMSE for this control run.Inventories are dimensionless (1 is control).For the Atlantic-Pacific FWF 1 is equivalent to 0.32 Sv.The alkalinity and carbon inventory experiments are very extreme and inconsistent with geologic evidence.

Figure 8 .
Figure 8.The output from sensitivity analysis of the relationship between pH gradient and δ 13 C gradient from the 13 500-run CYCLOPS ensemble (see text for model details).Panel (a) shows the mean gradient when the results from all 18 ocean boxes are included in the regression.Panel (b) shows only the boxes from the low-latitude ocean from all basins and (c) shows the regression from only North Atlantic low-latitude boxes.Note the lower pH / δ 11 B slope at the lower latitudes due to the effect of temperature.The 0.201 line in each panel is the mean gradient when all the ocean boxes are included in the regression.

Figure 12 .
Figure 12.(a) The δ 11 B sw curve from this study plotted with other trace element isotopic records.On the δ 11 B sw panel the darker shaded area highlights the 68 % confidence interval and the lighter interval highlights the 95 % confidence interval, (b) δ 26 Mg sw record from Pogge von Strandmann et al. (2014) (error bars are ±0.28 ‰ and include analytical uncertainty and scatter due to the spread in modern O. universa and the offset between the two analysed species), (c) δ 44/40 Ca sw record from Griffith et al. (2008) (error bars show 2σ uncertainty) and (d) δ 7 Li sw record from Misra and Froelich (2012) (error bars show 2σ uncertainty).Blue dashed lines show middle Miocene values and red dashed lines highlight the modern.

Figure 13 .
Figure 13.Crossplots of the records of δ 11 B sw using a variable pH gradient derived from δ 13 C (error bars show 2σ uncertainty), with δ 44/40 Ca sw from Griffith et al. (2008) (error bars show 2σ uncertainty), δ 7 Li sw from Misra and Froelich (2012) (error bars show 2σ uncertainty) and δ 26 Mg sw from Pogge von Strandmann et al. (2014) (error bars are ±0.28 ‰ and include analytical uncertainty and scatter due to the spread in modern O. universa and the offset between the two analysed species).The colour of the data points highlights the age of the data points where red = modern and blue = 23 Ma.
(Elderfield et al., 2006)o Fisher Scientific Element 2 XR.Al / Ca was also measured to assess the competency of the sample cleaning.Because of complications with the Mg / Ca temperature proxy in Cibicidoides species(Elderfield et al., 2006), bottom-water temperatures (BWTs) are estimated here by taking the global secular temperature change from the Mg / Ca temperature compilation ofCramer et  al. ( and (b) reconstructed secular changes in seawater[Ca](calcium concentration), carbonate compensation depth (CCD), weathering and atmospheric CO 2 (Table1).The following seven model parameters are systematically sampled to set the 13 500 model scenarios: (1) shallow-versus deep-Atlantic meridional overturning circulation represented by modern reference North Atlantic deep water (NADW) versus peak glacial North Atlantic intermediate water (GNAIW) circulation; (2) irondriven changes in nutrient drawdown in the subantarctic zone of the Southern Ocean; (3) changes in nutrient drawdown of the polar Antarctic; (4) changes in vertical exchange between the deep Southern Ocean and the polar Antarctic surface;

Table 1 .
CYCLOPS model parameter values defining the ensemble of 13 500 simulations * .

Table 2 .
Uncertainty inputs into the Monte Carlo simulations to calculate δ 11 B. The sources of uncertainty are also added.All uncertainty estimates are 2σ .
and based on the GEOCARB model, where the circularity prob-