Orbital CO₂ reconstruction using boron isotopes during the late Pleistocene, an assessment of accuracy

Boron isotopes in planktonic foraminifera are a widely used proxy to determine ancient surface seawater pH and by extension atmospheric CO₂ concentration and climate forcing on geological timescales. Yet, to reconstruct absolute values for pH and CO₂, we require a δ¹¹B_foram-borate to pH calibration and independent determinations of ocean temperature, salinity, a second carbonate parameter, and the boron isotope composition of seawater. Although δ¹¹B-derived records of atmospheric CO₂ have been shown to perform well against ice-core-based CO₂ reconstructions, these tests have been performed at only a few locations and with limited temporal resolution. Here we present two highly resolved CO₂ records for the late Pleistocene from Ocean Drilling Program (ODP) Sites 999 and 871. Our δ¹¹B-derived CO₂ record shows a very good agreement with the ice core CO₂ record with an average offset of 13±46 (2σ) and an RMSE of 26 ppm, with minor short-lived overestimations of CO₂ (of up to ∼50 ppm) occurring during some glacial onsets. We explore potential drivers of this disagreement and conclude that partial dissolution of foraminifera has a minimal effect on the CO₂ offset. We also observe that the general agreement between δ¹¹B-derived and ice core CO₂ is improved by optimising the δ¹¹B_foram-borate calibration. Despite these minor issues, a strong linear relationship between relative change in climate forcing from CO₂ (from ice core data) and pH change (from δ¹¹B) exists over the late Pleistocene, confirming that pH change is a robust proxy of climate forcing over relatively short (<1 million year) intervals. Overall, these findings demonstrate that the boron isotope proxy is a reliable indicator of CO₂ beyond the reach of the ice cores and can help improve determinations of climate sensitivity for ancient time intervals.

1 Introduction

The boron isotope composition of ancient planktonic foraminifera shells is widely used to reconstruct past concentrations of atmospheric CO₂ to understand the drivers and responses of climate change over orbital and geological timescales. Unlike many environmental proxies where it is difficult to assess the accuracy of the resulting reconstructions (e.g. for sea surface temperature), the boron isotope pH and CO₂ proxy can directly be compared with the ice core CO₂ records, i.e. the West Antarctic Ice Sheet divide (Ahn et al., 2012), the EPICA (European Project for Ice Coring in Antarctica) Dome Concordia ice core record (Siegenthaler et al., 2005; Lüthi et al., 2008; Bereiter et al., 2015), and the Vostok ice core record (Petit et al., 1999). This comparison of CO₂ over the last 800 kyr provides a very powerful test of proxy accuracy. Several past intervals have been studied to test the boron isotope proxy in this way (Sanyal et al., 1995; Foster, 2008; Hönisch and Hemming, 2005; Henehan et al., 2013; Raitzsch et al., 2018).

Given the success of these comparisons, the boron isotope proxy has been used to investigate the interaction between CO₂, the ocean carbon cycle, and climate beyond the reach of the ice cores, such as during the Mid-Pleistocene transition (Hönisch et al., 2009; Chalk et al., 2017; Dyez et al., 2018), the Pliocene (Martínez-Botí et al., 2015; De La Vega et al., 2020), the Miocene (Foster et al., 2012; Greenop et al., 2017; Guillermic et al., 2022), the Eocene (Anagnostou et al., 2016, 2020; Harper et al., 2020), the Paleocene–Eocene boundary (Penman et al., 2014; Gutjahr et al., 2017), and the Cretaceous–Palaeogene boundary (Henehan et al., 2019). Application of the boron isotope proxy is, however, complicated by the need for (i) an empirical species-specific calibration of δ¹¹B_foraminifera to δ¹¹B_borate in the pH expression (Henehan et al., 2013, 2016, hereafter δ¹¹B_foram-borate calibration), sometimes including extinct species for deep-time reconstruction; (ii) δ¹¹B of seawater (δ¹¹B_sw), temperature, and salinity in the past to calculate pH from δ¹¹B; and (iii) a second carbonate parameter (typically total alkalinity, dissolved inorganic carbon, DIC, or calcite saturation state) to convert pH to CO₂. While these variables do not influence the magnitude of uncertainty equally in all time intervals, assessment of the boron-based reconstructions against existing ice core records is a powerful test of the proxy's accuracy.

Recently, Hain et al. (2018) suggested that the radiative forcing from CO₂ change (ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$) is linearly related to pH change (ΔpH) of equilibrated water of the low-latitude surface ocean when the CO₂ change occurs faster than the residence time of carbon with respect to silicate weathering (e.g. ∼1 Myr). That is, glacial or interglacial CO₂ climate forcing could be estimated directly from reconstructed ΔpH. Given that one of the main priorities for accurate reconstructions of past CO₂ levels is to allow determinations of climate sensitivity, defined as the temperature response to a radiative forcing – typically a doubling of CO₂ with associated slow and fast feedbacks (e.g. Rohling et al., 2018; PALAEOSENS Project Members, 2012) – this recognition may provide a useful shortcut. Climate forcing is a perturbation of the planet's energy balance averaged over the planet (Hansen et al., 2008), and CO₂ forcing, ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ expressed in W m⁻², at a given time can be written as follows:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathrm{F}}_{{\mathrm{CO}}_{\mathrm{2}}}\cong {\mathit{\alpha }}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}\cdot {\displaystyle \frac{\mathrm{\Delta }{\log }_{\mathrm{10}}{\mathrm{CO}}_{\mathrm{2}}}{{\log }_{\mathrm{10}}\mathrm{2}}},\end{array}$$

where ${\mathit{\alpha }}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ is the sensitivity of the radiative balance per doubling of CO₂ and Δlog ₁₀CO₂ is the CO₂ change over time expressed in terms of how many 10-foldings of proportional (not absolute) CO₂ change (Hain et al., 2018).

By considering basic equilibrium reactions of carbon species, Δlog ₁₀CO₂ can be derived and expressed as follows:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{\log }_{\mathrm{10}}{\mathrm{CO}}_{\mathrm{2}}\cong \mathrm{\Delta }{\log }_{\mathrm{10}}\mathrm{DIC}+\mathrm{\Delta }p{K}_{\mathrm{0}}+\mathrm{\Delta }p{K}_{\mathrm{1}}-\mathrm{\Delta }\mathrm{pH}.\end{array}$$

Hain et al. (2018) showed that the terms Δlog₁₀DIC and ΔpK₀+ΔpK₁ are small and that ΔlogCO₂ can therefore simply be expressed as follows:

