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

Rosanna Greenop, Mathis Hain, Sindia Sosdian, Kevin Oliver, Philip Goodwin, Thomas B. Chalk, Caroline Lear, Paul Wilson, Gavin Foster

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

► To cite this version:

Rosanna Greenop, Mathis Hain, Sindia Sosdian, Kevin Oliver, Philip Goodwin, et al.. A record of Neogene seawater δ11B reconstructed from paired δ11B analyses on benthic and planktic foraminifera. Climate of the Past, 2017, 13 (2), pp.149-170. 10.5194/cp-13-149-2017. hal-04700470

HAL Id: hal-04700470

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

Submitted on 19 Sep 2024

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

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

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

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

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

Rosanna Greenop<fnref n="1" />,<fnref n="2" />, Mathis P. Hain<fnref n="1" />, Sindia M. Sosdian<fnref n="3" />, Kevin I. C. Oliver<fnref n="1" />, Philip Goodwin<fnref n="1" />, Thomas B. Chalk<fnref n="1" />,<fnref n="4" />,
Caroline H. Lear<fnref n="3" />, Paul A. Wilson<fnref n="1" />, and Gavin L. Foster<fnref n="1" />

*Correspondence to:* Rosanna Greenop (rg200@st-andrews.ac.uk)

Received: 27 November 2015 – Discussion started: 26 January 2016

Revised: 4 December 2016 – Accepted: 9 January 2017 – Published: 24 February 2017

**Abstract.** The boron isotope composition (δ¹¹B) of foraminiferal calcite reflects the pH and the boron isotope composition of the seawater the foraminifer grew in. For pH reconstructions, the δ¹¹B of seawater must therefore be known, but information on this parameter is limited. Here we reconstruct Neogene seawater δ¹¹B based on the δ¹¹B difference between paired measurements of planktic and benthic foraminifera and an estimate of the coeval water column pH gradient from their δ¹³C values. Carbon cycle model simulations underscore that the ΔpH–Δδ¹³C relationship is relatively insensitive to ocean and carbon cycle changes, validating our approach. Our reconstructions suggest that δ¹¹Bsw was ~37.5‰ during the early and middle Miocene (roughly 23–12 Ma) and rapidly increased during the late Miocene (between 12 and 5 Ma) towards the modern value of 39.61‰. Strikingly, this pattern is similar to the evolution of the seawater isotope composition of Mg, Li and Ca, suggesting a common forcing mechanism. Based on the observed direction of change, we hypothesize that an increase in secondary mineral formation during continental weathering affected the isotope composition of riverine input to the ocean since 14 Ma.

## 1 Introduction

Key to determining the relationship between CO₂ and climate in the geological past is the calculation of reliable estimates of absolute CO₂ through time. In recent years the boron isotope composition (δ¹¹B) of foraminiferal calcite has become a high-profile tool for reconstructing CO₂ beyond the last 800ky 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; Henahan et al., 2013; Greenop et al., 2014; Martínez-Botí et al., 2015a). Yet long-term change in the boron isotope composition of seawater (δ¹¹Bsw) is currently poorly constrained and represents a major source of the uncertainty associated with δ¹¹B-determined CO₂ 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 δ¹¹Bsw to be an important factor in CO₂ 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-

<footnote n="1">¹Ocean and Earth Science, National Oceanography Centre Southampton, University of Southampton, Waterfront Campus, European Way, Southampton SO14 3ZH, UK</footnote>

<footnote n="2">²School of Geography & Geosciences, Irvine Building, University of St Andrews, North Street, St Andrews, KY16 9AL, UK</footnote>

<footnote n="3">³School of Earth & Ocean Sciences, Cardiff University, Cardiff, CF10 3AT, UK</footnote>

