Reconstruction of Cenozoic δ11Bsw Using a Gaussian Process

The boron isotope ratio of seawater (δ11Bsw) is a parameter which must be known to reconstruct palaeo pH and CO2 from boron isotope measurements of marine carbonates. Beyond a few million years ago, δ11Bsw is likely to have been different to modern. Palaeo δ11Bsw can be estimated by simultaneously constraining the vertical gradients in foraminiferal δ11B (Δδ11B) and pH (ΔpH). A number of subtly different techniques have been used to estimate ΔpH in the past, all broadly based on assumptions about vertical gradients in oxygen, and/or carbon, or other carbonate system constraints. In this work we pull together existing data from previous studies, alongside a constraint on the rate of change of δ11Bsw from modeling. We combine this information in an overarching statistical framework called a Gaussian Process. The Gaussian Process technique allows us to bring together data and constraints on the rate of change in δ11Bsw to generate random plausible evolutions of δ11Bsw. We reconstruct δ11Bsw, and by extension palaeo pH, across the last 65Myr using this novel methodology. Reconstructed δ11Bsw is compared to other seawater isotope ratios, namely Sr87/86 ${}^{87/86}\mathrm{S}\mathrm{r}$ , Os187/188 ${}^{187/188}\mathrm{O}\mathrm{s}$ , and δ7Li, which we also reconstruct with Gaussian Processes. Our method provides a template for incorporation of future δ11Bsw constraints, and a mechanism for propagation of uncertainty in δ11Bsw into future studies.


Introduction
Boron and its two isotopes ( 10 B and 11 B) exist in seawater with speciation governed by acid-base equilibrium.There are four degrees of freedom in the boron isotope system, which we typically express with the following five parameters: δ 11 B 4 (the isotopic composition of borate), δ 11 B sw (the isotope composition of seawater), ɛ (the fractionation factor), pK * B (the apparent equilibrium constant for boron in solution), and pH.Knowing any four of these parameters allows us to calculate the fifth (Zeebe & Wolf-Gladrow, 2001).For a typical paleoclimatological application, we wish to determine palaeo pH, which is done by measuring δ 11 B calcite and translating this to δ 11 B 4 (potentially requiring a species specific calibration).δ 11 B 4 is then combined with the fractionation factor (ɛknown from Klochko et al. (2006); Nir et al. (2015)), the apparent equilibrium constant for boron in seawater (pK * B -which is estimated from temperature, pressure, and seawater composition (Dickson & Goyet, 1994)), and the boron isotope ratio of seawater (δ 11 B sw ).We refer the interested reader to Marschall and Foster (2018) for a full description of boron systematics in seawater, but for the purposes of this work it is sufficient to say that one of the key parameters that must be established in order to reconstruct palaeo ocean pH is δ 11 B sw .Modern day seawater has a δ 11 B sw of 39.61 ± 0.04 ‰ (Foster et al., 2010) -however δ 11 B sw is very likely to have been different in the past (Lemarchand et al., 2000).
Based on assessments of the relevant input and output fluxes to the ocean, boron is thought to have a relatively long residence time in seawater of roughly 10 million years (Broecker & Peng, 1982)-that is, the average amount of time an atom of boron spends in the ocean.It is estimated from geochemical box modeling of the boron cycle that δ 11 B sw might have changed by up to 0.1 ‰/Myr (Lemarchand et al., 2000).This means that beyond a few million years ago, ocean δ 11 B sw might have been appreciably different to modern day seawater, which must be taken into account when calculating palaeo pH.Unfortunately we currently have no direct archives of δ 11 B sw .Halites show some promise as a direct proxy for δ 11 B sw (Paris et al., 2010), but uncertainty remains about the observed variability in halite δ 11 B reconstructions, so δ 11 B sw is currently calculated using more indirect techniques.
The primary technique for estimating, or placing limits on, palaeo δ 11 B sw is to exploit vertical gradients in the ocean (Anagnostou et al., 2016;Foster et al., 2012;Greenop et al., 2017;Henehan et al., 2019).This works due to the sigmoidal relationship between δ 11 B 4 and pH (as shown in Figure 1a).The non-linearity of this relationship makes it possible to calculate δ 11 B sw given two δ 11 B 4 datapoints and an estimate of the change in pH across this δ 11 B 4 gradient.The vertical gradient in pH (ΔpH) is correlated to the vertical gradients in carbon isotopes (Δδ 13 C), carbon concentration (ΔDIC), and oxygen concentration (Δ[O 2 ], linked with Apparent Oxygen Utilisation-AOU) due to their shared controlling processes (see Figure 2).Vertical gradients in all of these parameters can therefore be related to δ 11 B sw .Paired measurements which allow the vertical gradient in δ 11 B 4 to be established are currently only available at select times within the Cenozoic-shown below in Table 1.
The simplest option for reconstructing palaeo δ 11 B sw from measured δ 11 B in carbonates is to assume a constant ΔpH through time equal to modern, an approach used in Raitzsch and Hönisch (2013) (which assumed a ΔpH of 0.3 units).By combining this with the sigmoidal relationship of δ 11 B 4 and pH, it is possible to estimate δ 11 B sw (see Figure 1).To improve upon the assumption of a constant ΔpH through time, one can use foraminiferal Δδ 18 O and Δδ 13 C to guide identification of palaeo depth habitats of particular species, and to place reasonable limits on ΔpH (Anagnostou et al., 2016;Greenop et al., 2017;Palmer et al., 1998).For δ 13 C, this works because Δδ 13 C and ΔpH are jointly controlled by subsurface remineralization of organic carbon, which releases isotopically light carbon into the water and causes acidification (decreasing pH)-as illustrated in Figure 2. The rationale for such constraints is broadly justified by the difficulty of reversing the sign of the gradients (which would require a fundamental change in biogeochemical dynamics), demonstrated by exploring a wide swath of model parameter space (Greenop et al., 2017) in carbon cycle models such as CYCLOPS (Hain et al., 2010), LOSCAR (Zeebe, 2012), and cGENIE (Ridgwell et al., 2007).That said, it is possible that models may not yet be capable of representing the full plethora of possible carbonate chemistry states that the ocean can truly inhabit.Despite recent progress (Caves Rugenstein et al., 2019;Derry, 2022) our understanding of Cenozoic carbon cycle fluxes, and their drivers, remains incomplete, as does our understanding of their relationship to oceanic properties such as DIC and pH, and by extension the vertical gradients in these ocean properties.There are circumstances which would result in the "normal" relationships between these parameters breaking down, for instance if organic carbon is predominantly respired by sulfate reducers (which would alter local alkalinity and pH).While these complications are rare, it is prudent to keep in mind their potential to impact estimates of δ 11 B sw .Despite these challenges, using model predicted gradients in δ 13 C and pH with paired planktic-benthic (or surface-subsurface) foraminiferal δ 11 B data remains attractive, as it requires only stable carbon isotope data (analyses which are routinely performed), and a model capable of providing estimates of feasible ocean profiles of δ 13 C) and pH.
δ 11 B sw from paired surface-deep samples can also be refined by estimation of AOU (Anagnostou et al., 2016).Surface ocean water oxygen saturation is mainly controlled by temperature (which is already a requirement for the δ 11 B to pH calculation) and it is assumed based on modern observations that samples from deeper dwelling planktic foraminifera can not have inhabited anoxic waters (Hull et al., 2011).Similar to pH and δ 13 C, the vertical gradient in oxygen concentration is primarily a function of how much remineralization of organic carbon has occurred (see Figure 2).Remineralization adds CO 2 to the water, increasing the DIC concentration.A Redfield Ratio can be used to convert constraints on ΔO 2 to ΔDIC, however the impact on ΔpH is complicated by the buffering influence of alkalinity (Subhas et al., 2022).The calculation therefore requires that we know (or assume) a second carbonate system parameter (in addition to pH) to make ΔDIC calculable.This is not a major imposition, as quantification of palaeo CO 2 from boron isotope derived pH already requires assumption or knowledge of a second carbonate system parameter.Such an estimate may itself be derived from the same type of carbon cycle models used to estimate Δδ 13 C and ΔpH.Carbon cycle models have been used to provide an estimate of saturation state for specific sites or times in the past, or to provide estimates of alkalinity (Anagnostou et al., 2016;Henehan et al., 2019), which simultaneously informs δ 11 B sw and palaeo CO 2 reconstructed from boron isotope derived pH.
Overall, it is therefore necessary to balance a number of requirements when reconstructing δ 11 B sw .The established δ 11 B sw must produce a reasonable: ΔpH, Δδ 13 C, saturation state and AOU.Finding the space in which all these parameters are viable allows quantification of δ 11 B sw with attendant uncertainty.The boron isotope proxy is fortunate however, in that almost all data and assumptions required to estimate δ 11 B sw are already a part of the δ 11 B 4 to CO 2 calculation.Oxygen and carbon dioxide concentrations are set in surface waters by their atmospheric concentrations and Henry's Law (which is temperature dependent).Surface ocean carbon is taken up by biomass to form organic tissues, which are preferentially enriched in the lighter carbon isotope.When that biomass is exported to depth and remineralized, organic carbon is oxidized (consuming oxygen) and releasing CO 2 , which drives acidification.These processes cause correlated gradients in pH, aqueous CO 2 , aqueous O 2 , saturation state, and δ 13 C.As δ 11 B sw is a well mixed signal, differences in δ 11 B 4 between the surface and subsurface are primarily driven by the gradient in pH.δ 11 B sw is homogenous throughout the water column due to the long 10 Myr residence time of boron (Foster et al., 2010;Lemarchand et al., 2000).
Most previous applications of the techniques described above have focused on individual time slices of the Cenozoic (as described in Section 2 and Table 1).By combining each of these individual studies, we are able to reach a critical mass of information, allowing us to produce an estimate for how δ 11 B sw evolved throughout the Cenozoic.Our study represents an advance compared to two particularly relevant previous works that provided Cenozoic timescale estimates of δ 11 B sw .Raitzsch and Hönisch (2013) predicted Cenozoic δ 11 B sw by assuming a constant ΔpH and a linear trend in deep ocean pH, meaning reconstructions of pH using this δ 11 B sw are only able to reconstruct the initial assumption (J.W. B. Rae, 2018) (further impacts of which are discussed in Section 5.3).J. W. Rae et al. (2021) compiled and reanalyzed marine palaeo CO 2 proxy data (including boron isotopes), and therefore required a curve for δ 11 B sw .However, generation of a Cenozoic δ 11 B sw curve was not the primary focus of that study and the authors noted the need to improve interpolation and constrain uncertainty.Here we provide a more robust estimate of the evolution of δ 11 B sw through time by bringing existing and updated δ 11 B sw constraints together into a rigorous statistical framework called a Gaussian Process.
The Gaussian Process is a flexible Bayesian non-parametric regression technique.We provide further details in the methods, a schematic illustration in Figure 3, and a full description in Supporting Information S1.Intuitively, we can consider the Gaussian Process as a method for creating stochastic smooth curves influenced to go through datapoints.For δ 11 B sw , we provide data constraints and prescribe the smoothness based on mechanistic modeling, which allows us to can draw many potential evolutions of δ 11 B sw .We are able to impose additional constraints (such as on the rate of change) by filtering the potential evolutions of δ 11 B sw to find those which are compatible with all available information-thus providing an improved estimate of δ 11 B sw .The realizations can be summarized to provide an estimate of δ 11 B sw at any specific time or used as inputs to other models as a method for propagating uncertainty in δ 11 B sw .