$$\begin{array}{}\text{(3a)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{\log }_{\mathrm{10}}{\mathrm{CO}}_{\mathrm{2}}\cong -\mathrm{\Delta }\mathrm{pH},\text{(3b)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{\mathrm{F}}_{{\mathrm{CO}}_{\mathrm{2}}}\cong -{\displaystyle \frac{{\log }_{\mathrm{10}}\mathrm{2}}{{\mathit{\alpha }}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}}\mathrm{\Delta }\mathrm{pH}\cong -\mathrm{12.3}\mathrm{\Delta }\mathrm{pH}.\end{array}$$

To assess the uncertainty of this approximate $-\mathrm{1}:\mathrm{1}$ Δlog ₁₀CO${}_{\mathrm{2}}/\mathrm{\Delta }$pH relationship, Hain et al. (2018) considered three different end-member causes to compute the accurate Δlog ₁₀CO${}_{\mathrm{2}}/\mathrm{\Delta }$pH relationship: (1) DIC addition or removal yields a slope of $-\mathrm{1.3}:\mathrm{1}$ (relative to the basic formalism), (2) CaCO₃ addition or removal (e.g. precipitation or dissolution, riverine input) yields a slope of $-\mathrm{0.9}:\mathrm{1}$, and (3) warming or cooling yields a slope of $-\mathrm{1.1}:\mathrm{1}$. That is, even if ΔpH was known exactly this range of plausible slopes results in estimated Δlog ₁₀CO₂ and ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ that are systematically biased by −10 % for change caused purely by CaCO₃ variations or +30 % for change purely caused by DIC variations relative to the approximate $-\mathrm{1}:\mathrm{1}$ Δlog ₁₀CO${}_{\mathrm{2}}/\mathrm{\Delta }$pH relationship. While introducing such structural uncertainty in the estimation of ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ is a concern, this approach eliminates the need to assume a second carbonate system parameter and the uncertainty incurred thereby. An estimate of δ¹¹B_sw is still needed to reconstruct pH based on the boron isotope proxy system (Foster and Rae, 2016) but estimated pH change (i.e. ΔpH) is much less sensitive to error in assumed δ¹¹B_sw than absolute pH (Hain et al., 2018). An important caveat to estimating ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ directly from ΔpH is that the intercept of the Δlog ₁₀CO${}_{\mathrm{2}}/\mathrm{\Delta }$pH relationship can change with silicate weathering carbon cycle dynamics thought to be important on a million-year timescale, such that the approach is applicable for orbital timescale variability and short-term shifts but not for long-term trends in ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$. Therefore, the orbital timescale ice age cycles of atmospheric CO₂ reconstructed from air occluded in Antarctic ice cores offer a unique opportunity to determine the Δlog ₁₀CO${}_{\mathrm{2}}/\mathrm{\Delta }$pH relationship observationally and compare to theory. Furthermore, given the principal drivers of the glacial–interglacial CO₂ cycles (e.g. change in water masses, sea ice cover, the soft tissue pump, the solubility pump, the CaCO₃ counter pump, and the disequilibrium pump; see Sigman et al., 2010; Hain et al., 2010, 2013, for a full review) will impact the Δlog ₁₀CO${}_{\mathrm{2}}/\mathrm{\Delta }$pH relationship in different ways, comparing the slope of the regressed ΔF$/\mathrm{\Delta }$pH line from data to theoretical end-members (temperature, DIC, CaCO₃) could allow the primary controlling mechanisms during glacial–interglacial (G-IG) cycles to be deciphered.

In light of these recent advances, our aims here are twofold. First, we extend previous ice core validation studies (Foster, 2008; Henehan et al., 2013; Chalk et al., 2017) and test the extent to which boron isotopes reconstruct CO₂ faithfully when current methods and assumptions are applied. In contrast to most previous studies, we use two deep-ocean sites and present δ¹¹B and CO₂ data at high temporal resolution (one sample every ∼3 to 6 kyr). This enables (i) a thorough test of the assumptions typically made including the central tenet of atmospheric CO₂ proxies that surface ocean CO₂ remains in equilibrium with the atmosphere over time at any given site; (ii) an evaluation of the overall uncertainty of the proxy; (iii) an evaluation of the influence of variable foraminiferal preservation on the accuracy of the CO₂ reconstructed; and (iv) a refinement of a number of the input assumptions and uncertainties, including the δ¹¹B_{foram–borate} calibration. Second, we evaluate the approach of Hain et al. (2018) and assess the robustness of pH change to provide insights into not only the magnitude of climate forcing from CO₂ change but also the ability of this approach to explain the causes of CO₂ change over glacial–interglacial cycles.

2 Methods

2.1 Core location and oceanographic setting

To accurately reconstruct atmospheric CO₂ with the δ¹¹B–CO₂ proxy, it is essential to measure δ¹¹B in foraminifera from locations where the CO₂ flux between the ocean and the atmosphere is in near equilibrium. We therefore target regions of the ocean where the water column is stratified and oligotrophic, as these regions are most likely to attain this condition (Takahashi et al., 2009). Here, following previous studies (Foster, 2008; Henehan et al., 2013; Chalk et al., 2017), we report and add new data from Ocean Drilling Program (ODP) Site 999 (Fig. 1, 12.75^∘ N, 78.73^∘ W, water depth 2827 m, sedimentation rate 3.7 cm kyr⁻¹) in the Caribbean and supplement this well-studied site with samples from ODP Site 871 in the western Pacific (5.55^∘ N, 172.35^∘ E, water depth 1255 m, sedimentation rate ∼1 cm kyr⁻¹). The sediments studied at ODP Site 871 are shallowly buried, and the site today features a deep thermocline and is located off the Equator; hence, these sediments are unlikely to be influenced by significant equatorial upwelling (Dyez and Ravelo, 2013, 2014). These two sites show a minor annual mean disequilibrium of +12 ppm (range ∼0 to ∼30 ppm, Takahashi et al., 2009) for ODP Site 871 and +21 ppm (Olsen et al., 2004; Foster, 2008) for ODP Site 999. These disequilibria are used to correct our CO₂ data derived from δ¹¹B and are assumed to be constant throughout the entire record presented here (with an uncertainty of ±10 ppm).

Figure 1 Map of air–sea CO₂ disequilibrium (seawater–air, in ppm) and location of ODP sites used in this study. CO₂ data are from Takahashi et al. (2009). The map was made with Ocean Data View (Schlitzer, 2023).

While we recognise that both sites have a minor disequilibrium, this is often a necessary compromise as areas of the ocean that are in strict equilibrium with the atmosphere are often located in the middle of oceanic gyres and tend to have deep sediments located under the lysocline, have a low sedimentation rate, and/or are outside the preferred geographic habitat of G. ruber. Furthermore, we present surface δ¹⁸O and δ¹³C (site 871) and temperature (both sites) from G. ruber that provide insight into the potential influence of upwelling (see Sect. 4.2.2) at these locations. Recent Earth system model (IPSL-CM5A-MR) outputs (Gray and Evans, 2019) also show that relative pH difference at our core sites between the Last Glacial Maximum (LGM) and the pre-industrial period (PI) compared to the ocean average pH difference are close to 0, giving confidence that changes in local disequilibrium are unlikely to drive large changes in our CO₂ reconstructions (at least during the last glacial period).

2.2 Samples

2.2.1 Sample selection and preparation

Samples of deep-sea sediment from our two study sites were taken at 6 cm (∼3 kyr) and 10 cm (∼6 kyr) resolution at ODP 871 and 999, respectively. Around 1–2 mg of the foraminifer (between 120 and 200 individuals) from the species Globigerinoides ruber sensu stricto white (hereafter G. ruber ss) were hand-picked from the size fraction 300–355 µm for a target of 10 to 20 ng of boron. G. ruber ss was chosen here because it is readily identified; is abundant throughout our chosen time interval; and a δ¹¹B_foram-borate calibration that accounts for vital effects is available from culture, plankton tows, and core-top samples (Henehan et al., 2013). It is also known to live in the upper surface of the ocean with a relatively small depth range (Rebotim et al., 2017), which prevents significant influence of deeper, more remineralised CO₂-rich waters on the measured δ¹¹B. The morphotype G. ruber sensu lato (hereafter G. ruber sl) has slightly different morphology (Aurahs et al., 2011; Carter et al., 2017) and is thought to live in deeper water compared to G. ruber ss (Wang, 2000). The morphotype G. ruber sl was also hand separated and analysed at lower resolution at ODP 871 to monitor any change over time in morphotype differences in δ¹¹B that could result from different habitats. For similar reasons, carbon and oxygen isotopes (δ¹⁸O and δ¹³C) were also measured on G. ruber ss and sl for comparison on the whole record at ODP 871. For this, around 10 individuals of G. ruber per sample were picked, their shells gently broken open and mixed and then a 100 µg aliquot of the homogenised carbonate was measured using a Thermo KIEL IV Carbonate device at the University of Southampton, Waterfront Campus. While this number of specimens is lower than classically done for δ¹⁸O and δ¹³C analysis, it provides power for the identification of species-specific preferential diagenetic alteration, which may have occurred in the sediment, and it was sometimes necessary due to the scarcity of some of the G. ruber spp. morphotypes.

2.2.2 Age constraints

Samples were taken from 1.5 to 5 m b.s.f. (metres below sea floor) for ODP 871 and from 9 to 21 m b.s.f. for ODP 999. Sample age at Site 871 was initially determined from sample depth using published age models (Dyez and Ravelo, 2013). At Site 999, the age was determined by developing a new Cibicidoides wuellerstorfi benthic δ¹⁸O record. The initial age model at Site 871 was refined by measuring δ¹⁸O on the benthic species Uvigerina peregrina (50 µg of 3–5 mixed, crushed, and homogenised specimens) measured on a Thermo KIEL IV Carbonate device at the University of Southampton, Waterfront Campus. These new δ¹⁸O data (Fig. 2) were then tuned to the benthic δ¹⁸O LR04 stack (Lisiecki and Raymo, 2005) using Analyseries (Paillard et al., 1996). A correction of +0.47 was applied to the δ¹⁸O Cibicidoides wuellerstorfi at ODP Site 999 following Marchitto et al. (2014).

Figure 2 $\mathrm{Mg}/\mathrm{Ca}$-derived temperature, coarse fraction (sand), fragmentation, and benthic δ¹⁸O and δ¹³C at ODP sites 999 and 871. (a, b) Temperature at ODP 999 (from G. ruber ss, black, Schmidt et al., 2006, and this study) and ODP 871 (G. ruber ss, black, G. ruber sl, red, 2 SD indicated by the grey error bar). (c, d) Fragmentation index (light grey, data from Schmidt et al., 2006, for ODP 999) and sand (black line). (e, f) Benthic C. wuellerstorfi (Site 999) and U. peregrina (Site 871) δ¹⁸O (blue) and LR04 benthic δ¹⁸O stack (black). A correction of +0.47 ‰ is applied to δ¹⁸O of C. wuellerstorfi data to adjust for species offset. (g, h) Benthic C. wuellerstorfi (Site 999) and U. peregrina (Site 871) δ¹³C (orange).

2.2.3 Fragment counts

Foraminifera fragment counts were conducted on ODP Site 871 to monitor variations in carbonate preservation. Samples were sub-sampled using a splitter (in order to maintain homogeneity) and poured onto a picking tray. The fragmentation index (FI) was calculated following the approach of Howard and Prell (1994) and Berger (1970) where percentage fragment is defined as follows:

$$\begin{array}{}\text{(4)}& \mathrm{FI}=\mathrm{100}\cdot {\displaystyle \frac{\mathrm{number}\mathrm{of}\mathrm{fragments}}{\mathrm{number}\mathrm{of}\mathrm{fragments}+\mathrm{number}\mathrm{of}\mathrm{whole}\mathrm{tests}}}.\end{array}$$

Counts of whole intact grains and fragments of grains were conducted three times and averaged. The standard deviation (1σ) of the fragmentation index is 1.69. This approach followed that used in an early study at ODP Site 999 (Schmidt et al., 2006), ensuring that the data sets between the two sites are comparable.

2.2.4 Boron separation

The hand-separated foraminifera tests for boron isotope analysis were broken open, detrital clay was removed, and oxidatively cleaned and leached in a weak acid to obtain a primary carbonate signal using established methods (Barker et al., 2003). Samples were then slowly dissolved in ∼100 µL 0.5M HNO₃ added to 200 µL of MQ water. Dissolved samples were then centrifuged for 5 min to separate any remaining undissolved contaminants (e.g. silicate grains, pyrite crystals) and transferred to screw top 5 mL Teflon pots for subsequent boron separation. An aliquot equivalent to 7 % of each sample was kept for elemental analysis and transferred to acid-cleaned plastic vials in 130 µL 0.5M HNO₃. Samples were purified for boron using anion exchange column chemistry method prior to isotope analysis as described elsewhere (Foster, 2008). A total procedure blank (TPB) was conducted for each batch of samples and typically ranged from 0 to 100 pg, which represents a blank contribution of up to 2.3 % (for samples containing ∼10–20 ng of boron). Most samples had a TPB below 40 pg and were not corrected. Two batches had a TPB of 70 and 100 pg, for which we corrected using a long-term median TPB δ¹¹B value of −7.27 ‰ from the University of Southampton. This represents a δ¹¹B correction of 0.1 ‰ to 0.7 ‰.

2.3 Analytical techniques

Boron isotope analyses were performed on a Thermo Scientific Neptune multi-collector inductively coupled plasma mass spectrometer (MC-ICPMS) with 10¹²Ω amplifier resistors using a standard-sample bracketing routine with NIST 951 boric acid standard (following Foster et al., 2013, and Foster, 2008). Elemental analysis was performed on each dissolved sample using a Thermo Scientific Element inductively coupled plasma mass spectrometer (ICPMS). All analyses were carried out at the University of Southampton, Waterfront Campus (following Foster, 2008, and Henehan et al., 2015). Element-to-calcium ratios were measured with ⁴³Ca and ⁴⁸Ca and measured against in house mixed element standards. Elemental ratios measured included $\mathrm{B}/\mathrm{Ca}$, $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Al}/\mathrm{Ca}$, $\mathrm{Mn}/\mathrm{Ca}$, and $\mathrm{Sr}/\mathrm{Ca}$. Based on the reproducibility of our in-house standards, the uncertainty for most elemental ratios is ∼5 % (at 95 % confidence).

2.4 Constraints on δ¹¹B-derived pH and CO₂

2.4.1 From δ¹¹B to pH

Seawater pH is related to the boron isotopic composition of dissolved borate ion by the following equation:

$$\begin{array}{}\text{(5)}& \mathrm{pH}=p{K}_{\mathrm{B}}-\log \left(-{\displaystyle \frac{{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{sw}}-{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{borate}}}{{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{sw}}-a_{\mathrm{B}}\cdot {\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{borate}}\left(a_{\mathrm{B}}-\mathrm{1}\right)}}\right),\end{array}$$

where the isotopic fractionation factor α_B between B(OH)₃ and B(OH)${}_{\mathrm{4}}^{-}$ is 1.0272 as determined by Klochko et al. (2006) and the δ¹¹B of seawater is 39.61 ‰ (Foster et al., 2010) for both sites and kept constant throughout the record due to the long residence time of boron (10–20 Myr, Lemarchand et al., 2002).

The sea surface temperature (SST) values necessary to calculate pK_B in Eq. (5) were determined at both sites using the $\mathrm{Mg}/\mathrm{Ca}$ of G. ruber (Dyez and Ravelo, 2013) including a depth-dependent dissolution correction for each site (following Dyez and Ravelo, 2013, for Site 871 and Schmidt et al., 2006, for Site 999) and a pH correction using the iterative approach of Gray and Evans (2019) to account for the observed pH effect on $\mathrm{Mg}/\mathrm{Ca}$ in G. ruber producing higher apparent sensitivity of $\mathrm{Mg}/\mathrm{Ca}$ during glacial cycles (Gray et al., 2018).

$\mathrm{Mg}/\mathrm{Ca}$ was corrected for depth-dependent dissolution at Site 871 using the following equation (Dyez and Ravelo, 2013):

$$\begin{array}{}\text{(6a)}& {\displaystyle \frac{\mathrm{Mg}}{\mathrm{Ca}}}\left(\mathrm{corrected}\right)={\displaystyle \frac{\mathrm{Mg}}{\mathrm{Ca}}}\left(\mathrm{measured}\right)+\mathrm{0.26}\cdot \mathrm{depth}+\mathrm{0.52}.\end{array}$$

$\mathrm{Mg}/\mathrm{Ca}$ from Site 999 was corrected following Schmidt et al. (2006):

$$\begin{array}{}\text{(6b)}& {\displaystyle \frac{\mathrm{Mg}}{\mathrm{Ca}}}\left(\mathrm{corrected}\right)={\displaystyle \frac{\mathrm{Mg}}{\mathrm{Ca}}}\left(\mathrm{measured}\right)+\mathrm{0.66}.\end{array}$$

To evaluate the effect of various $\mathrm{Mg}/\mathrm{Ca}$ treatment on temperature and calculated CO₂, we performed seven sensitivity tests (Table S1) with $\mathrm{Mg}/\mathrm{Ca}$-derived SST using the calibrations of (1) Gray et al. (2018) temperature-dependent only (global calibration), (2) Gray and Evans (2019) with a pH correction, (3) Gray et al. (2018) temperature-dependent with $\mathrm{Mg}/\mathrm{Ca}$ corrected for depth-dependent dissolution, (4) Gray and Evans (2019) with $\mathrm{Mg}/\mathrm{Ca}$ corrected for depth-dependent dissolution and pH correction, (5) Anand et al. (2003) with and without a depth correction, and (6) with temperature kept constant (26 ^∘C).

The differences in SST and resulting CO₂ can be substantial (Fig. S1, Table S2): up to 6^∘ and ∼50 ppm, respectively, between the Gray et al. (2018) calibration uncorrected for pH and the Anand et al. (2003) calibration corrected for dissolution. We have chosen the $\mathrm{Mg}/\mathrm{Ca}$ treatments that account for pH effect on $\mathrm{Mg}/\mathrm{Ca}$ and yield the closest agreement between coretop at both sites and modern temperature from Glodap v2 (Lauvset et al., 2022). This treatment is with a pH correction and $\mathrm{Mg}/\mathrm{Ca}$ corrected for depth-dependent dissolution. Choosing this approach is justified considering (1) the strong offset between Anand et al. (2003) multi-species $\mathrm{Mg}/\mathrm{Ca}$–temperature calibration and the more recent G. ruber compilation of Gray et al. (2018), (2) the effect of pH correction as shown in Gray et al. (2018) and Gray and Evans (2019), (3) the suggested influence of dissolution on $\mathrm{Mg}/\mathrm{Ca}$ (Dyez and Ravelo, 2013; Schmidt et al., 2006), and (4) the better agreement between coretop and modern SST at each site when using a pH and depth correction (Fig. S1).

The salinity (S) that is used in the expression of pK_B is kept constant for both sites (35) due to the very minor effect of salinity on calculated pH/CO₂ (1 salinity unit changes pH by 0.006).

2.4.2 From pH to CO₂

Calculating CO₂ from boron-isotope-derived pH is dependent on the determination of a second parameter of the carbonate system. Here we use the modern value of total alkalinity (TA) at each site: 2279 and 2350 µmol kg⁻¹ at ODP 871 and ODP 999, respectively (Shipboard Scientific Party, 1993; Takahashi et al., 2009). Following Chalk et al. (2017), these values were kept constant throughout the whole record. To account for any variations in alkalinity, a generous uniform (i.e. equal likelihood of values within the range of uncertainty) uncertainty of 175 µmol kg⁻¹, distributed equally on either side of the central value, is applied. This range in TA encompasses the likely range in this variable on glacial–interglacial (e.g. Toggweiler, 1999; Hain et al., 2010; Cartapanis et al., 2018) or longer timescales (Hönisch et al., 2009), and its adoption means the local TA record is not tied to a global sea level record as has been practised previously. We avoid drawing this link because the $\sim +\mathrm{3}$ % (+68 µmol kg⁻¹) concentration increase in solute alkalinity occurring from sea level lowering during the Last Glacial Maximum may not have been the dominant driver of ocean alkalinity change (Boyle, 1988a, b; Sigman et al., 1998; Toggweiler, 1999; Hain et al., 2010; Cartapanis et al., 2018). By assuming a uniform distribution for TA, we avoid imposing a temporal evolution to this variable because evolution of TA through a glacial cycle is uncertain and is unlikely to be simply a function of sea level or salinity (e.g. Dyez et al., 2018) due to the effect of carbonate compensation.

The surface water CO₂ is then calculated as follows (Zeebe and Wolf-Gladrow, 2001):

$$\begin{array}{}\text{(7)}& {\mathrm{CO}}_{\mathrm{2}}={\displaystyle \frac{\mathrm{TA}-\frac{K_{\mathrm{B}}\cdot B_{\mathrm{T}}}{K_{\mathrm{B}}+\left[{\mathrm{H}}^{+}\right]}-\frac{K_{\mathrm{W}}}{\left[{\mathrm{H}}^{+}\right]}+\left[{\mathrm{H}}^{+}\right]}{\frac{K_{\mathrm{1}}}{\left[{\mathrm{H}}^{+}\right]}+\frac{\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{{\left[{\mathrm{H}}^{+}\right]}^{\mathrm{2}}}}},\end{array}$$

where TA is the total alkalinity, K_B the equilibrium constant of boron species in seawater, B_T the concentration of boron in seawater (432.6 µmol kg⁻¹, Lee et al., 2010), [H⁺] the concentration of H⁺ determined from pH = −log[H⁺], K_W the dissociation constant of water (function of T, S, and pressure), and K₁ and K₂ the first and second dissociation constants of carbonic acid (function of T, S, and pressure, Luecker et al., 2000). The estimate of atmospheric CO₂ includes site-specific offsets relative to reconstructed surface water CO₂ to account for observed local disequilibrium (+21 and +12 ppm at ODP Sites 999 and 871, respectively).

2.5 Uncertainty

2.5.1 Analytical uncertainty

The uncertainty of the measured δ¹¹B is expressed as the external uncertainty, which includes instrumental error and chemical separation of the sample (see a detailed discussion in John and Adkins, 2010). This was determined empirically by long-term repeat measurements of JCp-1 subject to the same chemical purification as our foraminiferal samples. As discussed by Rae et al. (2011), this uncertainty is dependent on the intensity of the ¹¹B signal and is expressed here by the following relationship defined during the duration of this study at the University of Southampton (Anagnostou et al., 2019) for ¹¹B intensities < 0.54 V:

$$\begin{array}{}\text{(8)}& \mathrm{2}\mathit{\sigma }=\mathrm{129}\mathrm{600}{e}^{-\mathrm{212}}\times \left[{}^{\mathrm{11}}\mathrm{B}\right]+\mathrm{0.3385}{e}^{-\mathrm{1.544}}\times \left[{}^{\mathrm{11}}\mathrm{B}\right],\end{array}$$

where [¹¹B] is the intensity of ¹¹B signal in volts. The δ¹¹B uncertainty for ¹¹B intensities > 0.54 V is 0.15 ‰ (at 95 % confidence).

2.5.2 The pH and CO₂ uncertainty

The CO₂ uncertainty we report was calculated with a Monte Carlo simulation (10 000 realisations) in order to fully account for the uncertainty in all variables used in the calculation of pH and CO₂ (${\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}{\mathit{\delta }}^{\mathrm{11}}\mathrm{B}\text{-}\mathrm{derived}}$). The shape of the uncertainty distribution sampled is either normally distributed (for temperature, salinity, and δ¹¹B) or uniform (for alkalinity, as discussed above). The maximum probability of all realisations was used as the central value for CO₂, and an error envelope at 1 and 2σ was calculated based on the 68 % and 95 % distribution of the realisations.

2.5.3 Uncertainty of the CO₂ offset

To constrain the offset between δ¹¹B-derived CO₂ and ice core CO₂, each sediment age is compared to the ice core CO₂ record by interpolation of the record of highest resolution (in this case the δ¹¹B record onto the ice core compilation). To fully account for age uncertainty when interpolating the sediment age to the well-dated ice core record, a distribution of the ice core data was calculated within the 4σ uncertainty of the δ¹¹B age and weighed by the respective likelihood based on the age difference between ice core and sediment core (Hain et al., 2018).

The CO₂ offset (or residual) is defined by the following equation:

$$\begin{array}{}\text{(9)}& {\mathrm{Offset}}_{{\mathrm{CO}}_{\mathrm{2}}}={\mathrm{CO}}_{\mathrm{2}{\mathit{\delta }}^{\mathrm{11}}\mathrm{B}\text{-}\mathrm{derived}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{ice}}.\end{array}$$

The uncertainty of this offset (σ_offset) accounts for the uncertainty of the interpolated ice core CO₂ (${\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}\cdot \mathrm{interpol}}$) and the one of the δ¹¹B-derived CO₂ (${\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}\cdot {\mathit{\delta }}^{\mathrm{11}}\mathrm{B}\text{-}\mathrm{derived}}$), such as the following equation:

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{\mathrm{offset}}=\sqrt{{\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}\cdot \mathrm{interpol}}^{\mathrm{2}}+{\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}\cdot {\mathit{\delta }}^{\mathrm{11}}\mathrm{B}\text{-}\mathrm{derived}}^{\mathrm{2}}}.\end{array}$$