<footnote n="4">⁴Department of Physical Oceanography, Woods Hole Oceanographic Institution, Woods Hole, Massachusetts, USA</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The oceanic boron cycle. Fluxes are from Lemarchand et al. (2000) and Park and Schlesinger (2002). Isotopic compositions are from Lemarchand et al. (2000), Foster et al. (2010) and references therein.</figcaption></figure>

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 (δ¹¹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 (Spackach and Edmond, 1987) and co-precipitation into carbonates (Hemming and Hanson, 1992). In the case of all three outputs, the light ¹⁰B isotope is preferentially removed relative to ¹¹B, such that the seawater ¹¹B/¹⁰B ratio (δ¹¹B$_{sw}$, 39.61‰) is significantly greater than that of the cumulative inputs (δ¹¹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 δ¹¹B$_{sw}$ record we present.

Unlike many other isotope systems (e.g. δ⁷Li$_{sw}$, δ²⁶Mg$_{sw}$, δ⁴⁴/⁴⁰Ca$_{sw}$, ⁸⁷Sr/⁸⁶Sr) to date, no direct archive has been documented for δ¹¹B$_{sw}$. This is a result of the pH-dependent boron speciation in seawater upon which the δ¹¹B–pH proxy is based (Hemming and Hanson, 1992), which imparts a pH dependency on the δ¹¹B of all marine precipitates so far ex-

amined. Empirical reconstructions of δ¹¹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) δ¹¹B analysis of halites (Paris et al., 2010), (3) measurements of benthic foraminiferal δ¹¹B coupled to various assumptions about past changes in ocean pH (Raitzsch and Hönisch, 2013), and (4) measurements of δ¹¹B in surface- and thermocline-dwelling foraminifera coupled with additional information on the pH gradient of the surface ocean (Palmer et al., 1998; Pearson and Palmer 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 δ¹¹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 δ¹¹B$_{sw}$ (39.7‰) with no apparent fractionation, measurement of δ¹¹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, δ¹¹B$_{sw}$ is calculated from globally distributed benthic δ¹¹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
<figure><img src="image_5.png" /></figure>

**Figure 2.** A compilation of published δ¹¹Bsw 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 δ¹¹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 δ¹¹B measurements into δ¹¹Bsw. This approach yields broadly similar results to geochemical modelling (Fig. 2).

Approach 4 exploits the non-linear relationship between δ¹¹B and pH alongside estimated pH gradients in the ocean to constrain δ¹¹Bsw (Palmer et al., 1998; Pearson and Palmer 1999, 2000), and this is the basis of the approach used in this study. The advantage of this method is that δ¹¹Bsw can be reconstructed empirically without relying on a priori absolute-pH constraints. The non-linear relationship between δ¹¹B and pH means that the pH difference between two δ¹¹B data points varies as a function of δ¹¹Bsw (Fig. 3). Consequently, if the size of the pH gradient can be estimated, then there is only one δ¹¹Bsw value that is consistent with the foraminiferal δ¹¹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₂ to account for changes in δ¹¹Bsw (determined using δ¹¹B in surface- and thermocline-dwelling foraminifera) (Fig. 2) (Palmer et al., 1998; Pearson and Palmer 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-layer-dwelling species and thermocline dweller to calculate a value for δ¹¹Bsw (Pearson and Palmer, 1999). The resulting record

suggests that δ¹¹Bsw 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 δ¹¹B gradients to calculate δ¹¹Bsw, was recently applied to the middle Miocene where it yielded a δ¹¹Bsw of 37.6 ±0.4 ‰ (Foster et al., 2012). A further modification to the method of Pearson and Palmer (1999) was also proposed in that study wherein δ¹³C in foraminiferal calcite was used to estimate the surface-to-deep pH gradient (Foster et al., 2012). Here, we reconstruct δ¹¹Bsw 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 δ¹³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 δ¹¹Bsw. The implications of our work for understanding the evolution of Neogene ocean pH and atmospheric pCO₂ will be documented elsewhere.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Schematic diagram showing the change in pH gradient with a 3‰ change in δ¹¹B for δ¹¹Bsw of (a) 39.6 and (b) 37.5‰. Arrows highlight the different pH gradients. Note how a δ¹¹B difference of 3‰ is translated into different pH gradients depending on the δ¹¹Bsw. Calculated using B_T = 432.6 µmol kg⁻¹ (Lee et al., 2010) and α_B = 1.0272 (Klochko et al., 2006). (c) The pH change for a δ¹¹B change of 3‰ at a range of different δ¹¹Bsw.</figcaption></figure>

# 2 Methods

## 2.1 Site locations and age models

Foraminifera from four sites are used to construct the planktic-benthic δ¹¹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).

## 2.2 Boron isotope analysis and pH calculation

The boron isotope measurements (expressed in delta notation as δ¹¹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 sample-standard 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 δ¹¹B of CaCO₃ and pH is very closely approximated by the following equation:

$$ \text{pH} = pK_B^* - \log \left( - \frac{\delta^{11}B_{\text{SW}} - \delta^{11}B_{\text{CaCO}_3}}{\delta^{11}B_{\text{SW}} - \alpha_B \cdot \delta^{11}B_{\text{CaCO}_3} - 1000 (\alpha_B - 1)} \right), \quad (1) $$

where $pK_B^*$ is the equilibrium constant, dependent on salinity, temperature, pressure and seawater major ion composition (i.e. [Ca] and [Mg]); $\alpha_B$ is the fractionation factor be-
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Map of study sites and mean annual air-sea disequilibrium with respect to pCO₂. 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).</figcaption></figure>

tween the two boron species; and δ¹¹Bsw 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 δ¹¹B of foraminifera correlates well with pH and hence [CO₂]aq, the δ¹¹Bcalcite is often not exactly equal to δ¹¹Bborate (Sanyal et al., 2001; Foster, 2008; Henehan et al., 2013). The planktic species used to construct the benthic-planktic 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 praebullooides* 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. praebullooides* 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 δ¹¹B of both *Cibicidoides* species shows no offset from the theoretical δ¹¹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 δ¹¹Bcalcite, temperature, salinity, water depth (pressure) and seawater major ion composition are also needed to calculate pH from δ¹¹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 δ¹¹B sample (with a conservative 2σ uncertainty of 2 °C). Adjustments were made for changes in Mg/Casw using the records of Horita et al. (2002) and Brennan et al. (2013). We corrected for changes in dependence on Mg/Casw 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/Casw; Delaney et al., 1985; Hasiuk and Lohmann, 2010; Evans and Müller, 2012) using the equation

$$ \text{Mg/C}_{\text{sw,c}} = (\text{Mg/C}_{\text{sw,a}}/\text{Mg/C}_{\text{sw,m}})^{0.41}, \quad (2) $$

where Mg/Csw,c is the correction factor applied to the temperature equation for changing Mg/Casw, Mg/Csw,a is the estimated Mg/Casw for the age of the sample and Mg/Casw,m is modern Mg/Casw. Temperature is then calculated using the generic planktic foraminifera calibration of Anand et al. (2003) and including a correction factor for Mg/Casw.

$$ \text{Temperature} = \ln(\text{Mg/C}_{\text{test}}/(0.38 \times \text{Mg/C}_{\text{sw,c}}))/0.09 \quad (3) $$

Mg/Ca analysis was conducted on a small aliquot of the sample dissolved for isotope analysis at the University of
Southampton using a Thermo 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 of Cramer et al. (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 δ¹¹B~sw~ is dominated by the uncertainty in the δ¹¹B measurement and the estimate of the pH gradient.

The majority of the δ¹³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 δ¹³C from Site 999 were taken from the nearest sample in Haug and Tiedemann (1998). In almost all cases δ¹³C was analysed on the same foraminiferal species as δ¹¹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 δ¹³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 δ¹³C of dissolved inorganic carbon (DIC).

## 2.3 Carbon isotopes as a proxy for vertical ocean pH gradient

The use of δ¹³C in foraminiferal calcite to estimate the surface-to-deep pH gradient requires knowledge of the slope of the pH-δ¹³C relationship in the past. In this section we briefly outline the main factors that contribute to the pH-δ¹³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-δ¹³C organic carbon via the soft-tissue component of the biological pump leads to a broad correlation between δ¹³C, [CO₂⁻] and macronutrients in the ocean (e.g.

Hain et al., 2014a). The remineralization of this organic matter decreases δ¹³C and titrates [CO₃²⁻], thereby reducing pH, while nutrient concentrations are increased. In waters that have experienced more soft tissue remineralization, both pH and δ¹³C will be lower (Fig. 5a, b). This is the dominant reason for the positive slope between δ¹³C and pH in the modern ocean (e.g. Foster et al., 2012; Fig. 5c).

Another significant factor affecting the spatial distribution of both δ¹³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₃²⁻] and pH (Goodwin and Lauderdale, 2013), while warmer waters have relatively low equilibrium δ¹³C values (Lynch-Stieglitz et al., 1995). This means that a spatial gradient in temperature acts to drive δ¹³C and pH in opposite directions: warmer waters tend to have higher pH but lower δ¹³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 δ¹³C alone. In other words, when using δ¹³C differences to estimate the pH gradient between the warm low-latitude surface waters and cold deep waters, the appropriate ΔpH-Δδ¹³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-Δδ¹³C rather than the slope of the overall pH-δ¹³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 δ¹³C and pH (and [CO₃²⁻]). However, other processes will have had a minor effect on either pH or δ¹³C. For instance, the dissolution of CaCO₃ shells increases [CO₃²⁻] and pH (Broecker and Peng, 1982), but does not significantly affect δ¹³C (Zeebe and Wolf-Gladrow, 2001). Moreover, the long timescale of air-sea isotopic equilibration of CO₂ combined with kinetic isotope fractionation during net carbon transfer is an important factor in setting the distribution of δ¹³C on a global ocean scale (Galbraith et al., 2015; Lynch-Stieglitz et al., 1995). Meanwhile, the effect of CO₂ disequilibrium on [CO₃²⁻] and pH is modest (Goodwin and Lauderdale, 2013).

## 2.4 Modelling the pH-δ¹³C relationship

After correcting for the shift in δ¹³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 δ¹³CDIC, 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-
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Latitudinal cross section through the Atlantic showing (a) pH variations and (b) the δ¹³C composition. Data are plotted using Ocean Data View (Schlitzer 2016). pH data are from the CARINA data set (CARINA group, 2009) and the δ¹³C data are from the GLODAP data compilation (Key et al., 2004). (c) pH and δ¹³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).</figcaption></figure>

scribed above (i.e. warm surface to cold bottom-water ΔpH–Δδ¹³C). Here we use an ensemble approach with two independent carbon cycle models to investigate changes in the ΔpH–Δδ¹³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–Δδ¹³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 δ¹³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₃ 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₂. This enables us to assess which processes, if any, are capable of altering the oceanic relationship between ΔpH and Δδ¹³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 δ¹³C on surface disequilibrium (Galbraith et al., 2015); (4) inorganic-to-organic carbon rain ratio, controlling the relationship between DIC and alkalinity distributions; (5) "Antarctic shelf FWF", a freshwater flux adjustment (always switched off in control experiments with GENIE) 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-
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 δ¹³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-Δδ¹³C relationship we use the CYCLOPS bio-geochemical 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., 2010, 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) and (b) reconstructed secular changes in seawater [Ca] (calcium concentration), carbonate compensation depth (CCD), weathering and atmospheric CO₂ (Table 1). 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) iron-driven 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; (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.9 km via changes in the weathering flux, as per sedimentological evidence (Pälike et al., 2012); and (7) atmospheric CO₂ 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₃ 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₃ burial for the final solution (unlike the GENIE simulations, CaCO₃ 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₂ (Beerling and Royer, 2011) and records

of glacial-interglacial circulation change (Curry and Oppo, 2005).

## 2.5 Assessing uncertainty

δ¹¹Bsw 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 = 10000) 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⁻¹), [Mg] (±4.5 mmol kg⁻¹), δ¹¹Bplanktic (±0.15–0.42‰) and δ¹¹Bbenthic (±0.21–0.61‰). For the estimate of the surface-to-sea floor pH gradient, we use the central value of the ΔpH-Δδ¹³C relationship diagnosed from our CYCLOPS and GENIE sensitivity experiments (i.e. 0.175‰⁻¹; 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‰⁻¹ ΔpH / Δδ¹³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 δ¹³C output (see Sect. 3.2 below). The uncertainty in the δ¹¹B measurements is calculated from the long-term reproducibility of Japanese Geological Survey Porites coral standard (JCP; δ¹¹B = 24.3‰) at the University of Southampton using the equations

$$2\sigma = 2.25\exp^{-23.01^{[11B]}} + 0.28\exp^{-0.64^{[11B]}} \quad (4)$$

$$2\sigma = 33450\exp^{-168.2^{[11B]}} + 0.311\exp^{-1.477^{[11B]}}, \quad (5)$$

where [¹¹B] is the intensity of ¹¹B signal in volts and Eqs. (4) and (5), used with 10¹¹Ω and 10¹²Ω resistors respectively.

From the 10 000 Monte Carlo ensemble solutions of our 22 benthic-planktic pairs, we construct 10 000 randomized records of δ¹¹Bsw as a function of time. Each of these randomized δ¹¹Bsw 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 δ¹¹Bsw 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
<table><caption>**Table 1.** CYCLOPS model parameter values defining the ensemble of 13 500 simulations*.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Values assumed</th></tr></thead><tbody><tr><td>PAZ surface phosphate**</td><td>Unutilized polar nutrient</td><td>1, 1.25, 1.5, 1.75, 2 µM</td></tr><tr><td>PAZ vertical exchange**</td><td>Bottom water formation</td><td>2, 7.75, 13.5, 19.25, 25 Sv</td></tr><tr><td>SAZ surface phosphate**</td><td>Unutilized polar nutrient</td><td>0.7, 0.825, 0.95, 1.075, 1.2 µM</td></tr><tr><td>AMOC circulation scheme**</td><td>Deep vs. shallow overturning</td><td>NADW, GNAIW</td></tr><tr><td>Representative time slice***</td><td>Age ([Ca²⁺]/CCD); calcium set outright; CCD set via riverine CaCO₃ flux using inverse scheme</td><td>0 My (10.6 mM, 4.65 km), 9 My (12.89 mM, 4.4 km), 11 My (13.33 mM, 4.9 km), 16 My (14.28 mM, 4.7 km), 18 My (14.57 mM, 4.25 km), 20 My (14.86 mM, 4.7 km)</td></tr><tr><td>Atm. CO₂****</td><td>Set via silicate weatherability</td><td>200, 300, 400, 500, 600, 700, 800, 900, 1000 ppm</td></tr></tbody></table>

<table><caption>**Table 2.** Uncertainty inputs into the Monte Carlo simulations to calculate δ¹¹B. The sources of uncertainty are also added. All uncertainty estimates are 2σ.</caption><thead><tr><th>Input parameter</th><th>Uncertainty applied</th><th>Source of uncertainty estimate</th></tr></thead><tbody><tr><td>Surface-to-sea floor ΔpH</td><td>Uniform ±0.05 pH units</td><td>Plausible range of ΔpH / Δδ¹³C in CYCLOPS and GENIE sensitivity tests; prediction error of linear ΔpH / Δδ¹³C regression in GENIE</td></tr><tr><td>δ¹¹B measurement</td><td>0.15–0.61‰</td><td>Long-term external reproducibility</td></tr><tr><td>Temperature</td><td>±2 °C</td><td>Uncertainty in the Mg / Ca measurement and Mg / Ca temperature calibration</td></tr><tr><td>Salinity</td><td>±2 psu</td><td>In the absence of a salinity proxy this uncertainty is applied to cover variations through time.</td></tr><tr><td>Seawater [Mg]</td><td>±4.5 mmol kg⁻¹</td><td>Following Horita et al. (2002)</td></tr><tr><td>Seawater [Ca]</td><td>±4.5 mmol kg⁻¹</td><td>Following Horita et al. (2002)</td></tr></tbody></table>

input parameters), while not smoothing away the underlying signal. As a result of anomalously low δ¹¹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 δ¹¹Bsw curves represents the combination of the compounded, propagated uncertainties of the various inputs (i.e. Monte Carlo sampling) with the additional constraint of gradual δ¹¹Bsw 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 δ¹¹Bsw reconstruction were evaluated from the ensemble of smoothed δ¹¹Bsw records. Generally, for any given benthic–planktic pair the resulting δ¹¹Bsw 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).