Data
To create our curves, we combine observational data regarding levels of δ 11 B sw at specific time points from a range of previous studies.These observational estimates come with various forms of uncertainty, primarily due to how the raw (direct) measurements were processed/transformed into values for δ 11 B sw .Some of these previous studies report observed values that relate to central estimates of the underlying δ 11 B sw at the specific time of the sample (Gutjahr et al., 2017;Henehan et al., 2019Henehan et al., , 2020)).Others provide lower-and upper-bounds for the possible value of δ 11 B sw (Anagnostou et al., 2016).A further study provides complete probability distributions for possible values of δ 11 B sw given their observational constraints (Greenop et al., 2017).
Two of those studies reconstructed δ 11 B sw for a particular event (or short timeslice) (Gutjahr et al., 2017;Henehan et al., 2019), and two reconstructed δ 11 B sw over a wider time window, within which δ 11 B sw might have evolved (Anagnostou et al., 2016;Greenop et al., 2017).Data constraints on δ 11 B sw are varied.Some previous works presented central constraints of δ 11 B sw (Gutjahr et al., 2017;Henehan et al., 2019Henehan et al., , 2020)), while another provided lower or upper limits on δ 11 B sw (Anagnostou et al., 2016), and another provided full probability distributions for possible values of δ 11 B sw (Greenop et al., 2017).Two of those studies reconstructed δ 11 B sw for a particular event (or short timeslice) (Gutjahr et al., 2017;Henehan et al., 2019), and two reconstructed δ 11 B sw over a wider time window, within which δ 11 B sw might have evolved (Anagnostou et al., 2016;Greenop et al., 2017).Table 1 contains the compilation of datapoints used to create our estimate of δ 11 B sw using a Gaussian Process.Some of these data are as originally published, and others have been updated and refined.We use the δ 11 B sw values from Greenop et al. (2017) exactly as reported in the original study.The original published estimate from Henehan et al. (2019), however, is based on carbonate system calculations whose equilibrium constants use fitting parameters from the supplementary tables of Hain et al. (2015).These tables have since been found to contain inaccuracies (CenCO2PIP Consortium, 2023), and so here we updated carbonate system calculations, and hence δ 11 B sw constraints, based on corrected equilibrium constants packaged with Raitzsch et al. (2022) 2016) are integrated as full probability distributions instead of the lower/upper limits that were presented in the original study (see Table 1 and Data Set S1 in Supporting Information S2).
To supplement these existing δ 11 B sw constraints, we take planktic δ 11 B 4 data as presented in J. W.  2018), to calculate a range of valid δ 11 B sw by exploiting the sigmoidal shape of the relationship between pH and δ 11 B 4 (as seen in Figure 1).The maximum offset between δ 11 B 4 and δ 11 B sw is seen at low pH, and is described by the fractionation factor ɛ (or α).For any given estimate of δ 11 B 4 , the minimum valid δ 11 B sw is equal to δ 11 B 4 , and the maximum δ 11 B sw can be calculated by combination with the fractionation factor.The nature of the relationship between δ 11 B 4 and δ 11 B sw is examined in Figure 1a.For any given δ 11 B 4 , in order to sit somewhere on the sigmoidal curve, δ 11 B sw must be equal to or greater than δ 11 B 4 , but no more than δ 11 B 4 and the fractionation factor (the height of the sigmoid shape-strictly given by δ 11 B sw = δ 11 B 4 × α + ɛ).To create these additional constraints, we first bin the plantkic δ 11 B 4 data into 1 Myr intervals, then take the maximum measured δ 11 B 4 as the lower limit for δ 11 B sw , and use the minimum measured δ 11 B 4 to calculate the maximum δ 11 B sw .Uncertainty in measured δ 11 B foram is propagated to uncertainty in the possible δ 11 B sw , and the upper 99% quantile is used to estimate the maximum permissible δ 11 B sw .As described, lower limits can be calculated using this method, however these values are so low as to be uninformative across the Cenozoic.Upper limits created using this method are shown in Figure 4 (the horizontal pink bars).
In this work, we therefore have four potential categories of constraint: constraints where the uncertainty is well represented by a Gaussian distribution (Gaussian constraints), constraints where the uncertainty is not well represented by a Gaussian distribution (non Gaussian constraints), lower and upper limits, and limitations on the rate of change in δ 11 B sw from modeling.