2.6 The relationship between δ¹¹B-derived pH and ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$

The linear relationships between the relative CO₂ forcing ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ and pH are determined with a York regression (York et al., 2004) that accounts for the uncertainty in both the independent and dependent variable (i.e. x and y axes). The ice core CO₂ interpolation used to calculate ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ and uncertainty is determined as described in Sect. 2.6.3 (Hain et al., 2018).

2.7 Optimising the G. ruber δ¹¹B borate-foraminifera calibration

An optimised G. ruber calibration was obtained by minimising the root-mean-square error (RMSE) of the average offset between δ¹¹B-derived CO₂ and ice core CO₂. The steps are illustrated in Fig. S2. In order to optimise the calibration, 10 000 simulations of δ¹¹B_borate and δ¹¹B_foraminifera from the calibration of Henehan et al. (2013) were performed within their normally distributed uncertainty (1σ), from which we defined the same number of linear models each including their slope and intercept. We then calculate the equilibrium pH and resultant equilibrium δ¹¹B_borate from ice core CO₂ and the assumed constant TA at each core site. The δ¹¹B_borate from the 10 000 linear models is then calculated, and the difference from the ice core-derived δ¹¹B_borate is determined. The linear model calibration that yields the minimum RMSE between these two borate variables defines the new δ¹¹B_{foram–borate} calibration. To assess the effect of δ¹¹B records from different sites, we performed this exercise using the combined records (from both sites 999 and 871), 999 only, and 871 only (Fig. S3) and show that using a record from one particular site or the combination of sites yields similar CO₂ offsets (Table S3), and thus here we use the results from the combined sites. Unless indicated otherwise, to preserve a degree of independence, the pH results presented in this study are calculated with the published calibration (Henehan et al., 2013), and the results with the optimised calibration are presented in Sect. 4.2.6.