# 3 Results and discussion

## 3.1 δ¹¹B benthic and planktic data

Surface and deep-ocean δ¹¹B broadly show a similar but inverse pattern to δ¹³C and temperature throughout the Neogene (Fig. 6). The δ¹¹B benthic record decreases from

<footnote>* 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₃ cycle relevant to our study interval. **** These atmospheric CO₂ levels are chosen to span a range wider than expected for the study interval. Following the silicate-weathering-feedback paradigm, long-term CO₂ is fully determined by the balance of geologic CO₂ sources and silicate weathering, whereby faster acting processes of the open system CaCO₃ cycle compensate relative to that CO₂ level. All else equal, high CO₂ 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.</footnote>
<table><caption>**Table 3.** The average δ¹¹B, δ²⁶Mg, δ⁴⁴/⁴⁰Ca and δ⁷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 on Lemarchand et al. (2000) and references therein, lithium on Misra 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 δ²⁶Mgsw and δ¹¹Bsw are from Foster et al. (2010) and δ⁷Lisw is from Tomascak (2004). The δ⁴⁴/⁴⁰Ca presented here was measured relative to seawater and hence seawater has a δ⁴⁴/⁴⁰Ca sw of zero per mill by definition.</caption><thead><tr><th>Sources</th><th colspan="4">Isotopic ratio</th></tr></thead><tbody><tr><td><strong>Oceanic inputs</strong></td><td>δ<sup>11</sup>B<sub>sw</sub><br>39.61 ‰</td><td>δ<sup>7</sup>Li<sub>sw</sub><br>31 ‰</td><td>δ<sup>26</sup>Mg<sub>sw</sub><br>-0.83 ‰</td><td>δ<sup>44/40</sup>Ca<sub>sw</sub><br>0 ‰</td></tr><tr><td>Input from<br>hydrothermal</td><td>6.5<sup>a</sup></td><td>8.3<sup>b</sup></td><td>n/a</td><td>-0.96<sup>h</sup></td></tr><tr><td>Fluid from<br>accretionary prisms</td><td>25<sup>a</sup></td><td>15<sup>b</sup></td><td>n/a</td><td>n/a</td></tr><tr><td>Riverine inputs</td><td>10<sup>a</sup></td><td>23<sup>b</sup></td><td>-1.09<sup>d</sup></td><td>-1.28<sup>h</sup></td></tr><tr><td>Groundwater</td><td>n/a</td><td>n/a</td><td>-0.82<sup>d</sup></td><td>-1.02<sup>i</sup></td></tr><tr><td><strong>Outputs</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Precipitation into<br>carbonates</td><td>20<sup>a</sup></td><td>29<sup>c</sup></td><td>-3.5<sup>d,e,f</sup></td><td>-1.15<sup>h,j</sup></td></tr><tr><td>Ocean crust<br>alteration</td><td>4<sup>a</sup></td><td>15<sup>b</sup></td><td>-0.83<sup>d,g</sup></td><td>-1.2<sup>h</sup></td></tr><tr><td>Absorption onto<br>sediment</td><td>15<sup>a</sup></td><td>15<sup>b</sup></td><td>??</td><td>n/a</td></tr></tbody></table>

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.

~ 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 δ¹¹B are in keeping with previously published Neogene δ¹¹B benthic records measured using NTIMS (Raitzsch and Hönisch, 2013). This suggests that our deep-water δ¹¹B record is representative of large-scale pH changes in the global ocean. While the surface δ¹¹Bplanktic remained relatively constant between 24 and 11 Ma at ~ 16‰, there has been a significant increase in δ¹¹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 δ¹³C DIC both

broadly decrease through the Neogene and appear to co-vary on shorter timescales (Fig. 6e, f).

## 3.2 The relationship between δ¹³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 δ¹³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–δ¹³C relationship is
<figure><img src="image_9.png" /></figure>

**Figure 6.** δ¹¹B<sub>planktic</sub>, temperature and δ¹³C<sub>DIC</sub> estimates for the surface and deep ocean through the last 23 million years. (a) δ¹¹B<sub>planktic</sub> surface; (b) δ¹¹B<sub>borate</sub> 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 Raitzsch 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) δ¹³C<sub>DIC</sub> surface record; and (f) δ¹³C<sub>DIC</sub> 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. praebuloides*, dark blue *Cibicidoides* wuellerstorfi and light blue *Cibicidoides* *mundulus*. The two benthic–planktic pairs that were removed prior to smoothing are highlighted with arrows.

slightly greater than 0.2‰⁻¹ for the full three-dimensional

data regression (not shown) and about 0.175‰⁻¹ for the warm surface–cold deep ΔpH–Δδ¹³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‰⁻¹ 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 δ¹³C (see below).