Methods
Having assembled the data described in Section 2, we reconstruct δ 11 B sw using a Gaussian Process (GP).The Gaussian Process is a statistical technique that allows us to generate smooth time series conditioned to match data constraints (see Figure 3, or for a fuller, more rigorous description of a Gaussian Process-see Supporting Information S1).In the case of δ 11 B sw , we have an expectation of smoothness from the modeling work of Lemarchand et al. (2000), which we combine with data constraints as described in Section 2 and shown in Figure 4.A Gaussian Process works by using a kernel function, which encodes structure into the reconstruction, and hyperparameters which tune the behavior.Here we use the squared exponential kernel (also known as the Radial Basis Function), which expresses that nearby points are more likely to be similar to each other than distant points.There are two controlling hyperparameters-one expresses the length scale over which there is significant covariance, which is related to the expected rate of change in a signal.The other parameter is the noise scalewhich manifests as how much uncertainty is expected at times where no data constraints are available.The Gaussian Process can be used to generate random smooth lines with the prescribed characteristics even in the absence of data constraints (Figure 3a.), however data can be incorporated by adapting the covariance structure such that the statistical samples it generates will go through datapoints where permitted by the chosen hyperparameters (Figure 3b.).Uncertainty in data constraints can be incorporated directly into the Gaussian Process if uncertainty in the estimates is Gaussian in nature, and other types of constraint can be incorporated by adapting the approach.
The Gaussian Process is able to provide a number of equally likely, independent, stochastic time series, which are useful in reconstructing time series of palaeo data with an estimate of the uncertainty.Each obeys the smoothness constraint encapsulated by the kernel, while simultaneously attempting to match any available data.The result is that where a data constraint with low uncertainty is available, time series will be strongly influenced to go through this datapoint.Where data are sparse, or data constraints are uncertain, statistical replicates diverge to represent increasing uncertainty (see Figure 3b or Supporting Information S1 for a demonstration of this behavior).Because each sample drawn from the Gaussian process is independent and equally likely, we can straightforwardly apply filters to reject any sample that has undesirable properties (Figure 3c).As mentioned above, data constraints where uncertainties can be well represented by a Gaussian distribution can be directly assimilated into the Gaussian Process.Here we also wish to enforce three other types of constraint, lower/upper limitations on δ 11 B sw , limitations on the rate of change in δ 11 B sw , and data constraints with a non Gaussian uncertainty structure.Limitations (either on the value of δ 11 B sw or the rate of change in δ 11 B sw ) are relatively straightforward to enforce -we can compare each Gaussian Process sample to the limitation and simply reject those which are outside the established limits.Non-Gaussian constraints are slightly more difficult to integrate, but can also be done within a rejection framework by approximating each non-Gaussian constraint with a Gaussian distribution.Here we do this by creating Gaussian constraints with the same mean as their non-Gaussian counterpart, but with higher standard deviation.Samples can be then be strategically rejected to tweak the Gaussian Process such that it has effectively sampled from non-Gaussian constraints (the mechanics of this type of rejection sampling are clarified in Supporting Information S1).
One difficulty of rejection sampling, however, is dealing with simultaneous but mutually exclusive data constraints.Gaussian distributions are nonzero across their domain, meaning no value (however unlikely) is truly