3 Results

3.1 Temperature and fragment counts

The SST at ODP Sites 999 and 871 show a cyclicity that agrees with the well-known glacial–interglacial cycles of the late Pleistocene (Fig. 2). The SST determined from G. ruber sl (filled red circles, Fig. 2b) at Site 871, show equal or cooler temperatures (by 1–2 ^∘C) than G. ruber ss (filled black circles). The fragmentation index (Fig. 2d) at ODP 871 ranges from 20 % to 50 % and follows the well-documented “Pacific style” dissolution cycles (Sexton and Barker, 2012) with well-preserved carbonate (low fragments) during glacials and less well-preserved carbonates (higher fragments) during interglacials. The percentage sand typically anti-correlates with fragmentation counts at both sites, although it is less clear at ODP 999, perhaps due to the shorter record available. Fragmentation counts reach maxima at ODP 999 of 20 % during interglacials and up to 50 % during marine isotope stage (MIS) 11 at 400 ka, which is concomitant with the Mid-Brunhes Dissolution Interval (MBDI, Barker et al., 2006). The fragmentation counts at ODP 871 show no substantive anomaly during the MBDI.

3.2 The pH and CO₂ reconstructions

The δ¹¹B, pH, and δ¹¹B-derived absolute CO₂ (Fig. 3) from Sites 871 and 999 show clear cyclicity related to glacial–interglacial cycles. The CO₂ values carry an average uncertainty of ±48 ppm, and the mean offset from the ice core CO₂ for a combination of the two records is 13±46 (2σ) ppm, showing that there is a minor overestimation of CO₂ using the boron method, but it agrees well on average within uncertainty. The RMSE of the CO₂ offset for the combined record is 26 ppm.