In our CYCLOPS model ensemble, the central values of the slopes of the full three-dimensional pH–δ¹³C regressions and of the warm surface-to-cold deep ΔpH–Δδ¹³C are 0.2047‰⁻¹ (1σ of 0.0196‰⁻¹; Fig. 8a) and 0.1797‰⁻¹ (1σ of 0.0213‰⁻¹; Fig. 8b) respectively. If we restrict our analysis of the CYCLOPS ensemble to only the Atlantic basin warm surface-to-cold deep ΔpH–Δδ¹³C, where most of our samples come from, we find a relationship of only 0.1655‰⁻¹ (1σ of 0.0192‰⁻¹; 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–Δδ¹³C central value of 0.175‰⁻¹ to estimate the surface–sea floor pH difference from the planktic–benthic foraminifera δ¹³C difference. To account for our ignorance as to the accurate value of ΔpH–Δδ¹³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–Δδ¹³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 planktic–planktic gradients are employed, the plausible ΔpH–Δδ¹³C range should be significantly higher than applied here to account for the relatively lower temperature difference. Based on the appropriate ΔpH–Δδ¹³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–Δδ¹³C relationship, we point to changes in that relationship that arise in our GENIE 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–Δδ¹³C remaining inside that range. While CYCLOPS simulates the balance between weathering and CaCO₃ 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-
<figure><img src="image_10.png" /></figure>

**Figure 7.** The output from GENIE sensitivity analysis showing the warm surface-to-cold deep ΔpH–Δδ¹³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–Δδ¹³C relationship (where the colours reflect the CO₂ level of each experiment) and red open circles show the RMSE. The green stars are the ΔpH–Δδ¹³C relationship for the control experiment conducted at 292.67 ppm CO₂. 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.

val and likely for most of the Phanerozoic. The differing outputs from CYCLOPS and GENIE in the DIC and ALK experiments shows that ΔpH–Δδ¹³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–Δδ¹³C for at least the Neogene, if not the entire Cenozoic, outside of extreme events such as the Palaeocene–Eocene Thermal Maximum.

### 3.3 δ¹¹Bsw record through the Neogene

Using input parameter uncertainties as described in section 2.5 yields individual Monte Carlo member δ¹¹Bsw 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. However, for each planktic–benthic time point most individual Monte Carlo δ¹¹Bsw estimates fall into a much narrower central range (~1 to 4‰; thick black line showing interquartile range in
<figure><img src="image_11.png" /></figure>

**Figure 8.** The output from sensitivity analysis of the relationship between pH gradient and δ¹³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 / Δδ¹¹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.

Fig. 10a). The δ¹¹B<sub>sw</sub> for Plio-Pleistocene time points clus-
ter 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 ob-
served clustering is strong evidence for an underlying long-
term signal in our data, albeit one that is obscured by the un-
certainties involved in our individual δ¹¹B<sub>sw</sub> estimates. The
same long-term signal is also evident when pooling the indi-
vidual Monte Carlo member δ¹¹B<sub>sw</sub> estimates into 8-million-
year 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 δ¹¹Bsw appears to also have been significantly lower than
present day during the early Miocene.

3.3.1 Data smoothing

The ~1 to 4‰ likely ranges for δ¹¹B<sub>sw</sub> would seem to
be rather disappointing given the goal to constrain δ¹¹B<sub>sw</sub>
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 con-
siderations that δ¹¹B<sub>sw</sub> of seawater should not change by
more than ~0.1‰ per million years (Lemarchand et al.,
<figure><img src="image_12.png" /></figure>

**Figure 9.** The pH gradient between surface and deep water through time calculated from the $\delta^{13}$C gradient and using a flat probability derived from the low-latitude ensemble regressions from the CYCLOPS 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 $\delta^{11}$B$_{sw}$ estimate is equally likely and the smoothing should therefore target randomly selected individual Monte Carlo $\delta^{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 $\delta^{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 $\delta^{11}$B$_{sw}$ change to be dictated by the data, with only the most extreme shifts in $\delta^{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 $\delta^{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).

## 3.4 Comparison to other $\delta^{11}$B$_{sw}$ records

The comparison of our new $\delta^{11}$B$_{sw}$ record to those previously published reveals that despite the differences in methodology, the general trends in the records show excel-

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

**Figure 10.** The calculated $\delta^{11}$B$_{sw}$ from the benthic–planktic $\delta^{11}$B pairs using a pH gradient derived from $\delta^{13}$C. The uncertainty on each data point is determined using a Monte Carlo approach, including uncertainties in temperature, salinity, $\delta^{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 $\delta^{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 $\delta^{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 $\delta^{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 $\delta^{11}$B$_{sw}$ evolution is an increase through time from the middle Miocene to the Plio-Pleistocene (Fig. 11). While the model-based $\delta^{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-
<figure><img src="image_14.png" /></figure>

**Figure 11.** The δ¹¹Bsw curve calculated using the variable pH gradient derived from δ¹³C. The median (red line), 68% (dark red band) and 95% (light red band) confidence intervals are plotted. Plotted with a compilation of published δ¹¹Bsw 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 δ¹¹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 δ¹¹Bsw 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 reconstruc-
tion. 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 δ¹¹Bsw rise following the middle Miocene
is in almost perfect agreement. We are encouraged by these
agreements resulting from approaches based on very differ-
ent underlying assumptions and techniques, which we take
as indication for an emerging consensus view of δ¹¹Bsw evo-
lution over the last 25 My and as a pathway towards recon-
structing δ¹¹Bsw further back in time.

The record by Pearson and Palmer (2000) is well corre-
lated to our reconstruction, but especially during the early
Miocene there is a notable ~0.5‰ offset (Fig. 11). This dis-
crepancy could be due to a number of factors. Firstly, the
applicability of this δ¹¹Bsw record (derived from δ¹¹B data
measured using NTIMS) to δ¹¹B records generated using the
MC-ICPMS is uncertain (Foster et al., 2013). In addition, this
δ¹¹Bsw record is determined using a fractionation factor of
1.0194 (Kakihana et al., 1977), whereas recent experimen-
tal 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 under-
standing of the δ¹¹B difference between species–size frac-
tions (Foster, 2008; Henehan et al., 2013), the mixed species
and size fractions used to make the δ¹¹B measurements in

that study may have introduced some additional uncertainty
in the reconstructed δ¹¹Bsw. Conversely, there is a substan-
tial spread between our three time points during the earli-
est 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 bi-
ased during the early Miocene.

The δ¹¹Bsw record calculated using benthic δ¹¹B and as-
sumed deep-ocean pH changes (Raitzsch and Hönisch, 2013)
is also rather similar to our δ¹¹Bsw reconstruction. The dis-
crepancy 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 assump-
tion 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₂ 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) and
based on the GEOCARB model, where the circularity prob-
<figure><img src="image_15.png" /></figure>

**Figure 12.** (a) The $\delta^{11}$B$_{sw}$ curve from this study plotted with other trace element isotopic records. On the $\delta^{11}$B$_{sw}$ panel the darker shaded area highlights the 68% confidence interval and the lighter interval highlights the 95% confidence interval, (b) $\delta^{26}$Mg$_{sw}$ record from Pogge von Strandmann et al. (2014) (error bars are $\pm$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) $\delta^{44/40}$Ca$_{sw}$ record from Griffith et al. (2008) (error bars show $2\sigma$ uncertainty) and (d) $\delta^{7}$Li$_{sw}$ record from Misra and Froelich (2012) (error bars show $2\sigma$ uncertainty). Blue dashed lines show middle Miocene values and red dashed lines highlight the modern.

This study 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 $\delta^{11}$B$_{sw}$ reconstructions.

Our new $\delta^{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 $\delta^{11}$B$_{sw}$ signal in the fluid inclusions (Paris et al., 2010) is most likely a combination of the $\delta^{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-$\delta^{11}$B isotope values such that the fluid inclusion record likely provides a lower limit for the $\delta^{11}$B$_{sw}$ through time (i.e. $\delta^{11}$B$_{sw}$ is heavier than the halite fluid inclusions of Paris et al., 2010). For this halite record to be interpreted directly as $\delta^{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.

## 3.5 Common controls on the seawater isotopic ratios of B, Mg, Ca and Li

Our new record of $\delta^{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 ($\delta^{7}$Li$_{sw}$; Misra and Froelich, 2012) and the calcium isotopic composition of seawater as recorded in marine barites ($\delta^{44/40}$Ca$_{sw}$; Griffith et al., 2008) both increase through the Neogene, whereas the magnesium isotopic composition of seawater ($\delta^{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 $\delta^{11}$B$_{sw}$, we interpolate and cross plot $\delta^{11}$B$_{sw}$ and the $\delta^{7}$Li$_{sw}$, $\delta^{44/40}$Ca$_{sw}$ and $\delta^{26}$Mg$_{sw}$ records. This analysis suggests that the isotopic composition of $\delta^{11}$B$_{sw}$, $\delta^{7}$Li$_{sw}$, $\delta^{26}$Mg$_{sw}$ and $\delta^{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 $\delta^{11}$B$_{sw}$, $\delta^{7}$Li$_{sw}$, $\delta^{26}$Mg$_{sw}$ and $\delta^{44/40}$Ca$_{sw}$, given that the factors fractionating these stable isotopic systems are similar (see below).

To better constrain the controls on $\delta^{11}$B$_{sw}$, $\delta^{7}$Li$_{sw}$, $\delta^{26}$Mg$_{sw}$ and $\delta^{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,
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Crossplots of the records of δ¹¹B~sw~ using a variable pH gradient derived from δ¹³C (error bars show 2σ uncertainty), with δ⁴⁴/⁴⁰Ca~sw~ from Griffith et al. (2008) (error bars show 2σ uncertainty), δ⁷Li~sw~ from Misra and Froelich (2012) (error bars show 2σ uncertainty) and δ²⁶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.</figcaption></figure>

2012). The three dominant controls on magnesium concen-
tration 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 bal-
ance between riverine and hydrothermal inputs and removal
through CaCO₃ deposition and alteration of oceanic crust
(Fantle and Tipper, 2014; Griffith et al., 2008). Dolomitiza-
tion has also been cited as playing a potential role in con-
trolling δ⁴⁴/⁴⁰Caₗ, 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 δ¹¹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 δ¹¹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 ap-
pears to have limited influence on the δ⁷Li composition of
rivers, and riverine δ⁷Li varies primarily with weathering in-
tensity (Kisakūrek et al., 2005; Millot et al., 2010). The river-
ine input of calcium to the oceans is controlled by the compo-
sition of the primary continental crust (dominated by carbon-
ate weathering) and a recycled component, although the rel-
ative influence of these two processes is not well understood
(Tipper et al., 2006a). In addition, vegetation may also play
a significant role in the δ⁴⁴/⁴⁰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 δ⁷Liₗ has already been attributed to an increase in
the δ⁷Liₗ composition of the riverine input (Hathorne and
James, 2006; Misra and Froelich, 2012). The causes of the
shift in δ⁷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 up-
lift (Misra and Froelich, 2012; Li and West, 2014), (2) a re-
duction in weathering intensity (Hathorne and James, 2006;
Froelich and Misra, 2014; Wanner et al., 2014), (3) an in-
crease in silicate weathering rate (Liu et al., 2015), (4) an in-
crease in the formation of floodplains and the increased for-
mation of secondary minerals (Pogge von Strandmann and
Henderson, 2014), and (5) a climatic control on soil produc-
tion 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 miner-
als on land is also consistent across Mg, Ca and B isotope
systems. For instance, clay minerals are preferentially en-
riched 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), ade-
quately explaining the inverse relationship between $\delta^{11}$B$_{sw}$
and $\delta^{26}$Mg$_{sw}$. Consequently, the increased formation or re-
tention on land of secondary minerals would alter the iso-
topic composition of the riverine input to the ocean in the
correct direction to explain the trends in all four isotope sys-
tems through the late Neogene (Fig. 13). While the relation-
ships 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 iso-
tope systems (Tipper et al., 2006b; Fantle and Tipper, 2014).
Consequently, a future model of seawater chemistry evolu-
tion 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 chem-
istry. Nonetheless, given the similarities between the geo-
chemical cycles of B and Li, and despite the large differ-
ence in residence time (Li = 1 million years, B = 11–17 mil-
lion years), the correlation between these two records is com-
pelling and would no doubt benefit from additional study.

## 4 Conclusions

Here we present a new $\delta^{11}$B$_{sw}$ record for the Neogene (0–23 Ma) based on paired planktic–benthic $\delta^{11}$B measurements. Our new record suggests that $\delta^{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 pub-
lished data-based reconstructions capture the first-order late
Miocene $\delta^{11}$B$_{sw}$ rise suggests that consensus is building for
the $\delta^{11}$B$_{sw}$ evolution through the Neogene. This emerging
view on $\delta^{11}$B$_{sw}$ change provides a vital constraint required to
quantitatively reconstruct Neogene ocean pH, ocean carbon
chemistry and atmospheric CO$_2$ using the $\delta^{11}$B–pH proxy.
When our new $\delta^{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 $\delta^{26}$Mg$_{sw}$, $\delta^{44/40}$Ca$_{sw}$, $\delta^{11}$B$_{sw}$ and $\delta^{7}$Li$_{sw}$
records are all consistent with an increase in secondary min-
eral formation through time. While a more quantitative treat-
ment of these multiple stable isotope systems is required, the
$\delta^{11}$B$_{sw}$ record presented here provides additional constraints
on the processes responsible for the evolution of ocean chem-
istry through time.

## 5 Data availability

Data are available at this link: https://doi.pangaea.de/10.
1594/PANGAEA.871900.

**The Supplement related to this article is available online**
**at doi:10.5194/cp-13-149-2017-supplement.**

**Competing interests.** The authors declare that they have no con-
flict of interest.

**Acknowledgements.** This work used samples provided by (I)ODP, which is sponsored by the US National Science Foundation, and participating countries under the management of Joint Oceano-graphic Institutions, Inc. We thank Walter Hale and Alex Wuelbers of the Bremen Core Repository for their kind assistance. The work was supported by NERC grants NE/I006176/1 (Gavin L. Foster and Caroline H. Lear), NE/H006273/1 (Gavin L. Foster), NE/I006168/1 and NE/K014137/1 and a Royal Society Research Merit Award (Paul A. Wilson), a NERC Independent Research Fellowship NE/K00901X/1 (Mathis P. Hain) and a NERC studentship (Rosanna Greenop). Matthew Cooper, J. Andy Milton and the B-team are acknowledged for their assistance in the laboratory. We thank the two anonymous reviewers and Philip Pogge von Strandmann for their helpful suggestions that improved the paper.

Edited by: A. Sluijs

Reviewed by: P. A. E. Pogge von Strandmann and two anonymous
referees
## References

Al-Rousan, S., Pätzold, J., Al-Moghrabi, S., and Wefer, G.: Invasion of anthropogenic CO₂ recorded in planktonic foraminifera from the northern Gulf of Aqaba, Int. J. Earth Sci., 93, 1066–1076, 2004.

Anagnostou, E., John, E. H., Edgar, K. M., Foster, G. L., Ridgewell, A., Inglis, G. N., Pancost, R. D., Lunt, D. J., and Pearson, P. N.: Changing atmospheric CO₂ concentration was the primary driver of early Cenozoic climate, Nature, 533, 380–384, 2016.

Anand, P., Elderfield, H., and Conte, M. H.: Calibration of Mg/Ca thermometry in planktonic foraminifera from a sediment trap time series, Paleooceanography, 18, 1–15, doi:10.1029/2002PA000846, 2003.

Badger, M. P. S., Lear, C. H., Pancost, R. D., Foster, G. L., Bailey, T. R., Leng, M. J., and Abels, H. A.: CO₂ drawdown following the middle Miocene expansion of the Antarctic Ice Sheet, Paleooceanography, 28, 42–53, doi:10.1002/palo.20015, 2013.

Bartoli, G., Hönisch, B., and Zeebe, R. E.: Atmospheric CO₂ decline during the Pliocene intensification of Northern Hemisphere glaciations, Paleooceanography, 26, PA4213, doi:10.1029/2010PA002055, 2011.

Beerling, D. J. and Royer, D. L.: Convergent Cenozoic CO₂ history, Nat. Geosci., 4, 418–420, 2011.

Berner, R. A. and Kothavala, Z.: GEOCARB III: A revised model of atmospheric CO₂ over Phanerozoic time, Am. J. Sci., 301, 182–204, 2001.

Brennan, S. T., Lowenstein, T. K., and Cendón, D. I.: The major ion composition of Cenozoic seawater: the past 36 million years from fluid inclusions in marine halite, Am. J. Sci., 313, 713–775, 2013.

Broecker, W. S. and Peng, T. H.: Tracers in the Sea, Lamont-Doherty Earth Observatory, Palisades, NY, 1982.

Burton, K. W. and Vigier, N.: Lithium isotopes as tracers in Marine and terrestrial environments, Handbook of Environmental Isotope Geochemistry, Springer, Berlin, Heidelberg, 41–59, 2012.

CARINA Group: Carbon in the Atlantic Ocean Region – the CARINA project: Results and Data, Version 1.0: Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, US Department of Energy, Oak Ridge, Tennessee, doi:10.3334/CDIAC/otg.CARINA.ATL.V1.0, 2009.

Catanzaro, E. J., Champion, C., Garner, E., Marinenko, G., Sappenfield, K., and W., S.: Boric Acid: Isotopic and Assay Standard Reference Materials NBS (US) Special Publications. National Bureau of Standards, Institute for Materials Research, Washington, DC, 1970.

Cramer, B., Miller, K., Barrett, P., and Wright, J.: Late Cretaceous–Neogene trends in deep ocean temperature and continental ice volume: Reconciling records of benthic foraminiferal geochemistry (δ¹⁸O and Mg/Ca) with sea level history, J. Geophys. Res.–Oceans, 116, C12023, doi:10.1029/2011JC007255, 2011.

Curry, W. B. and Oppo, D. W.: Glacial water mass geometry and the distribution of δ¹³C of ΣCO₂ in the Western Atlantic Ocean. Paleooceanography, 20, PA1017, doi:10.1029/2004PA001021, 2005,

Delaney, M. L., Be, A. W. H., and Boyle, E. A.: Li, Sr, Mg and Na in foraminiferal calcite shells from laboratory culture, sediment traps and sediment cores, Geochim. Cosmochim. Ac., 49, 1327–1341, 1985.

Deyhle, A. and Kopf, A.: Possible influence of clay contamination on B isotope geochemistry of carbonaceous samples, Appl. Geochem., 19, 737–745, 2004.

Edwards, N. R. and Marsh, R.: Uncertainties due to transport-parameter sensitivity in an efficient 3-D ocean-climate model, Clim. Dynam., 24, 415–433, doi:10.1007/s00382-004-0508-8, 2005.

Elderfield, H., Yu, J., Anand, P., Kiefer, T., and Nyland, B.: Calibrations for benthic foraminiferal Mg/Ca paleothermometry and the carbonate ion hypothesis, Earth Planet. Sc. Lett., 250, 633–649, 2006.

Evans, D. and Muller, W.: Deep time foraminifera Mg/Ca paleothermometry: Nonlinear correction for secular change in seawater Mg/Ca, Paleooceanography, 27, PA4205, doi:10.1029/2012PA002315, 2012.

Fantle, M. S. and Tipper, E. T.: Calcium isotopes in the global biogeochemical Ca cycle: Implications for development of a Ca isotope proxy, Earth-Sci. Rev., 129, 148–177, 2014.

Foster, G., Lear, C. H., and Rae, J. W. B.: The evolution of pCO₂, ice volume and climate during the middle Miocene, Earth Planet. Sc. Lett., 341–344, 243–254, 2012.

Foster, G., Hönisch, B., Paris, G., Dwyer, G., Rae, J., Elliott, T., Gaillardet, J., Hemming, N., Louvat, P., and Vengosh, A.: Interlaboratory comparison of boron isotope analyses of boric acid, seawater and marine CaCO₃ by MC-ICPMS and NTIMS, Chem. Geol., 358, 1–14, 2013.

Foster, G. L.: Seawater pH, pCO₂ and [CO₃²⁻] variations in the Caribbean Sea over the last 130 kyr: A boron isotope and B/Ca study of planktic foraminifera, Earth Planet. Sc. Lett., 271, 254–266, 2008.

Foster, G. L., Pogge von Strandmann, P. A. E., and Rae, J. W. B.: Boron and magnesium isotopic composition of seawater, Geochem. Geophy. Geosy., 11, Q08015, doi:10.1029/2010GC003201, 2010.

Froelich, F. and Misra, S.: Was the late Paleocene-early Eocene hot because Earth was flat? An ocean lithium isotope view of mountain building, continental weathering, carbon dioxide, and Earth's Cenozoic climate, Oceanography, 27, 36–49, 2014.

Galbraith, E. D., Kwon, E. Y., Bianchi, D., Hain, M. P., and Sarmiento, J. L.: The impact of atmospheric pCO₂ on carbon isotope ratios of the atmosphere and ocean, Global Biogeochem. Cy., 9, 307–324, doi:10.1002/2014GB004929, 2015.

Goodwin, P., and Lauderdale, J. M.: Carbonate ion concentrations, ocean carbon storage, and atmospheric CO₂, Global Biogeochem. Cy., 27, 882–893, doi:10.1002/gbc.20078, 2013.

Gradstein, F. M., Ogg, J. G., Schmitz, M., and Ogg, G.: The Geologic Time Scale 2012, Boston, Elsevier, 1144 pp., doi:10.1016/B978-0-444-59425-9.00004-4, 2012.

Greenop, R., Foster, G. L., Wilson, P. A., and Lear, C. H.: Middle Miocene climate instability associated with high-amplitude CO₂ variability, Paleooceanography, 29, 845–853, 2014.

Greenop, R., Hain, M. P., Sosdian, S. M., Oliver, K. I. C., Goodwin, P., Chalk, T. B., Lear, C. H., Wilson, P. A., and Foster, G. L.: Planktic and benthic boron and carbon isotopes from ODP sites 121–758, 165–999, 154–926 and 122–761, Dataset #871900, https://doi.pangaea.de/10.1594/PANGAEA.871900, 2017.

Griffith, E., Paytan, A., Caldeira, K., Bullen, T., and Thomas, E.: A Dynamic Marine Calcium Cycle During the Past 28 Million Years, Science, 322, 1671–1674, 2008.
Hain, M. P., Sigman, D. M., and Haug, G. H.: Carbon dioxide effects of Antarctic stratification, North Atlantic IntermEDIATE Water formation, and subantarctic nutrient drawdown during the last ice age: Diagnosis and synthesis in a geochemical box model, Global Biogeochem. Cy., 24, GB4023, doi:10.1029/2010GB003790, 2010.

Hain, M. P., Sigman, D. M., and Haug, G. H.: The Biological Pump in the Past, Treatise on Geochemistry, 2nd Edn., 8, 485–517, doi:10.1016/B978-0-08-095975-7.00618-5, 2014a.

Hain, M. P., Sigman, D. M., and Haug, G. H.: Distinct roles of the Southern Ocean and North Atlantic in the deglacial atmospheric radiocarbon decline, Earth Planet. Sc. Lett., 394, 198–208, doi:10.1016/j.epsl.2014.03.020, 2014b.

Hain, M. P., Sigman, D. M., Higgins, J. A., and Haug, G. H.: The effects of secular calcium and magnesium concentration changes on the thermodynamics of seawater acid/base chemistry: Implications for Eocene and Cretaceous ocean carbon chemistry and buffering, Global Biogeochem. Cy., 29, 517–533, doi:10.1002/2014GB004986, 2015.

Hasiuk, F. and Lohmann, K.: Application of calcite Mg partitioning functions to the reconstruction of paleocean Mg/Ca, Geochim. Cosmochim. Ac., 74, 6751–6763, 2010.

Hathorne, E. C. and James, R. H.: Temporal record of lithium in seawater: A tracer for silicate weathering?, Earth Planet. Sc. Lett., 246, 393–406, 2006.

Haug, G. H. and Tiedemann, R.: Effect of the formation of the Istmus of Panama on Atlantic Ocean thermohaline circulation, Nature, 393, 673–676, 1998.

Hemleben, Ch., Spindler, M., Breitinger, I., and Ott, R.: Morphological and physiological responses of *Globigerinoides sacculifer* (Brady) under varying laboratory conditions, Mar. Micropaleontol., 12, 305–324, 1987.

Hemming, N. G. and Hanson, G. N.: Boron isotopic composition and concentration in modern marine carbonates, Geochim. Cosmochim. Ac., 56, 537–543, 1992.

Henehan, M. J., Rae, J. W. B., Foster, G. L., Erez, J., Prentice, K. C., Kucera, M., Bostock, H. C., Martínez-Botí, M. A., Milton, J. A., Wilson, P. A., Marshall, B. J., and Elliott, T.: Calibration of the boron isotope proxy in the planktonic foraminifera *Globigerinoides ruber* for use in palaeo-CO₂ reconstruction, Earth Planet. Sc. Lett., 364, 111–122, 2013.

Hindshaw, R. S., Bourdon, B., Pogge von Strandmann, P. A. E., Vigier, N., and Burton, K. W.: The stable calcium isotopic composition of rivers draining basaltic catchments in Iceland, Earth Planet. Sc. Lett., 374, 173–184, 2013.

Hodell, D. A., Mueller, P. A., and Garrido, J. R.: Variations in the strontium isotopic composition of seawater during the Neogene, Geology, 11, 24–27, 1991.

Holbourn, A., Kuhnt, W., Simo, J., and Li, Q.: Middle Miocene isotope stratigraphy and paleoceanographic evolution of the northwest and southwest Australian margins (Wombat Plateau and Great Australian Bight), Palaeogeogr. Palaeocl., 208, 1–22, 2004.

Holden, P. B., Edwards, N. R., Müller, S. A., Oliver, K. I. C., Death, R. M., and Ridgwell, A.: Controls on the spatial distribution of oceanic δ¹³C_DIC, Biogeosciences, 10, 1815–1833, doi:10.5194/bg-10-1815-2013, 2013.

Hönisch, B., Hemming, N. G., Archer, D., Siddall, M., and McManus, J. F.: Atmospheric Carbon Dioxide Concentration Across the Mid-Pleistocene Transition, Science, 324, 1551–1554, 2009.

Horita, J., Zimmermann, H., and Holland, H. D.: Chemical evolution of seawater during the Phanerozoic: Implications from the record of marine evaporites, Geochim. Cosmochim. Ac., 66, 3733–3756, 2002.

Kaczmarek, K., Nehrk, G., Misra, S., Bijma, J., and Elderfield, H.: Investigating the effects of growth rate and temperature on the B/Ca ratio and δ¹¹B during inorganic calcite formation, Chem. Geol., 421, 81–92, 2016.

Kakihana, H., Kotaka, M., Satoh, S., Nomura, M., and Okamoto, M.: Fundamental studies on ion-exchange separation of boron isotopes, Bull. Chem. Soc. Jap., 50, 158–163, 1977.

Keeling, C. D.: The Suess effect: ¹³Carbon-¹⁴Carbon interrelations, Environ. Int., 2, 229–300, 1979.

Key, R. M., Kozyr, A., Sabine, C. L., Lee, K., Wanninkhof, R., Bullister, J. L., Feely, R. A., Millero, F. J., Mordy, C., and Peng, T. H.: A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP), Global Biogeochem. Cy., 18, GB4031, doi:10.1029/2004GB002247, 2004.

Kisakürek, B., James, R. H., and Harris, N. B. W.: Li and δ⁷Li in Himalayan rivers: Proxies for silicate weathering?, Earth Planet. Sc. Lett., 237, 387–401, 2005.

Klochko, K., Kaufman, A. J., Yao, W. S., Byrne, R. H., and Tossell, J. A.: Experimental measurement of boron isotope fractionation in seawater, Earth Planet. Sc. Lett., 248, 276–285, 2006.

Lear, C. H., Mawbey, E. M., and Rosenthal, Y.: Cenozoic ben-thic foraminiferal Mg/Ca and Li/Ca records: Toward unlock-ing temperatures and saturation states, Paleceanography, 25, PA4215, doi:10.1029/2009PA001880, 2010.

Lee, K., Kim, T. W., Byrne, R. H., Millero, F. J., Feely, R. A., and Liu, Y. M.: The universal ratio of boron to chlorinity for the North Pacific and North Atlantic oceans, Geochim. Cosmochim. Ac., 74, 1801–1811, 2010.

Lemarchand, D. and Gaillardet, J.: Transient features of the erosion of shales in the Mackenzie basin (Canada), evidences from boron isotopes, Earth Planet. Sc. Lett., 245, 174–189, 2006.

Lemarchand, D., Gaillardet, J., Lewin, E., and Allegre, C. J.: The influence of rivers on marine boron isotopes and implications for reconstructing past ocean pH, Nature, 408, 951–954, 2000.

Li, G.-J. and West, A. J.: Evolution of Cenozoic seawater lithium isotopes: coupling of global denudation regime and shifting seawater sinks, Earth Planet. Sc. Lett., 401, 284–293, 2014.

Liu, W. G., Xiao, Y. K., Peng, Z. C., An, Z. S., and He, X. X.: Boron concentration and isotopic composition of halite from experiments and salt lakes in the Qaidam Basin, Geochim. Cosmochim. Ac., 64, 2177–2183, 2000.

Liu, X.-M., Wanner, C., Rudnick, R. L., and McDonough, W. F.: Processes controlling δ⁷Li in rivers illuminated by study of streams and groundwaters draining basalts, Earth Planet. Sc. Lett., 409, 212–224, 2015.

Lynch-Steiglitz, J., Stocker, T. F., Broecker, W. S., and Fairbanks, R. G.: The influence of air-sea exchange on the isotopic composition of oceanic carbon: Observations and modeling, Global Biogeochem. Cy., 9, 653–665, 1995.

Martínez-Botí, M. A., Foster, G. L., Chalk, T. B., Rohling, E. J., Sexton, P. F., Lunt, D. J., Pancost, R. D., Badger, M. P. S., and
Schmidt, D. N.: Plio-Pleistocene climate sensitivity from on a new high-resolution CO₂ record, Nature, 518, 49–54, 2015a.

Martinez-Botí, M. A., Marino, G., Foster, G. L., Ziveri, P., Henehan, M. J., Rae, J. W. B., Mortyn, P. G., and Vance, D.: Boron isotope evidence for oceanic CO₂ leakage during the last deglaciation, Nature, 518, 219–222, 2015b.

McCorkle, D. C., Corliss, B. H., and Farnham, C. A.: Vertical distributions and stable isotopic compositions of live (stained) benthic foraminifera from the North Carolina and California continental margins, Deep-Sea Res. Pt. I, 44, 983–1024, 1997.

Millot, R., Vigier, N., and Gaillardet, J.: Behaviour of lithium and its isotopes during weathering in the Mackenzie Basin, Canada, Geochim. Cosmochim. Ac., 74, 3897–3912, 2010.

Misra, S. and Froelich, P.: Lithium Isotope History of Cenozoic Seawater: Changes in Silicate Weathering and Reverse Weathering, Science, 335, 818–823, 2012.

Ockert, C., Gussone, N., Kaufhold, S., and Teichert, B. M. A.: Isotope fractionation during Ca exchange on clay minerals in a marine environment, Geochim. Cosmochim. Ac., 112, 374–388, 2013.

Olsen, A. and Ninneman, U. S.: Large δ¹³C gradients in the preindustrial North Atlantic revealed, Science, 330, 658–659, 2010.

Pälike, H., Lyle, M. W., Nishi, H., Raffi, I., Ridgwell, A., Gamage, K., Klaus, A., Acton, G., Anderson, L., Backman, J., Baldauf, J., Beltran, C., Bohaty, S. M., Bown, P., Busch, W., Channell, J. E. T., Chun, C. O. J., Delaney, M., Dewangan, P., Dunkley Jones, T., Edgar, K. M., Evans, H., Fitch, P., Foster, G. L., Guscone, N., Hasegawa, H., Hathorne, E. C., Hayashi, H., Herrle, J. O., Holbourn, A., Hovan, S., Hyeong, K., Iijima, K., Ito, T., Kamikuri, S., Kimoto, K., Kuroda, J., Leon-Rodriguez, L., Malinverno, A., Moore, T. C., Murphy, B. H., Murphy, D. P., Nakamura, H., Ogane, K., Ohneiser, C., Richter, C., Robinson, R., Rohling, E. J., Romero, O., Sawada, K., Scher, H., Schneider, L., Sluijs, A., Takata, H., Tian, J., Tsujimoto, A., Wade, B. S., Westerhold, T., Wilkens, R., Williams, T., Wilson, P. A., Yamamoto, Y., Yamamoto, S., Yamazaki, T., and Zeebe, R. E.: A Cenozoic record of the equatorial Pacific carbonate compensation depth, Nature, 488, 609–614 2012.

Palmer, M. R., Pearson, P. N., and Cobb, S. J.: Reconstructing past ocean pH-depth profiles, Science, 282, 1468–1471, 1998.

Paris, G., Gaillardet, J., and Louvat, P.: Geological evolution of seawater boron isotopic composition recorded in evaporites, Geology, 38, 1035–1038, 2010.

Park, H. and Schlesinger, W. H.: Global biogeochemical cycle of boron, Global Biogeochem. Cy., 16, 1072, doi:10.1029/2001GB001766, 2002.

Pearson, P. N. and Palmer, M. R.: Middle Eocene seawater pH and atmospheric carbon dioxide concentrations, Science, 284, 1824–1826, 1999.

Pearson, P. N. and Palmer, M. R.: Atmospheric carbon dioxide concentrations over the past 60 million years, Nature, 406, 695–699, 2000.

Pearson, P. N. and Wade, B. S.: Taxonomy and Stable Isotope Palaeoecology of Well-Preserved Planktonic Foraminifera from the Uppermost Oligocene of Trinidad, J. Foramin. Res., 39, 191–217, 2009.

Pearson, P. N., Foster, G. L., and Wade, B. S.: Atmospheric carbon dioxide through the Eocene-Oligocene climate transition, Nature, 461, 1110–1113, 2009.

Pistiner, J. S. and Henderson, G. M.: Lithium-isotope fractionation during continental weathering processes, Earth Planet. Sc. Lett., 214, 327–339, 2003.

Pogge von Strandmann, P. A. E., Burton, K. W., James, R. H., van Calsteren, P., Gislason, S. R., and Sigfússon, B.: The influence of weathering processes on riverine magnesium isotopes in a basaltic terrain, Earth Planet. Sc. Lett., 276, 187–197, 2008.

Pogge von Strandmann, P. A. E., Forshaw, J., and Schmidt, D. N.: Modern and Cenozoic records of seawater magnesium from foraminiferal Mg isotopes, Biogeosciences, 11, 5155–5168, doi:10.5194/bg-11-5155-2014, 2014.

Pogge von Strandmann, P. A. E. and Henderson, G. M.: The Li isotope response to mountain uplift, Geology, 43, 67–70, doi:10.1130/G36162.1, 2014.

Rae, J. W. B., Foster, G. L., Schmidt, D. N., and Elliott, T.: Boron isotopes and B / Ca in benthic foraminifera: Proxies for the deep ocean carbonate system, Earth Planet. Sc. Lett., 302, 403–413, 2011.

Raitzsch, M. and Hönisch, B.: Cenozoic boron isotope variations in benthic foraminifers, Geology, 41, 591–594, 2013.

Ridgewell, A.: A mid Mesozoic revolution in the regulation of ocean chemistry, Mar. Geol., 217, 339–357, 2005.

Ridgwell, A., Hargreaves, J. C., Edwards, N. R., Annan, J. D., Lenton, T. M., Marsh, R., Yool, A., and Watson, A.: Marine geochemical data assimilation in an efficient Earth System Model of global biogeochemical cycling, Biogeosciences, 4, 87–104, doi:10.5194/bg-4-87-2007, 2007.

Rose, E. F., Chaussidon, M., and France-Lanord, C.: Fractionation of boron isotopes during erosion processes: the example of Hi-malayan rivers, Geochim. Cosmochim. Ac., 64, 397–408, 2000.

Rose-Koga, E. F., Sheppard, S. M. F., Chaussidon, M., and Carignan, J.: Boron isotopic composition of atmospheric precipitations and liquid-vapour fractionations, Geochim. Cosmochim. Ac., 70, 1603–1615 2006.

Sanyal, A., Hemming, N. G., Hanson, G. N., and Broecker, W. S.: Evidence for a higher pH in the glacial ocean from boron isotopes in foraminifera, Nature, 373, 243–236, 1995.

Sanyal, A., Bijma, J., Spero, H., and Lea, D. W.: Empirical relationship between pH and the boron isotopic composition of Globigerinoides sacculifer: Implications for the boron isotope paleo-pH proxy, Paleooceanography, 16, 515–519, 2001.

Schlitzer, R.: Ocean Data View, http://www.awi-bremerhaven.de/GEO/ODV, 2016.

Seki, O., Foster, G. L., Schmidt, D. N., Mackensen, A., Kawamura, K., and Pancost, R. D.: Alkenone and boron-based Pliocene $pCO_2$ records, Earth Planet. Sc. Lett., 292, 201–211, 2010.

Shipboard Scientific Party: Proc. ODP, Init. Repts., edited by: Peirce, J., Weissel, J., Taylor, E., Dehn, J., Driscoll, M., Farrell, J., Fourtanier, E., Frey, F., Garrison, P. D., Gee, J. S., Gibson, I.L., Janecek, T., Klootwijk, C., Lawrence, J. R., Litteke, R., Newman, J. S., Nomura, R., Owen, R. M., Pospichal, J. J., Rea, D. K., Resiwati, P., Saunders, A. D., Smit, J., Smith, G. M., Tamaki, K., Weis, D., and Wilkinson, K., College Station, TX (Ocean Drilling Program), 121, 359–453, doi:10.2973/odp.proc.ir.121.112.1989, 1989.

Shipboard Scientific Party: Proc. ODP, Init. Repts., edited by: Curry, W. B., Shackleton, N. J., Richter, C., Backman, J. E., Bassinot, F., Bickert, T., Chaisson, W. P., Cullen, J. L., deMenocal, P., Dobson, D. M., Ewert, L., Qrutzner, J., Hagelberg, T. K., Hampt, Q.,
Harris, S. E., Herbert, T. D., Moran, K., Murayama M., Murra, D. W., Pearson, P. N., Raffi, I., Schneider, D. A., Tiedemann, R., Valet, J.-P., Weedon, G. P., Yasuda, H., and Zachos, J. C., College Station, TX (Ocean Drilling Program), 154, 153–232, doi:10.2973/odp.proc.ir.154.105.1995, 1995.

Shipboard Scientific Party: Proc. ODP, Init. Repts., edited by: Sigurdsson, H., Leckie, R. M., Acton, G. D., Abrams, L. J., Bralower, T. J., Carey, S. N., Chaisson, W. P., Cotillon, P., Cuningham, A. D., D'Hondt, S. L., Droxler, A. W., Qalbrun, B., Gonzalez, J., Haug, G., Kameo, K., King, J., Lind, I. L., Louvel, V., Lyons, T. W., Murray, R. W., Mutti, M., Myers, G., Pearce, R. B., Pearson, D. G., Peterson, L. C., and Röhl, U., College Station, TX (Ocean Drilling Program), 165, 131–230, doi:10.2973/odp.proc.ir.165.104.1997, 1997.

Sigman, D. M., McCorkle, D. C., and Martin, W. R.: The calcite lysocline as a constraint on glacial/interglacial low-latitude production changes, Global Biogeochem. Cy., 12, 409–427, 1998.

Simon, L., Lecuyer, C., Marechal, C., and Coltice, N.: Modelling the geochemical cycle of boron: Implications for the long-term $\delta^{11}$B evolution of seawater and oceanic crust, Chem. Geol., 225, 61–76, 2006.

Smith, H. J., Spivack, A. J., Staudigel, H., and Hart, S. R.: The boron isotopic composition of altered oceanic crust, Chem. Geol., 126, 119–135, 1995.

Sosdian, S. M., Greenop, R., Hain, M. P., Foster, G. L., Pearson, P. N., and Lear, C. H.: Constraining the evolution of Neogene ocean carbonate chemistry using the boron isotope pH proxy, in preparation, 2017.

Spero, H., Mielke, K., Kalve, E., Lea, D., and Pak, D.: Multispecies approach to reconstructing eastern equatorial Pacific thermocline hydrography during the past 360 kyr, Paleooceanography, 18, 1022, doi:10.1029/2001GC000200, 2003.

Spezzaferri, S., Kucera, M., Pearson, P. N., Wade, B. S., Rappo, S., Poole, C. R., Morard, R., and Stalder, C.: Fossil and genetic evidence for the polyphyletic nature of the planktonic foraminifera "Globigerinoides", and description of the new Genus *Trilobatus*, PLoS ONE, 10, e0128108, doi:10.1371/journal.pone.0128108, 2015.

Spivack, A. J. and Edmond, J. M.: Boron isotope exchange between seawater and the oceanic crust, Geochim. Cosmochim. Ac., 51, 1033–1043, 1987.

Takahashi, T., Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D. C. E., Schuster, U., Metzl, N., Yoshikawa-Inoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Kortzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T. S., Tilbrook, B., Johannessen, T., Olsen, A., Bellerby, R., Wong, C. S., Delille, B., Bates, N. R., and de Baar, H. J. W.: Climatological mean and decadal change in surface ocean $pCO_2$, and net sea-air $CO_2$ flux over global oceans, Deep-Sea Res. Pt. II, 56, 554–557, 2009.

Tipper, E. T., Galy, A., and Bickle, M. J.: Riverine evidence for a fractionated reservoir of Ca and Mg on the continents: Implications for the oceanic Ca cycle, Earth Planet. Sc. Lett., 247, 267–279, 2006a.

Tipper, E. T., Galy, A., Gaillardet, J., Bickle, M. J., Elderfield, H., and Carder, E. A.: The magnesium isotope budget of the modern ocean: Constraints from riverine magnesium isotope ratios, Earth Planet. Sc. Lett., 250, 241–253, 2006b.

Tomascak, P. B.: Developments in the Understanding and Application of Lithium Isotopes in the Earth and Planetary Sciences, Rev. Mineral. Geochem., 55, 153–195, 2004.

Tyrrell, T. and Zeebe, R. E.: History of carbonate ion concentration over the last 100 million years, Geochim. Cosmochim. Ac., 68, 3521–3530, 2004.

Vengosh, A., Starinsky, A., Kolodny, Y., Chivas, A. R., and Raab, M.: Boron isotope variations during fractional evaporation of sea-water – new constraints on the marine vs. nonmarine debate, Geology, 20, 799–802, 1992.

Vigier, N. and Goddéris, Y.: A new approach for modeling Cenozoic oceanic lithium isotope paleo-variations: the key role of climate, Clim. Past, 11, 635–645, doi:10.5194/cp-11-635-2015, 2015.

Wanner, C., Sonnenthal, E. L., and Liu, X.-M.: Seawater $\delta^7$Li: a direct proxy for global $CO_2$ consumption by continental silicate weathering?, Chem. Geol., 381, 154–167, 2014.

Wimpenny, J., Colla, C. A., Yin, Q.-Z., Rustad, J. R., and Casey, W. H.: Investigating the behaviour of Mg isotopes during the formation of clay minerals, Geochim. Cosmochim. Ac., 128, 178–194, 2014.

Wombacher, F., Eisenhauer, A., Böhm, F., Gussone, N., Regenberg, M., Dullo, W. C., and Rüggeberg, A.: Magnesium stable isotope fractionation in marine biogenic calcite and aragonite, Geochim. Cosmochim. Ac., 75, 5797–5818, 2011.

You, C.F., Spivack, A. J., Smith, J. H., and Gieskes, J. M.: Mobilization of boron in convergent margins: Implications for the boron geochemical cycle, Geology, 21, 207–210, 1993.

Zeebe, R. E. and Wolf-Gladrow, D. A.: $CO_2$ in seawater, equilibrium, kinetics, isotopes IN Elsevier oceanography series, Amsterdam, PAYS-BAS, Elsevier, XIII, 346 pp., 2001.

Zeeden, C., Hilgen, F., Westerhold, T., Lourens, L., Röhl, U., and Bickert, T.: Revised Miocene splice, astronomical tuning and calcareous plankton biochronology of ODP Site 926 between 5 and 14.4 Ma, Palaeogeogr. Palaeocl., 369, 430–451, 2013.