Paleoceanography and Paleoclimatology
10.1029/2023PA004769 impossible.The nature of the non-Gaussian constraints in this work though, suggests that at any given time some δ 11 B sw values are impossible to reconcile with the data estimates.Given multiple constraints of this type, in close proximity to one another, it becomes exceedingly unlikely to draw samples from the Gaussian Process which have the requisite smoothness and pass through the possible region of every datapoint.We overcome this difficulty by giving each non-Gaussian constraint the possibility of being an outlier value.In summary, Gaussian Process samples are not always rejected completely where in disagreement with data, but the probability of acceptance is highest where in agreement with the data constraints, and lower (but non zero) further from the data constraints.
We reconstruct δ 11 B sw by conditioning a Gaussian Process on the few available central estimates of δ 11 B sw (shown in Table 1).Additional data constraints are integrated by filtering the generated statistical samples using rejection sampling (as described above and in Supporting Information S1).Any curves which fall outside the lower or upper limits are rejected, and some curves are rejected to adjust the Gaussian Process posterior such that it appropriately incorporates non Gaussian constraints.We use hyperparameters of 2 ‰ noise scale, and 10 Myr length scale.The Gaussian Process methodology does not enforce a specific rate of change, but the chosen noise scale and length scale enable an approximation of the rate of change, in this case roughly equivalent to a 0.2 ‰/Myr rate of change.This is on the faster end of agreement with the residence time and maximum rate of change from box modeling of the geochemical cycle of boron (Lemarchand et al., 2000).Lemarchand et al. (2000) suggest an upper limit to the rate of change of δ 11 B sw of 0.1 ‰/Myr.This estimate is based on modern day fluxes, so we allow greater uncertainty in the past when the boron fluxes in and out of the ocean may have been different, increasing the permissible maximum rate of change linearly to 0.7 ‰/Myr at 70 Ma (discussed further in Section 5.2).Output samples are filtered such that those with a rate of change greater than the values described above are rejected.This, alongside the rejection strategy used to integrate non-Gaussian constraints, leaves a subset of the originally generated samples.Approximately 10 in each 10,000 are accepted, and we run the algorithm until 10,000 samples have been accepted.