Figure 3 The δ¹¹B, pH, and boron-derived CO₂ at site 999 and 871. The δ¹¹B of G. ruber ss and sl (a, b), boron-derived pH (c, d), and CO₂ (e, f) reconstruction from two core locations: ODP 999 (blue, this study and published data, Foster, 2008; Henehan et al., 2013; Chalk et al., 2017) and ODP 871 (orange, this study). The black line in the CO₂ panels is the composite Antarctic ice core CO₂ record (Bereiter et al., 2015). All δ¹¹B-derived data points are from G. ruber ss except the black dots at ODP Site 871 measured on G. ruber sl. Numbers at the bottom of the CO₂ records represent marine isotope stages (black boxes for glacials and white boxes for interglacials). Note the age scale is different between Sites 999 and 871.

Despite the overall close agreement between δ¹¹B-derived CO₂ and ice-core-derived CO₂, each of our δ¹¹B–CO₂ records exhibit some short-lived intervals where the offsets from the ice core record are larger. This is further revealed by the residual CO₂ and the identification of the data above the upper quartile (i.e. the upper 25 % of the data, Fig. S4). Those data do not appear to be randomly distributed and instead occur at ∼100, ∼220–290, and ∼390 ka at ODP Site 999, in all three cases during the early stages of the glaciation (except for the MIS 8 glacial at 280 ka, Fig. S4). The mismatches with the ice core at ODP Site 871 show a similar temporal pattern occurring at ∼220 and ∼300 and ∼350–390 ka (i.e. at glacial inceptions).

3.3 Contrasting δ¹¹B between morphotypes

Within error, the few measurements of δ¹¹B G. ruber sl at ODP 871 all agree with δ¹¹B G. ruber ss (Fig. 3), although the δ¹¹B data of G. ruber sl are higher than G. ruber ss for all four data pairs available. The CO₂ derived from G. ruber sl (Fig. 3) is on average 22 ppm lower than the one derived from G. ruber ss; though the much lower resolution (n=4) impedes a thorough comparison at this stage. The δ¹⁸O and δ¹³C of both morphotypes were compared for the whole records at ODP 871 (Fig. S5), and a cross-plot shows a moderate to good agreement between G. ruber ss and G. ruber sl (r²=0.55 and 0.22 for δ¹⁸O and δ¹³C, respectively, Fig. S6). This is in contrast to other studies (e.g. Wang, 2000; Steinke et al., 2005) that show δ¹⁸O in G. ruber sl to be systematically higher.

3.4 Relationship between δ¹¹B-derived pH and CO₂ forcing from the ice core

A cross plot of δ¹¹B-derived pH CO₂ forcing from the ice core record for each of our marine core study sites is shown in Fig. 4 and is compared to the theoretically derived approximate ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}/\mathrm{\Delta }$pH relationships as adopted by Hain et al. (2018): $-\mathrm{1}:\mathrm{1}$ W m⁻² (dashed black line), CaCO₃ addition or removal ($-\mathrm{0.9}:\mathrm{1}$ W m⁻² plain yellow line), DIC addition or removal ($-\mathrm{1.3}:\mathrm{1}$ W m⁻² dashed–dotted blue line), and warming or cooling temperature forcing ($-\mathrm{1.1}:\mathrm{1}$ W m⁻², dashed red line). Our analysis includes full propagation of uncertainty in pH, in contrast to Hain et al. (2018), who considered only the reported uncertainty of δ¹¹B_borate in their validation exercise. In both cases the uncertainty in ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ accounts for the error in interpolation arising when comparing age-uncertain δ¹¹B-derived pH with ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ from the well-dated and high-resolution ice core CO₂ record (see Sects. 2.7 and 2.6 for details). This treatment of ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ uncertainty is dominated by the spread of ice core CO₂ data points within the δ¹¹B age uncertainty. The data are fitted with a York-type regression (thin black line; York et al., 2004) where the grey envelope represents the uncertainty of the linear relationship that best represents the data (i.e. the envelope is not the prediction interval), considering the uncertainty in pH and ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$. The regressed slope is ΔF$/\mathrm{\Delta }$pH = $-\mathrm{17.2}\pm \mathrm{1}$ W m⁻² ($-\mathrm{1.4}:\mathrm{1}$ relative to basic formalism) and shows a good agreement with the theoretical temperature- and DIC-driven relationships.

Figure 4 Ice-core-based ΔFCO₂ (CO₂ forcing) vs. δ¹¹B-based pH for ODP 999 (filled blue circles, this study and published data from Foster, 2008; Henehan et al., 2013; Chalk et al., 2017) and 871 (filled orange circles). The lines show the relationship between ΔFCO₂ and pH for the simplified formalism (see Sect. 2) ΔFCO${}_{\mathrm{2}}=-\mathrm{12.3}$ ΔpH (dashed black line) and when driven by changes in DIC only (dashed–dotted blue line, ΔF$/\mathrm{\Delta }$pH = −16 W m⁻²), CaCO₃ (yellow line, ΔF$/\mathrm{\Delta }$pH = −11.1 W m⁻²), and temperature T (dashed red , ΔF$/\mathrm{\Delta }$pH = −14.1 W m⁻²). The York regressed line (thin black line and grey shading) falls close to the DIC-driven line (dashed–dotted blue line).

The effect of the uncertainty assigned to pH (fully propagated or using the measurement uncertainty of the boron isotope) on the regressed slope is shown in Fig. S7. The slope of the York regression when using the uncertainty from δ¹¹B only, as in Hain et al. (2018), shows a close agreement with the basic formalism, with a slope of ΔF$/\mathrm{\Delta }$pH = $-\mathrm{13.7}\pm \mathrm{0.3}$ W m⁻², ($-\mathrm{1.1}:\mathrm{1}$ relative to the basic formalism) but with a unsatisfactory goodness of fit (mean squared weighted deviation, MSWD) of 5.3, whereas propagating the full pH uncertainty based on our iterative Monte Carlo simulations improves goodness of fit to ∼0.9 at a Δlog ₁₀CO${}_{\mathrm{2}}/\mathrm{\Delta }$pH of $-\mathrm{1.4}:\mathrm{1}$ (Fig. 4).

4 Discussion

4.1 Cyclicity in foraminifera preservation

Percentage fragments and sand fraction (>63 µm) at both studied core sites are anticorrelated and show a clear cyclicity, with better preservation of carbonates during glacial periods (Fig. 2). The anticorrelation is clearer at ODP Site 871, where we have the longest record (Fig. 2). Preservation in the Pacific (Farrell and Prell, 1989) shows improved (poorer) preservation during glacial (interglacial), and this pattern seems to have originated after the Mid-Pleistocene Transition (Sexton and Barker, 2012). The origin of these cycles could be a combination of enhanced ventilation during glacials in the Pacific (Sexton and Barker, 2012) or increased burial due to enhanced global alkalinity following a decrease in burial in the Atlantic (Cartapanis et al., 2018). However, glacial periods seem to have been accompanied by a diminution in oxygenation in the deep Pacific (Anderson et al., 2019) that may have also impacted preservation.

The observation that the fragmentation records of sites 999 and 871 covary is likely attributable to the different water masses that fill the Caribbean basin relative to the rest of the Atlantic basin. During glacials, the deep Atlantic is filled by nutrient- and carbon-rich corrosive southern-sourced waters (Antarctic Bottom Water) with a reduced contribution from the less corrosive, nutrient-poor North Atlantic Deep Water (Oppo and Lehman, 1993) causing calcareous sediments in the deep Atlantic Ocean > 2500 m to be less well preserved during glacials than interglacials. The opposite pattern of dissolution is seen in the Caribbean because shoaling of the northern-sourced waters during glacials produces a mid-depth well-ventilated water mass that feeds into the Caribbean through its deepest sill (∼1900 m, Johns et al., 2002). Thus, the deep Caribbean is filled with less corrosive waters during glacials than interglacials, improving the preservation of carbonate during glacials in a similar pattern to a Pacific-styled dissolution cycle, albeit in response to Atlantic circulation changes. During interglacials, the northern-sourced waters are mixed with corrosive southern-sourced waters (Antarctic Intermediate Water and Upper Circumpolar Deep Water), leading to sediments that less well preserved.

4.2 Causes of offset between δ¹¹B-derived and ice core CO₂

The δ¹¹B-derived CO₂ record from both of our study sites is in very good agreement with the ice core record, with an average offset for combined both cores of 13±46 (2σ) ppm and corresponding RMSE of 26 ppm. However, the minor CO₂ offsets observed in both records do not appear to be random and tend to fall during the first half of each glacial cycle (Fig. S4). In order to have the highest confidence in CO₂ reconstructions using δ¹¹B, this pattern warrants further investigation (see below).

4.2.1 Comparison between morphotypes of G. ruber

If, as others suggested (e.g. Wang, 2000; Steinke et al., 2005; Numberger et al., 2009), G. ruber sl and G. ruber ss occupied different depth habitats, then inadvertent sampling of the cryptic G. ruber sl morphotype might conceivably produce the biases we observe between δ¹¹B-derived CO₂ and atmospheric CO₂ from the ice cores. However, while our $\mathrm{Mg}/\mathrm{Ca}$-derived temperatures for G. ruber sl and G. ruber ss display variable offsets, they are within uncertainty (Fig. 2) and our δ¹⁸O and δ¹³C data for the two morphotypes at ODP 871 show a good agreement with no consistent differences (Fig. S5). Thus, while the water column profile of δ¹⁸O and δ¹³C can be affected by factors other than temperature, salinity, and biological productivity (e.g, carbonate ion effect, Spero et al., 1997), our data overall suggest that the two morphotypes we analysed shared similar depth habitat preferences.

Henehan et al. (2013) found that G. ruber ss and G. ruber sl record similar δ¹¹B in core-top sediments and through necessity used mixed morphotypes in their culture study. The δ¹¹B-derived pH and CO₂ for G. ruber sl examined here are consistently higher and lower than G. ruber ss by around +0.02 pH units and −22 ppm CO₂ on average, respectively (Fig. 3). This is contrary to expectation if G. ruber sl lived in deeper and more acidic waters as suggested by other studies (Wang, 2000; Steinke et al., 2005) but consistent with some data sets that show that the habitat of G. ruber ss and G. ruber sl can vary by location and seems to be dependent on local productivity (Numberger et al., 2009). Other data sets from the Atlantic and Indian oceans nevertheless show similar $\mathrm{Mg}/\mathrm{Ca}$ between both morphotypes (Gray et al., 2018). We acknowledge that the scarcity of G. ruber sl in our samples means that our data set for this morphotype is too small to draw firm conclusions, and this warrants further investigation at other study sites. Nonetheless, the closeness of the morphotypes in terms of δ¹¹B and depth habitat throughout our record implies any inadvertent sampling of G. ruber sl in the G. ruber ss fraction in this study and location would not significantly bias our reconstructions.

4.2.2 Change in upwelling and CO₂ disequilibrium

ODP sites 871 and 999 are today both located in stratified oligotrophic environments with a deep modern thermocline (the base of the thermocline is at ∼200 and 400 m at ODP 871 and 999, respectively; Olsen et al., 2016). It should be noted, however, that both sites are situated relatively close to regions displaying ΔpCO₂>40 ppm (Fig 1). However, if local upwelling occurred over the study interval or if these areas of upwelled water expanded, we would expect these periods to be characterised by relatively low SST, high surface δ¹⁸O, and low surface δ¹³C due to an increased influence of colder and more remineralised deep waters. The identified anomalous intervals in residual CO₂ at ODP 871 (e.g. at ∼210, ∼290 ka, Fig. 5) show no particular anomaly in planktonic C and O isotopes (Fig. S5) or in SST (Figs. 2 and S8), ruling out significant variations in upwelling at that site. The $\mathrm{Mg}/\mathrm{Ca}$-derived SST record of nearby Site MD97-2140 (Fig. S8) from the Western Pacific Warm Pool (de Garidel-Thoron et al., 2005), a location outside of the upwelling from the Pacific cold tongue, confirms this view in that the periods of high CO₂ offset at Site 871 are not associated with relatively cold periods at site MD97-2140. Equally, no SST anomaly was identified at ODP 999 to be coincident with the intervals of high residual CO₂. Foster and Sexton (2014) have also reconstructed CO₂ zonally across the equatorial Atlantic and the Caribbean and showed that while enhanced disequilibrium was detected in the eastern Atlantic, Site 999 has remained in equilibrium with the atmosphere for the last 30 kyr at least. This suggests the CO₂ anomalies revealed in Fig. 5 are not the result of enhanced local disequilibrium via sub-surface water mixing. Whilst SST is a first-order constraint on upwelling, we acknowledge future constraints are needed using paired proxies of local CO₂, temperature, and productivity to evaluate changes in local CO₂ fluxes.

Figure 5 (a, b) CO₂ offset (defined as offset being equal to CO_{2_δ_11B-derived}–CO_{2_ice}) for ODP Sites 999 (this study and Chalk et al., 2017) and 871. See the text for error bar calculations. (c, d) Fragmentation index at Site 999 (Schmidt et al., 2006) and 871 (this study). Red dots in (c, d) are the fragments above the upper quartile (and corresponding CO₂ in a, b, red dots). Green dots represent periods of low fragments below the upper quartile (and corresponding CO₂ in a, b, green dots).

4.2.3 Partial dissolution

The CO₂ derived from G. ruber δ¹¹B at ODP 999 and 871 appears to show, at the first order at least, positive CO₂ offset during periods of high fragmentation (∼100, ∼210, ∼400 ka, filled red circles in Fig. 5, defined by the upper 25 % quantile of fragments) following a “Pacific style” dissolution cycle (better preservation and lower fragmentation during glacial periods). Periods of high fragmentation at ODP Site 999 and 871 correspond to a positive CO₂ offset 65 % and 75 % of the time, respectively, and 35 % and 25 % of the time to a negative or no (i.e. ±10 ppm) CO₂ offset (note that values ± 10 ppm were omitted in the criteria for positive or negative offset). We also note that almost all CO₂ offset uncertainties (2σ) overlap with the 0 line; hence, the percentage of CO₂ offsets that are above or below the 0 line should be interpreted with caution.