Results
Our reconstruction of δ 11 B sw suggests rather muted change across the Cenozoic relative to the change in δ 11 B 4 .δ 11 B 4 has changed by roughly 8 ‰, and we find that δ 11 B sw has been responsible for perhaps 2 ‰ of that change, leaving a 6 ‰ change driven predominantly by pH.The pattern of change in δ 11 B sw is non-linear, with a decline of roughly 1 ‰ from the early to mid Cenozoic, followed by an increase of approximately 1.5 ‰ from 30 Ma to 10 Ma, and a relatively stable value between those time periods and from 10 Ma until the present day.At times there is large uncertainty in the absolute value of δ 11 B sw , in particular around 30 Ma, however the relative change in δ 11 B sw is more certain and the rise in δ 11 B sw between 40 and 10 Ma-also illustrated in Figure 4-appears robust.Our reconstruction generally matches the evolution curve used by J. W. Rae et al. (2021), except in the interval between 20 Ma and 10 Ma.Within this interval both curves are guided predominantly by the data from Greenop et al. (2017), but J. W. Rae et al. (2021) used Greenop et al. (2017)'s binned values whereas our approach uses the probability distribution of each individual estimate.Our Gaussian Process methodology favors the higher δ 11 B values in this interval due to their lower reported uncertainty (relative to the lower δ 11 B sw estimates from Greenop et al. (2017)).This results in an overall higher estimate of δ 11 B sw during the Miocene, and implies pH was potentially lower than previously calculated at this time.However, we note the uncertainties in both our reconstruction and the individual δ 11 B sw constraints are high in the Neogene, and that some of the data constraints suggest rates of change in δ 11 B sw exceeding those compatible with geochemical box model predictions (as discussed further in Section 5.2).Uncertainty on δ 11 B sw is also particularly high in data gaps where δ 11 B 4 has yet to be measured.In particular we note a large window during the Oligocene within which no boron isotope data has yet been published, and which correspondingly has high uncertainty in δ 11 B sw .
We provide illustrative curves for δ 11 B 4 and pH (see Figure 4) by using a Gaussian Process to interpolate δ 11 B 4 , then combining this with our predicted δ 11 B sw and other ancillary parameters as in J. W. Rae et al. (2021).The Gaussian Process used to reconstruct δ 11 B 4 uses a length scale of 2 Myr to target relatively long term changes in pH, rather than individual paleoclimatic events.
The pattern of change we find in δ 11 B sw is broadly similar to temporal trends in a number of other related biogeochemical signals, such as oceanic 87/86 Sr, 187/188 Os, and δ 7 Li.Here we take published 87/86 Sr and δ 7 Li from Misra and Froelich (2012) Robinson et al. (2009).For each signal, we assume that data constraints have Gaussian uncertainty (with magnitude given with original estimates, except in the case of 87/86 Sr where an illustrative uncertainty is used), which allows us to perform a straightforward Gaussian Process reconstruction, with the residence time of each element informing the length scale used ( 87/86 Sr-5.1 Myr (Broecker & Peng, 1982), δ 7 Li-2.8Myr (Stoffynegli & Mackenzie, 1984), 187/188 Os -1 Myr).
Interpolated data products for both 87/86 Sr and δ 7 Li have been previously published (for instance in Misra and Froelich (2012)), but to our knowledge this is the first interpolated data product available for 187/188 Os (though an interpolation is shown in Torres et al. (2014).There remains uncertainty in the residence time of osmium, but it is thought to be extremely short relative to the other signals shown here (all estimates are <100kyr-see e.g.Oxburgh (2001)).We use a length scale of 1 Myr in the reconstruction to balance the available data density against the short residence time.Overall, our reconstructions use a more sophisticated fitting strategy than previous incarnations, which is guided primarily by the data and integrates information on residence time of each element.This allows us to provide a robust quantification of the uncertainty in each signal, within the limitations of currently available data density.
Each of seawater isotope signals we examine, like δ 11 B sw , shows a long term increase over the Cenozoic.Oceanic strontium isotopes appear to show most similarity with our reconstructed δ 11 B sw , while we note resemblance between the evolutions of pH, lithium isotopes, and osmium isotopes (as shown in Figure 5).We are able to provide an extremely narrow uncertainty window on our 87/86 Sr reconstruction due to strontium's long residence time and relatively high data density (and data quality) over the Cenozoic.By contrast, 187/188 Os which has greater data density and comparable scale of uncertainty in datapoints, has more uncertainty in our reconstruction.This is due to osmium's extremely short residence time, and persists despite us using an artificially long length scale to reconstruct this signal.δ 7 Li by contrast has quite large uncertainties, due in large part to the uncertainty in the individual data constraints.We are unable to reconstruct the rapid shift in δ 7 Li at the K-Pg boundary, which occurs faster than the modern day residence time of lithium (see Section 5.2 for further discussion).Despite these challenges we believe reconstruction of these signals benefits from the Gaussian Process approach, in that uncertainties on our reconstruction are more representative, and the timescale of change in each signal aligns more closely to our expectations based on their residence times.That is roughly equivalent to an increase in pH of 0.4 units across the Cenozoic, or a decrease in hydrogen ion availability by 60%.By comparison, anthropogenic CO 2 release has driven a surface ocean pH change of roughly 0.15 units, equivalent to a 41% increase in hydrogen ion availability (Findlay et al., 2022) relative to preindustrial conditions.Relative to their respective initial conditions, the anthropogenic surface ocean [H + ] perturbation is therefore approximately two thirds the magnitude of long term change seen within the last 65 million years.

δ 11 B sw Rate of Change
A key feature of the Gaussian Process methodology used here is the ability to draw smooth potential evolutions of δ 11 B sw , which we augment by filtering samples to reject those with an unfeasibly high temporal gradient.The feasibility of the rate of change in δ 11 B sw (and the hyperparameters which tune the Gaussian Process smoothness) are primarily based on the geochemical box modeling work of Lemarchand et al. (2000).Lemarchand et al. (2000) suggest that the likely maximum rate of change in δ 11 B sw is 0.1 ‰/Myr.We acknowledge that this is based on modern fluxes-fluxes which are not exceptionally well constrained (Park & Schlesinger, 2002)-so (as described in Section 3) we allow increasing maximum rates of change further back into the past -up to 0.7 ‰/Myr at 70 Ma.We choose the value of 0.7 ‰/Myr at 70 Ma because it allows us to enforce a rate limit near 0.1 ‰/Myr when we are most sure of that value (close to modern), without enforcing a strong limit earlier in the Cenozoic when we are much less certain about the residence time of boron and potential rate of change in δ 11 B sw .0.7 ‰/Myr is large enough that it does not result in rejection of potential evolutions in the earlier half of the Cenozoic (see Figure S3 in Supporting Information S1), meaning the influence of the rate of change limit is mostly constrained to the rejection of samples as a result of their temporal gradient during the Neogene.
Suggestions of a higher rate of change in δ 11 B sw than previously recognized are present in the data set of Greenop et al. (2017).The surface-deep δ 11 B 4 pairs of Greenop et al. (2017) (as described in Section 2) appear to record a  Misra and Froelich (2012).Data constraining 87/86 Sr and δ 7 Li is sourced from Misra and Froelich (2012), while data constraining 187/188 Os is compiled from a number of sources (listed in the main text).We interpolate each of these signals with a Gaussian Process with hyperparameters guided by the residence time of each signal.bimodal distribution, with some indicating δ 11 B sw similar to modern, and others (mostly in the middle Miocene climatic optimum) suggesting δ 11 B sw roughly 2 ‰ lower.These would appear to suggest oscillations in δ 11 B sw more rapid than would be consistent with a rate of change of 0.1 ‰/Myr.We find that in order to match these data constraints, our Gaussian Process would need to have a length scale of approximately 3 Myr (which is very close to Park and Schlesinger (2002)'s estimate of the residence time of boron at 3.3 Myr).Our understanding of the primary difference between the estimate of Park and Schlesinger (2002) and Lemarchand et al. (2000)'s is that is driven by different estimate of the atmospheric fluxes of boron, which are particularly difficult to constrain.We draw evolutions consistent with the rate of change suggested by Lemarchand et al. (2000), which is the current paradigm for palaeo pH estimation from boron isotopes-however we note the possibility for future changes in our understanding of the rate at which δ 11 B sw may have changed.Our method requires that we give each non-Gaussian constraint (shown in Table 1) the possibility of being an outlier.The Gaussian Process then forges a path roughly through the center of these constraints, drawn slightly high by the purportedly lower uncertainties in the higher δ 11 B sw estimates.

Paleoceanography and Paleoclimatology
In summary, we proceed here under the existing paradigm of a long residence time for boron in seawater, and a slow rate of change in δ 11 B sw .However, we acknowledge the possibility that the rate of change in δ 11 B sw is faster than currently appreciated, and look to further modeling efforts and an improved understanding of the boron cycle to guide future interpretation of the data collated here.Faster rates of change may be explored with the algorithm presented here, but this increases the requirement for higher resolution data constraining δ 11 B sw .Indeed our work underscores the value of, and need for more, data-derived δ 11 B sw constraints outside of individual paleoclimatic events. 87/86 Sr, 187/188 Os, and δ 7 Li sw Figure 5 shows that there is a broad scale similarity between the temporal evolution of δ 11 B, δ 7 Li, 87/86 Sr and 187/188 Os.In particular though, we observe greater similarity between δ 11 B sw and 87/86 Sr, and also between pH, δ 7 Li and 187/188 Os.The similarity between 87/86 Sr and δ 11 B sw is present in the overall shape of both signals, and in both showing an inflection point at approximately the same time, around 35 Ma.Given 87/86 Sr and δ 11 B sw share many controls (such as weathering and hydrothermal influx), it is not surprising that their signals share some similarities.However we note that 87/86 Sr has been increasing continuously since 40 Ma, while our reconstruction of δ 11 B sw begins to decline at approximately 10 Ma.

Comparison to
We note an striking correspondence between δ 11 B 4 and δ 7 Li (shown in Figures 4 and 5, or together in Figure S1 in Supporting Information S1), in agreement with previous findings (Greenop et al., 2017;Raitzsch & Hönisch, 2013) as well as similarity between δ 11 B 4 and 187/188 Os.Previously the correspondence between δ 11 B 4 and δ 7 Li has been attributed to the overlap in drivers of δ 11 B sw and δ 7 Li.Under the two component model described in Section 5.1, we propose that similarity of δ 11 B 4 and δ 7 Li could either be due to a relationship between δ 7 Li and δ 11 B sw (as previously suggested by Raitzsch and Hönisch (2013)), or due to a relationship between δ 7 Li and pH, or a combination of both.δ 11 B sw and δ 7 Li share drivers, as do pH and δ 7 Li, which are linked to weathering, clay formation, and seafloor spreading.These confounding factors make it difficult to ascribe this correlation to either the δ 11 B sw or pH signals unambiguously (see Figure S1 in Supporting Information S1).If pH is prescribed to change slowly and linearly (as in Raitzsch and Hönisch (2013)), then, by necessity, δ 11 B sw will reflect δ 11 B 4requiring relatively rapid fluctuations in δ 11 B sw .However, as discussed in Section 5.2, we proceed here under the assumption that the rate of change in δ 11 B sw is as calculated in Lemarchand et al. (2000), removing the possibility of fast changes in δ 11 B sw , and necessitating that δ 11 B 4 and pH are tightly correlated.Therefore here, it is pH and δ 7 Li which have similar trajectories across the Cenozoic.
Notwithstanding the broad scale correlation between δ 7 Li and δ 11 B 4 , we note an interesting divergence between these signals at the K-Pg boundary.At this time, there is a large excursion in δ 7 Li of approximately 5 ‰ in scale (almost as large as all the change which occurs during the rest of the Cenozoic).This excursion occurs at the same time as a large excursion in 187/188 Os.While there is also a large perturbation in foraminiferal boron isotopes at approximately this time (Henehan et al., 2019), the temporal agreement is poor, and the nature of the two signals is different.δ 11 B 4 undergoes an excursion but rapidly recovers to near pre-perturbation levels, whereas δ 7 Li values stay low for the next ∼15 million years.This indicates that the drivers of pH and δ 7 Li can be decoupled.However, the recorded change in δ 7 Li at the K-Pg is extremely rapid, faster even than the modern residence time of lithium in seawater (approximately 1.2 Myr (Misra & Froelich, 2012)).It seems equally plausible that there is alternative, non-seawater, driver of foraminiferal δ 7 Li at this time.If this alternative control is linked to the carbonate system (as suggested by Vigier et al. (2015); Roberts et al. (2018)), similar effects could conceivably be influencing the correlation between δ 11 B 4 and δ 7 Li more broadly during the Cenozoic.At present, however, studies disagree as to the nature of carbonate system control on foraminiferal δ 7 Li (Roberts et al., 2018;Vigier et al., 2015).

Palaeo pH and CO 2
pH is linked to the ocean carbonate system, with a particularly close relationship to atmospheric CO 2 concentrations (Hain et al., 2018).The maximum offset between central estimates of the long term trajectory of the Cenozoic pH reconstruction of J. W. Rae et al. (2021) and our pH (shown in Figure 4) occurs during the Miocene (at approximately 12 Ma) and is 0.15 units in scale.All other factors being equal, a fall in pH of 0.15 units would suggest an increase in atmospheric CO 2 concentration of roughly 50%.Naïvely scaling Miocene CO 2 estimates from J. W. Rae et al. (2021) results in CO 2 concentrations of approximately 750 ppm after the Miocene Climatic Optimum (MCO).During the MCO, the offset between our predicted pH and previous work is slightly smaller, which would suggest (again, assuming all other factors are equal) an increase in atmospheric CO 2 of 45% relative to J. W. Rae et al. (2021) -suggesting an MCO CO 2 of approximately 860ppm.These suggestions of increased CO 2 relative to previous work come with the important caveat that, as mentioned above, the constraints on δ 11 B sw in this interval are highly variable, perhaps suggesting additional complicating factors which make δ 11 B sw difficult to estimate at this time.We note that some recent paleotemperature and modeling studies (e.g., Steinthorsdottir et al. (2021) and Herbert et al. (2022)) have made the suggestion that Miocene CO 2 might be higher than estimated in earlier reconstructions, with Burls et al. (2021) finding aspects of improved model-data agreement at CO 2 of around 850 ppm.Our reconstruction of δ 11 B sw (and consequently pH) has high overall uncertainties during the Neogene to reflect this.We encourage future work to provide additional constraints in this interval and across the Cenozoic.
Second, uncertainties remain in estimates of the second carbonate system parameter in the Miocene and throughout the Cenozoic.The difference in estimated CO 2 from this work and J. W. Rae et al. (2021) may be partially or wholly ameliorated by changing our expectations of the second carbonate system parameter.In this case, the suggested higher atmospheric CO 2 could be averted by a reduction in DIC relative to previous estimates -or the reality could lie somewhere in between, with atmospheric CO 2 concentration mildly elevated compared to previous reconstructions, and DIC mildly lowered.Understanding of possible ocean DIC at this time is mostly derived from carbon cycle box models, supplemented by suggestions from the B/Ca proxy (Sosdian et al., 2018).The range in DIC estimates is from roughly 1,200 to 2500 μmol/kg, meaning we are unable to disambiguate whether this record is truly indicative of higher CO 2 concentrations at this time, or lower DIC concentrations, or a combination of both.

Benefits of the Gaussian Process Approach
The Gaussian Process methodology allows us to integrate data constraints with limitations on the rate of change in δ 11 B sw from modeling, and (as described in Section 3) by tweaking the standard approach we are able to incorporate constraints with non-Gaussian uncertainty structures.Uncertainty in the reconstruction itself behaves intuitively, such that the spread in the reconstruction is guided by uncertainties in the data where available, and grows larger with increasing separation from data constraints.The shape of the reconstruction is not specified a priori, as would be the case with parametric fits.Instead, the reconstruction can take on almost any shape as guided by the data constraints and chosen hyperparameters, allowing us to model arbitrary shapes in the evolution of δ 11 B sw -and the other isotope systems reconstructed here (Figure 5).
While the standard approach to fitting a Gaussian process (constrained by data with Gaussian uncertainties) can directly predict the mean and variance of the signal being reconstructed, our approach is reliant on drawing samples from the Gaussian Process and then filtering them to adjust the posterior prediction.Once this process is complete we have 10,000 possible time series which are plausible evolutions in δ 11 B sw over the last 65 Myr.The 10,000 possible evolutions can be summarized by their mean, median, and/or 95% confidence interval, but retaining each of the possibilities makes it possible to propagate uncertainties into future data products.For instance, uncertainty in palaeo pH can be propagated by using a Monte Carlo approach whereby each of these evolutions is sampled alongside other required parameters to provide 10,000 possible evolutions of pH.Keeping each possible evolution maintains the embedded covariance structure, meaning that it is possible to calculate derivative properties (such as the change in δ 11 B sw or pH between one time and another), which would not be possible to do from metrics such as the mean and standard deviation.This is particularly beneficial in the context of reconstructing the palaeo carbonate system, as trends in parameters are often more robust or important than their absolute value.For instance, for short time windows (relative to the residence time of boron), although we may not know the absolute value of δ 11 B sw , we believe that it can not have changed substantially.Combining the δ 11 B sw statistical samples generated here with a Monte Carlo approach allows uncertainty to be propagated in such a way as to explore the full range of absolute values for δ 11 B sw while each sample preserves a reasonable Δδ 11 B sw (see Tierney et al. (2022) for an example of how an analogous approach was used to constrain change in atmospheric CO 2 concentration).As the 10,000 possible evolutions of δ 11 B sw that we provide are considered equally likely, it also allows further filtering.For instance if looking at a particular time period, or when new information comes to light, which means we are able to be more certain about the rate of change in δ 11 B sw , then the samples provided here can be refiltered to enforce the more restrictive condition-though naturally this will result in a decreased number of valid samples and weakened statistical power.
In summary, we believe the Gaussian Process approach provides benefits both in being able to provide a holistic representation of our understanding of δ 11 B sw from data and modeling, and also in terms of producing results which facilitate more sophisticated forms of onward uncertainty propagation.

Limitations
The Gaussian Process methodology used in this work has many properties that make it ideal for geochemical data interpolation, but also has a few caveats.The first is that, in the standard approach, input data constraints are expected to be Gaussian distributions.As discussed above and in Supporting Information S1, we adjust the standard Gaussian Process fitting methodology by using a rejection sampling strategy to incorporate other types of constraint.
The second limitation of the Gaussian Process approach is more fundamental.The two hyperparameters which tune the fit describe the length scale and noise scale, as described above (Section 3).Here, we use the residence time of each element as the length scale of the Gaussian Process, however those concepts are not identical.In particular, as most signals we are reconstructing here are isotope ratios, the concept of elemental residence time is not necessarily directly applicable.Nonetheless, we believe that using the residence time as a guide for the rate of change in these signals is an improvement over using parametric methods, or non-parametric methods with arbitrary smoothing parameters.In particular, for δ 11 B sw , we are not strongly reliant on the assumption that the Gaussian Process length scale is equivalent to the residence time because Lemarchand et al. (2000) provide a direct rate of change estimate which we use as a constraint.Given a residence time for boron of 10 Myr, and a rate of change of δ 11 B sw of 0.1 ‰/Myr, this implies a noise scale of 1 ‰.However, 1 ‰ is a small range for uncertainty where there are no data constraints.Increasing the noise scale results in rates of change incompatible with rate of change presented in Lemarchand et al. (2000) unless the length scale is commensurately increasedhowever this is then inconsistent with our understanding of the residence time of boron in seawater.Our solution to this is to use values which permit slightly faster changes than suggested by Lemarchand et al. (2000)-a length scale of 10 Myr, but a noise scale of 2 ‰, then filtering out results which are incompatibly fast.
Using the Gaussian Process with the methodology described here is highly computationally intensive.Each generated sample has only a small chance of being accepted, and it is necessary to try tens of millions of possibilities to achieve 10,000 viable statistical samples.It is inherent in the rejection sampling methodology to be inefficient in this way, and the more criteria that are used (or the more restrictive those criteria are) the less efficient this method becomes.For this use case, our approach takes approximately 40 hr of computational time on a standard desktop machine to generate the requisite 10,000 samples.In future, we may look to alternative statistical techniques to increase our efficiency and allow us to explore a greater range of possibilities, in particular with respect to gradient limitations and signal smoothness.
Our reconstruction of δ 11 B sw is also limited by our understanding of past ocean conditions.As discussed in Section 1, most constraints on δ 11 B sw are at some level dependent on models.Typically, a wide range of model conditions were used to predict vertical gradients in δ 13 C and pH, and an even wider range was used in uncertainty propagation to calculate these estimates (see for instance Greenop et al. (2017)).Nonetheless, it would be remiss not to acknowledge that models are a simplification of reality, meaning it is possible that the ocean occupied a Paleoceanography and Paleoclimatology 10.1029/2023PA004769 different mode in the past where the vertical gradients in δ 13 C and pH were decoupled or otherwise difficult to predict.Overall, the data constraints in this work are contingent upon carbon cycle model simulations producing realistic ranges for the δ 13 C versus pH gradient, and the validity of assumptions which determine AOU, while our reconstruction itself is dependent on these data and the rate of change determined from modeling of the boron cycle.

Conclusions
We provide Cenozoic reconstructions of δ 11 B sw by collating existing data constraints and integrating these into a Gaussian Process based statistical approach.This allows us to bring together varying types of constraint (including central estimates, lower and upper limits, and other forms of distribution) while rigorously propagating uncertainties.Our results suggest that δ 11 B sw was slightly higher than previously thought during the Miocene, but are generally in agreement with previous estimates of δ 11 B sw during the remainder of the Cenozoic.Generally speaking, uncertainties on δ 11 B sw are approximately 1 ‰, apart from during the Neogene where large uncertainties on data constraints propagate to large uncertainties in our reconstruction.Our higher estimated δ 11 B sw is indicative of lower Miocene pH than previously thought, though we acknowledge high uncertainties during this time.
We see a notable correspondence between δ 11 B 4 and δ 7 Li, something which has previously been used to infer that the evolution of δ 11 B sw was likely similar to that of δ 7 Li (Raitzsch & Hönisch, 2013).However the evolution of δ 11 B sw constrained here shows that the majority of the Cenozoic change in δ 11 B 4 was driven by pH, indicating that links between the controls on pH and δ 7 Li are perhaps the more important control on the similar trajectories of δ 11 B 4 and δ 7 Li.
Uncertainties in our reconstruction are largest where boron isotope data are sparse, such as during the Oligocene, or where datapoints are in disagreement with one another, such as during the Neogene, and we encourage the generation of future records to target these intervals.Looking forward, it would undoubtedly be helpful to find more direct proxies for δ 11 B sw , and to improve constraints on the various models which guide both the data constraints and limits on the rate of change in δ 11 B sw .However, using this method we are able to constrain δ 11 B sw to a range of ±1 ‰ across most of the Cenozoic, improving current estimates.Results from this study can be used to propagate uncertainties in δ 11 B sw into future reconstructions of palaeo pH from boron isotopes and, by extension, palaeo CO 2 .We provide both the metrics which describe our reconstruction of δ 11 B sw (the median, and 95% confidence interval), and also 10,000 statistical samples of possible evolutions of δ 11 B sw .Uncertainty in δ 11 B sw can then be propagated into future data products using a Monte Carlo approach as described in Section 6.In addition we provide analogous information for all signals reconstructed here using the Gaussian Process (δ 11 B 4 , pH, 87/86 Sr, 187/188 Os, and δ 7 Li).

Data Availability Statement
All code used in and produced by this project is stored within our GitHub Repository: https://github.com/St-Andrews-Isotope-Geochemistry/d11Bsw-Gaussian-Processand archived using Zenodo (Whiteford, 2024b).Data files associated with this project are archived using figshare (Whiteford, 2024a).Three data files are provided which include: • Data output for δ 11 B sw in the form of a .xlsxfile, which contains both summary metrics and all 10,000 time series realizations.• Original data, metrics summary of the reconstruction, and 10,000 individual statistical samples for 87/86 Sr, 187/188 Os, and δ 7 Li in a .xlsxfile.
• δ 11 B 4 data input to the Gaussian Process, metrics summary of the reconstruction, and individual statistical samples for δ 11 B 4 and pH in a .xlsxfile.
Two forms of output for δ 11 B sw are given, because while the median and 95% confidence interval give a sense of reasonable values and allow easy plotting, they are unable to convey the covariance embedded in each Gaussian Process sample.Thus for uncertainty propagation in future calculations, we recommend using the time series contained within the .xlsxfile.This methodology allows propagation of both uncertainties in δ 11 B sw and also the rate of change of δ 11 B sw as described in Section 6.Data for δ 11 B 4 can be found as a supplementary table to J. W.

Paleoceanography and Paleoclimatology
10.1029/2023PA004769 WHITEFORD ET AL.Rae et al. (2021).Software used in the project is written in Python, and is available on our GitHub repository: https://github.com/St-Andrews-Isotope-Geochemistry/d11Bsw-Gaussian-Process(Whiteford, 2024b).It uses a software package we've written to perform the statistical calculations, in particular representing distributions, drawing samples, and performing the Gaussian Process interpolation.Scripts to perform the calculation, analyze the output, and display the results are also included.

Figure 1 .
Figure 1.The left panel (a) shows a graph of pH versus δ 11 B 4 for three different δ 11 B sw 's (indicated by lines and corresponding shaded windows).The y axis here can be used for both δ 11 B 4 and δ 11 B sw .A hypothetical 2‰ excursion from 22 to 20‰ is depicted by the black dotted lines.The equivalent pH change for the same excursion at variable δ 11 B sw is displayed in the shaded region.Note this curve is symmetrical about the value of pK * B (roughly 8.6 for standard modern open ocean conditions).The right panel (b.) shows (for the same 22-20‰ excursion as shown in the left panel) ΔpH as a function of δ 11 B sw (pink line).A hypothetical constraint on ΔpH of 0.2 ± 0.04 (at 2σ) is shown by the black dashed lines, with uncertainty in the black dotted lines.Where this hypothetical constraint intersects with the pink line gives the region of possible δ 11 B sw 's.We illustrate that a Monte Carlo approach, which samples possible ΔpH from the hypothetical constraint to give probability distributions (shown in pink shaded regions) for δ 11 B sw , aligns with the emplaced constraint.Due to the symmetry of the curves shown in the left panel, there are often two δ 11 B sw 's compatible with a ΔpH constraint.Typically the lower window is rejected as it would resolve to an unreasonably high pH.

Figure 2 .
Figure 2. Schematic of the processes which cause vertical gradients in oceanic properties as harnessed in reconstructions of δ 11 B sw .Oxygen and carbon dioxide concentrations are set in surface waters by their atmospheric concentrations and Henry's Law (which is temperature dependent).Surface ocean carbon is taken up by biomass to form organic tissues, which are preferentially enriched in the lighter carbon isotope.When that biomass is exported to depth and remineralized, organic carbon is oxidized (consuming oxygen) and releasing CO 2 , which drives acidification.These processes cause correlated gradients in pH, aqueous CO 2 , aqueous O 2 , saturation state, and δ 13 C.As δ 11 B sw is a well mixed signal, differences in δ 11 B 4 between the surface and subsurface are primarily driven by the gradient in pH.δ 11 B sw is homogenous throughout the water column due to the long 10 Myr residence time of boron(Foster et al., 2010;Lemarchand et al., 2000).

Figure 3 .
Figure 3.An illustrative example of the Gaussian Process approach applied to reconstruct a hypothetical time series (shown in black).Each subplot shows the mean reconstruction (colored line), 95% confidence interval (colored window), and three random Gaussian Process samples (gray lines).Panel (a) shows a Gaussian Process prior with set mean, kernel function, and hyperparameters but without any data constraints.Panel (b) shows a Gaussian Process posterior with the same prescribed mean and hyperparameters as panel a., but with three noisy data constraints (shown with 2σ uncertainties in blue).Panel (c) shows the same Gaussian Process posterior as panel (b), where, in addition to the three noisy data constraints, we have applied additional constraints on the acceptable rate of change of the value through time.

Figure 4 .
Figure 4. Our reconstruction of δ 11 B sw is shown in pink, with a central line depicting the median, and a window representing the 95% confidence interval.Central data constraints are shown by vertical pink bars representative of 95% confidence intervals, and upper limits (obtained based upon the planktic δ 11 B 4 measurements) are shown by horizontal pink bars.The gray line depicts δ 11 B sw reconstructed by J. W. Rae et al. (2021) for comparison.δ 11 B 4 datapoints from J. W. Rae et al. (2021) are shown as purple dots (interpolated in the purple shaded region using a Gaussian Process), which are converted to pH using our δ 11 B sw (red points, line, and shaded region) and the δ 11 B sw from J. W.Rae et al. (2021) (gray points and line) for comparison.We note that our Gaussian Process interpolation of δ 11 B 4 (and therefore also pH) has difficulty during large data gaps.In particular during the Palaeogene (which we have faded out), the large data gap is bounded by events with rapid changes in δ 11 B 4 , and our reconstruction would be tempered by filling in this region with δ 11 B 4 data.

5. 1 .
Decomposition of δ 11 B 4 δ 11 B 4 is an integrated signal, combining the effects of δ 11 B sw , pH, pK * B , and ɛ.Of these, only a few are timevariable: δ 11 B sw , pH, and some of the factors which influence pK * B -in particular temperature, and seawater elemental composition.The sensitivity of δ 11 B 4 to these factors (illustrated in Figure S2 in Supporting Information S1) is such that only changes in δ 11 B sw or pH would be of sufficient magnitude to drive the observed 8 ‰ change in δ 11 B 4 across the Cenozoic.Temperature may also have driven some portion of the change, though estimating how much is complicated because individual sample sites may deviate significantly from global average temperature.All other factors being equal, a temperature change of +10°C drives approximately 1.5 ‰ change in δ 11 B 4 .Temperature change is taken into account in all calculations presented here,.and if we adjust for any secular evolution of temperature we might therefore simplify the δ 11 B 4 record into two components, δ 11 B sw and pH.Any change in δ 11 B 4 not explained by changes in δ 11 B sw must be the result of pH, and vice versa.Given the observed 8 ‰ change in δ 11 B 4 , our reconstructed 2 ‰ change in δ 11 B sw , and temperature change across the Cenozoic of approximately 12°C driving approximately 1.8‰ of change, this suggests roughly 4 ‰ of change has been driven by changing pH.

Figure 5 .
Figure 5.Our reconstruction of δ 11 B sw is shown pink, with a thick central line depicting the median, and a window representing the 95% confidence interval.Thinner pink lines give a sense of the covariance structure in our posterior estimate of the δ 11 B sw .Data constraints are shown by the pink dots (here displayed without uncertainty for clearer comparison of signal trends).Other related geochemical signals δ 7 Li, 87/86 Sr, and 187/188 Os are shown in yellow, green, and blue respectively afterMisra and Froelich (2012).Data constraining 87/86 Sr and δ 7 Li is sourced fromMisra and Froelich (2012), while data constraining187/188 Os is compiled from a number of sources (listed in the main text).We interpolate each of these signals with a Gaussian Process with hyperparameters guided by the residence time of each signal.

Table 1 Table of Data
Constraints With Representative Central Values and 95% Confidence Intervals as an Indicator of Uncertainty Anagnostou et al. (2016)em calculations and CO 2 estimates for this time are provided in Supporting Information S1, and are plotted inHenehan and Witts (2023).Datapoints fromAnagnostou et al. (2016)are updated by modifying the calculation of vital effects to account for changing seawater chemistry.With these recalculations, Henehan et al. (2019)'s estimate of δ 11 B sw is updated from the originally published value of 39.05-39.85‰ to 39.05-39.55‰.Finally, estimates fromAnagnostou et al. (