In detail, however, a cross-plot of fragment counts and CO₂ offset (Fig. S9) fitted with a linear regression shows no significant correlation for both core site 999 (r²=0.06, p=0.03) and 871 (r²=0.002, p=0.77). Although it should be noted that this simple linear regression presupposes a linear relationship between the variables and does not account for the significant uncertainty in both CO₂ offset and fragmentation index. In particular, the CO₂ offset carries the uncertainty from the interpolated ice core CO₂ (see Sect. 2). Fragment counts at ODP 999 also come with the additional uncertainty related to the interpolation of the record of Schmidt et al. (2006), whereas fragments counts and δ¹¹B-derived CO₂ at 871 are measured on the same samples. A cross-correlation function also shows no correlation between CO₂ offset and fragmentation (Fig. S10).

While it seems unlikely the small offsets observed are fully explained by partial dissolution, the positive CO₂ offsets observed during some periods of high fragmentation index (Fig. 5) are in line with trends observed in other species like T. sacculifer (sacc). For instance, field studies observed lower δ¹¹B in T. sacculifer for core-top samples from deeper ocean sites bathed by waters with a low calcite saturation state (Hönisch and Hemming, 2004; Seki et al., 2010). Tests of T. sacculifer can contain a significant proportion of gametogenic calcite (ranging 30 % to 75 % of the weight of pregametogenic calcite, Bé, 1980; Caron et al., 1990), which forms at the end of the life cycle in deeper lower pH cold waters. It has been suggested that δ¹¹B is lower in gametogenic calcite than in the primary test (Ni et al., 2007), reflecting the digestion and expulsion of symbionts (Bé et al., 1983) before gametogenesis and driving a relative acidification of the micro-environment (no CO₂ uptake by photosynthesis) around the foraminifera (Zeebe et al., 2003; Hönisch et al., 2003; Henehan et al., 2016) and movement to deeper more acidic waters during that life stage. It has been shown that this gametogenic calcite is more resistant to dissolution (Hemleben et al., 1989; Wycech et al., 2018), resulting in partial dissolution acting preferentially on ontogenic calcite and driving δ¹¹B in the residual test to lower isotopic composition.

However, while the decrease in δ¹¹B in tests of T. sacculifer found in corrosive waters is well explained by the lighter isotopic composition of gametogenic calcite, G. ruber tests do not contain such gametogenic calcite (Caron et al., 1990). Hence, if the observed occasional decrease in δ¹¹B (low pH, high CO₂) was caused by partial dissolution, it needs to be explained by other processes. It should also be considered that alternative measures and proxies of dissolution (e.g. benthic $\mathrm{B}/\mathrm{Ca}$ as an indicator of bottom water carbonate ion concentration) may yield more quantitative constraints on the importance of dissolution in generating our observed CO₂ offsets. Some studies have shown that laboratory-dissolved specimens of T. sacculifer (Sadekov et al., 2010) and naturally dissolved specimens of G. ruber (Iwasaki et al., 2019) undergo targeted partial preferential dissolution of the shell. However, variations in intra-shell δ¹¹B are currently unknown due to limitations in laser ablation techniques that currently impede a direct evaluation of δ¹¹B heterogeneity in foraminifera chambers. Future studies are needed to constrain the δ¹¹B spatial distribution in foraminiferal shells caused by potential variations in δ¹¹B from dissolution, ontogeny (e.g. Meilland et al., 2021), and/or vital effects (e.g. change in photosymbiotic activity throughout the life cycle, Lombard et al., 2009; Henehan et al., 2013; Takagi et al., 2019). In the absence of these constraints, we conclude that partial dissolution is unlikely to be a significant driver of the δ¹¹B–CO₂ records we present here. Even though it was thought to be a species susceptible to dissolution (Berger, 1970), we confirm that the δ¹¹B of G. ruber appears more resistant to dissolution-driven modification than T. sacculifer.

4.2.4 Effect of dissolution on $\mathrm{Mg}/\mathrm{Ca}$ and calculated CO₂

The direction of change in $\mathrm{Mg}/\mathrm{Ca}$ with partial dissolution is towards lower ratios in partially dissolved foraminifera (e.g. Brown and Elderfield, 1996; Dekens et al., 2002; Fehrenbacher and Martin, 2014). If the $\mathrm{Mg}/\mathrm{Ca}$ is impacted during periods of high fragmentation, the lower ratio would result in lower temperatures, leading to lower calculated CO₂ values (Eq. 7). This effect is opposite to the occasional positive deviation of CO₂ observed during intervals of high fragmentation at ODP Site 999. While the weak correlation between fragmentation and CO₂ precludes a firm interpretation of dissolution effect, we conclude that the effect of partial dissolution on $\mathrm{Mg}/\mathrm{Ca}$ ratio and resulting CO₂ (if any) are negligible and not responsible for the CO₂ offsets observed during intervals of high fragmentation.

4.2.5 Change in the second carbonate parameter: alkalinity

Past changes in TA are poorly constrained, although some constraints are starting to emerge for the late Quaternary (e.g. Cartapanis et al., 2018). However, since pH is directly determined by δ¹¹B, pH defines the ratio of alkalinity to DIC (see Sect. S11). Hence, at any given pH, any change in alkalinity must be counteracted by a change in DIC, which has the opposing effect on CO₂. This is demonstrated by the tight relation between pH and CO₂ highlighted by our data (Fig. 4). The largest residual CO₂ is ∼50 ppm at ODP 999. To produce an effective alkalinity-driven change in CO₂ of this magnitude at a given pH requires an alkalinity reduction of about ∼300 to 500 µmol mol⁻¹ (Fig. S12). This is far larger than any expected change over a glacial cycle (Cartapanis et al., 2018; Hönisch et al., 2009). We therefore rule out varying TA as the cause of the minor CO₂ offsets observed (Fig. 5).

4.2.6 Improving the δ¹¹B-pH G. ruber calibration

A further potential cause for the minor offsets observed between δ¹¹B-derived and ice core CO₂ could be a small inaccuracy in the calibration between δ¹¹B of foraminifera and borate for G. ruber (Henehan et al., 2013). Having the ice core data to compare with δ¹¹B-derived CO₂ offers an opportunity to explore the effect of altering the input variables of the pH-CO₂ calculation to see if doing so improves the fit to ice core values. Note that such an exercise is for illustrative purposes only because we seek to retain the independence offered by the δ¹¹B-calibrated data in the context of CO₂ forcing (Sect. 4.3). Nonetheless, in future work we suggest this calibration can be applied in tandem to the empirical relationship of Henehan et al. (2013). The published (Henehan et al., 2013) and obtained optimised calibration (Fig. S13) are as follows:

$$\begin{array}{}\text{(11)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{borate}}=& {\displaystyle \frac{{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{foram}}-\mathrm{8.87}(\pm \mathrm{1.52})}{\mathrm{0.60}(\pm \mathrm{0.09})}}\\ & (\mathrm{Henehan}\mathrm{et}\mathrm{al}.,\mathrm{2013}),\end{array}\text{(12)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{borate}}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{foram}}-\mathrm{6.49}}{\mathrm{0.71}}}\left(\mathrm{optimised}\mathrm{calibration}\right).\end{array}$$

The newly calculated CO₂ with the updated calibration shows an improved average CO₂ offset (Fig. 6) of $-\mathrm{4}\pm \mathrm{36}\left(\mathrm{2}\mathit{\sigma }\right)$ ppm (vs. 13±46 (2σ) ppm with the calibration of Henehan et al., 2013) and an RMSE of 18 ppm (vs. 26 ppm with the published calibration).

Figure 6 Composite δ¹¹B-derived CO₂ from both core sites 999 and 871 using the published δ¹¹B_{foram–borate} calibration (grey points, Henehan et al., 2013) and the optimised calibration (red points). The black line is the Antarctic composite ice core CO₂ record (Bereiter et al., 2015).

When analysing the CO₂ offset using the optimised G. ruber calibration and the fragmentation index at each core location (the same approach as Fig. 5), we observe that intervals of high fragments (defined as values above the upper quartile) are no longer preferentially associated with positive CO₂ offset (Fig. S14). Intervals of high fragments occur 5 % and 33 % of the time at Site 999 and Site 871, respectively, during positive CO₂ offsets (and 95 % and 67 % of the time during negative or no offset to the ice cores, respectively).

This analysis shows that a small change in the borate G. ruber δ¹¹B calibration is enough to improve the fit to the ice core and diminishes the apparent correlation between high fragmentation and CO₂ offset (Fig. S14) and that uncertainty in the δ¹¹B_foram-borate calibration of Henehan et al. (2013) can – at least partly – explain the minor discrepancies we observe between δ¹¹B-derived and ice core CO₂.

4.3 Relative CO₂ forcing and pH

Our new pH data, added to the existing compilation, show a good agreement with the formalism defined by Hain et al. (2018; Fig. 4). It should be noted that CO₂ in this case is provided by the ice core directly and is not estimated from the δ¹¹B-derived pH. As discussed above, because these two proxies are independent of one another, the slope of their relationship may be used to interrogate the mechanisms of CO₂ change. Our data fall between the CaCO3 (plain yellow line) and the DIC (dotted-dashed blue line) end-members, suggesting that the CO₂ change observed on glacial–interglacial timescales was driven by a mix of mechanisms rather than a single cause. This is in line with studies that require a number of mechanisms to explain glacial–interglacial CO₂ change such as the soft tissue pump, carbonate compensation pump, solubility pump (e.g. Brovkin et al., 2007; Kohfeld and Ridgwell, 2009; Hain et al., 2010; Chalk et al., 2019; Sigman et al., 2021), and disequilibrium pump (Eggleston and Galbraith, 2018). We note that this is a preliminary interpretation because of the sensitivity of our finding to pH uncertainty (Sect. 3.4, Fig. S7). To overcome this ambiguity in estimating past ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ and to better deconvolve the driving mechanisms of glacial–interglacial CO₂ change, we recommend that future studies collect pH data at higher temporal resolution to examine the change in slope through a glacial cycle and strive to further quantify and reduce uncertainties related to pH determination.

The close agreement of the pH and ice core CO₂ data with the theoretical relationships has a number of consequences for the reconstruction of CO₂ change during periods of Earth's history beyond the ice core CO₂ and climate records where constraints on δ¹¹B_sw and the second carbonate parameter and temperature are uncertain. The ΔpH formalism still requires an estimation of δ¹¹B_sw and temperature (for the pK_B term, Eq. 5); however, as discussed in Hain et al. (2018), while absolute reconstruction of pH is significantly influenced by estimates of δ¹¹B_sw and temperature, reconstruction of relative pH change (ΔpH) is inherently much less sensitive to these input variables.

Reconstructing ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ from ΔpH is ideally applicable only on relatively short timescales less than 1 Myr, when δ¹¹B_sw is likely to be invariant given the multi-million-year residence time of boron in the ocean (Lemarchand et al., 2002; Greenop et al., 2017). Furthermore, to reconstruct ΔF${}_{{\mathrm{CO}}_{\mathrm{2}}}$ (and thus climate sensitivity to CO₂), the formalism can be applied as long as ΔpH remains the overwhelming control in Eq. (2). This is dependent on the residence time of carbon in the ocean with respect to silicate weathering – approximately a million years (Hain et al., 2018), such that net carbon addition to or removal from the Earth system through volcanic outgassing or silicate weathering is likely to be minor over the million-year timescale. However, during some short events, for instance the Paleocene–Eocene Thermal Maximum, considerable carbon was added to the system in <200 kyr (e.g. Gutjahr et al., 2017), invalidating the formulation described in Eq. (2) on these intervals. We also emphasise that this formalism is only valid as long as core sites remain in equilibrium with the atmosphere.

4.4 Caveats and future studies

The aim of this study is to evaluate the capacity of the δ¹¹B–pH proxy in G. ruber to accurately reconstruct atmospheric CO₂ in the past. The overall agreement with the high-confidence ice core CO₂ (e.g. Bereiter et al., 2015) is very promising and gives confidence to δ¹¹B-derived CO₂ reconstructions beyond the ice core record (>800 ka). We have, however, identified occasional minor offsets between the two records and explored potential drivers (partial dissolution, δ¹¹B_{foram–borate} calibration, local air–sea disequilibrium). It is likely that the minor disagreement observed (Fig. 5) has a combination of drivers and that a single mechanism is not solely responsible for the CO₂ offsets observed. To confirm these trends, we recommend future work to focus on the following points.

1. The improved δ¹¹B calibration approach should be tested at more core locations. We note that the improved calibration to the ice core records reported here was achieved using data from two sites. While care is taken in the choice of study site to minimise air–sea CO₂ disequilibrium and sediment dissolution, the newly defined improved δ¹¹B_{foram–borate} calibration should be seen as an exercise that is tailored to the available data in this study, and future high-resolution studies can apply the method used here (Sect. 4.2.6) to further test how the G. ruber calibration changes if CO₂ offsets occur in a similar fashion (i.e. at a particular time in each glacial cycles). We note the importance of high-resolution (at least 3 kyr) sampling in future studies because most CO₂ offsets observed are short lived.

2. A multiproxy approach is ideally needed. In particular, reliable indicators of temperature and productivity are needed to assess change in upwelling and foraminifera ecology. We encourage future studies to expand high-resolution boron-derived CO₂ record and ancillary data (C and O isotopes, proxy of carbonate preservation and bottom water corrosiveness, biological productivity) to further constrain the capacity of the boron isotope pH/CO₂ proxy to generate reliable CO₂ records. As more recent IODP expeditions include porewater data, constraints on bottom water conditions and the degree of corrosiveness at a given site will become available to evaluate the impact on δ¹¹B signals in foraminifera.

3. Efforts should continue to decrease the analytical uncertainty associated with a δ¹¹B measurement by MC-ICPMS because this still accounts for ∼40 % of the total uncertainty associated with each δ¹¹B-derived CO₂ estimate.

4. We find little evidence to suggest that partial dissolution of foraminiferal tests (G. ruber) is a major driver of uncertainty in δ¹¹B-derived CO₂ estimates, but well-constrained dissolution experiments are desirable because of site-to-site differences in foraminifera taphonomy.

5 Conclusion

We carried out the most thorough test to date of the δ¹¹B–pH (CO₂) proxy by comparing new high-resolution (3 to 6 kyr per sample) boron-isotope-based pH and CO₂ at two locations with CO₂ from the ice core record. Results suggest that the boron isotope proxy is robust and suited to reconstructing CO₂ to a precision of ±46 ppm (2σ, RMSE = 26 ppm) over this interval, with little or no systematic bias shown by a mean residual of 13±46 (2σ) ppm. This provides high confidence to the application of the proxy beyond the reach of the ice core records.

Despite the overall good agreement, there are some minor short-lived CO₂ offsets that appear to have some temporal structure and we explored a number of possible drivers. A visual correlation between CO₂ offset and fragmentation index at core site 999 is observed (Fig. 5) but is not statistically significant. The effect of partial dissolution on δ¹¹B in G. ruber appears to be negligible in our record, but the possible heterogeneity of δ¹¹B within shells and variable susceptibly to dissolution of the different parts of the foraminifera encourage further exploration.

A revised δ¹¹B borate–foram calibration was calculated by minimising the offset between δ¹¹B-derived CO₂ and ice core CO₂ using published calibrations (Henehan et al., 2013). While the new calibration improves the fit to the ice core records, we caution against its use to estimate CO₂ given that it is no longer independent of the ice core or the assumptions we make here to calculate CO₂ (i.e. that TA is constant).

The formalism established by Hain et al. (2018) is robust, showing that relative CO₂ forcing in the past can be determined from pH change alone, even in the face of significant uncertainty in δ¹¹B of seawater and without the need to determine a second carbonate parameter. This will not only be of great interest to determine CO₂ forcing in ancient geological times where δ¹¹B of seawater and a second carbonate parameter are poorly constrained, but the nature of the observed relationship over the last 400 kyr confirms that multiple drivers are likely responsible for glacial–interglacial CO₂ change.