No detectable influence of the carbonate ion effect on changes in stable carbon isotope ratios (δ13C) of shallow dwelling planktic foraminifera over the past 160 kyr

Laboratory experiments showed that the isotopic fractionation of δ13 C and of δ18 O during calcite formation of planktic foraminifera are species-specific functions of ambient CO 32- concentration. This effect became known as the carbonate ion effect (CIE), whose role for the interpretation of marine sediment data will be investigated here in an in-depth analysis of the 13 C cycle. For this investigation, we constructed new 160 kyr long mono-specific stacks of changes in both δ13 C and δ18 O from either the planktic foraminifera Globigerinoides ruber (rub) or Trilobatus sacculifer (sac) from 112 and 40 marine records, respectively, from the wider tropics (latitudes below 38°). Both mono-specific time series Δ(δ13Crub) and Δ(δ13Csac) are very similar to each other, and a linear regression through a scatter plot of both data sets has a slope of ∼ 0.99 – although the laboratory-based CIE for both species differs by a factor of nearly 2, implying that they should record distinctly different changes in δ13C , if we accept that the carbonate ion concentration changes on glacial–interglacial timescales. For a deeper understanding of the 13 C cycle, we use the Solid Earth version of the Box model of the Isotopic Carbon cYCLE (BICYLE-SE) to calculate how surface-ocean CO 32- should have varied over time in order to be able to calculate the potential offsets which would by caused by the CIE quantified in culture experiments. Our simulations are forced with atmospheric reconstructions of CO 2 and δ13 CO 2 derived from ice cores to obtain a carbon cycle which should at least at the surface ocean be as close as possible to expected conditions and which in the deep ocean largely agrees with the carbon isotope ratio of dissolved inorganic carbon (DIC), δ13CDIC , as reconstructed from benthic foraminifera. We find that both Δ(δ13Crub) and Δ(δ13Csac) agree better with changes in simulated δ13CDIC when ignoring the CIE than those time series which were corrected for the CIE. The combination of data- and model-based evidence for the lack of a role for the CIE in Δ(δ13Crub) and Δ(δ13Csac) suggests that the CIE as measured in laboratory experiments is not directly transferable to the interpretation of marine sediment records. The much smaller CIE-to-glacial–interglacial-signal ratio in foraminifera δ18 O, when compared to δ13 C, prevents us from drawing robust conclusions on the role of the CIE in δ18 O as recorded in the hard shells of both species. However, theories propose that the CIE in both δ13 C and δ18 O depends on the pH in the surrounding water, suggesting that the CIE should be detectable in neither or both of the isotopes. Whether this lack of role of the CIE in the interpretation of planktic paleo-data is a general feature or is restricted to the two species investigated here needs to be checked with further data from other planktic foraminiferal species.

For a reconstruction of past changes in the ocean and in the carbon cycle, various variables are measured in microfossils obtained from marine sediment cores. Among the most widely used are the stable carbon and oxygen isotope ratios, ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathit{\delta }}^{\mathrm{18}}$ O, from hard shells of planktic and benthic foraminifera. Since the publication of the first stable isotope time series [42] a vast number of stable isotope records have been published and to a large part compiled in the World Atlas of late Quaternary Foraminiferal Oxygen and Carbon Isotope Ratios [119]. One of the fundamental problems with the interpretation of foraminiferal isotope ratios is how and why a stable isotope signal was altered on its way from the seawater to the shell of living foraminifera. Are there vital and other effects necessary to be considered when interpreting the paleo-records [14], [197], [67]?

The carbonate ion effect (CIE) is one of these potentially important effects that might alter the isotopic signal. The CIE implies that both ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathit{\delta }}^{\mathrm{18}}$ O measured in hard shells of marine organisms undergo isotopic fractionation during calcite formation with the amplitude of the fractionation, among other factors, being a function of the carbonate ion concentration ([CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ]) of the surrounding seawater [168]. The CIE has been found to be species-specific [169], ranging from $-\mathrm{4.7}\times {\mathrm{10}}^{-\mathrm{3}}$  ‰ to $-\mathrm{13.0}\times {\mathrm{10}}^{-\mathrm{3}}$  ‰ per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ of [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] for ${\mathit{\delta }}^{\mathrm{13}}$ C and between $-\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{3}}$  ‰ and $-\mathrm{4.5}\times {\mathrm{10}}^{-\mathrm{3}}$  ‰ per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ of [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] for ${\mathit{\delta }}^{\mathrm{18}}$ O in four planktic foraminifera. The CIE for ${\mathit{\delta }}^{\mathrm{13}}$ C has been explained for Orbulina universa, a spinose, symbiont-bearing species, by the pH-related distribution of dissolved inorganic carbon (DIC) into its three species: CO ${}_{\mathrm{2}}$ , CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ and HCO ${}_{\mathrm{3}}^{-}$ [[190], [196]]. The CIE on ${\mathit{\delta }}^{\mathrm{18}}$ O is also explained by the CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ -related varying pH [194]. These theories, however, were unable to base the full amplitudes found in experiments solely on this effect. The CIE is maybe the most prominent isotopic fractionation effect which has to be considered when interpreting the paleo-records, but others, e.g vital effects and dependency on light, temperature, pressure and shell size, have been put forward [165], [166], [167], [163], [164], [125]. The CIE is found to play a minor role when comparing late Holocene deep-ocean ${\mathit{\delta }}^{\mathrm{13}}$ C in benthic foraminifera, with ${\mathit{\delta }}^{\mathrm{13}}$ C of DIC ( ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ ) [152] being responsible for $-\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{3}}$  ‰ per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ of [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] disturbance in the recorded signal. In a recent study focusing on the benthic species Cibicidoides wuellerstorfi, $-\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{3}}$  ‰ per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ of [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] has been obtained for the late Holocene [123]. Both studies also found in addition to the CIE that ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{benthic}}$ was also partly controlled by other variables, mainly pressure (water depth) and temperature.

The CIE in planktic foraminifera is one of the reasons why the interpretation of the whole ${\mathit{\delta }}^{\mathrm{13}}$ C cycle over glacial–interglacial timescales is still challenging. In a compilation of foraminiferal ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurements covering the past 150 kyr, [124] find relatively large disagreements between different planktic ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ records within a region compared to benthic records, consistent with large uncertainty attributed to the estimation of ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ from planktic species. Since benthic compilations are less affected by the CIE, they should, however, robustly constrain deep-ocean changes in ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ . A more recent compilation of benthic ${\mathit{\delta }}^{\mathrm{13}}$ C was given in [89]. Furthermore, ${\mathit{\delta }}^{\mathrm{13}}$ C of atmospheric CO ${}_{\mathrm{2}}$ ( ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ ) is now available over the last 155 kyr [39] from ice cores. So far, tight constraints on the change in surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ are missing from our understanding, but, in principle, this information should be recorded in the hard shells of planktic foraminifera, even if hidden under the CIE.

In this study, we therefore aim to construct a robust time series of orbital changes in surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ based on planktic foraminifera data. We compiled ${\mathit{\delta }}^{\mathrm{13}}$ C data largely based on [119] covering the last 160 kyr. In order to be able to apply any species-specific CIE corrections, we compile mono-specific isotope records on the widely abundant shallow-dwelling planktic foraminifera species Globigerinoides ruber (G. ruber or rub) and Trilobatus sacculifer (T. sacculifer or sac) into stacks. Due to their spatial distribution [44], this species selection leads effectively to the construction of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ stacks based on sediment core data from latitudes smaller than 40°, potentially informing us about mean changes of ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ on orbital timescales in the surface of the wider tropical ocean. Accompanied stacks of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}\right)$ from the same cores will add further information on the CIE in ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ .

A first surface-ocean ${\mathit{\delta }}^{\mathrm{13}}$ C stack based on data from T. sacculifer obtained from five equatorial Atlantic records has been constructed by [25] without any knowledge on the CIE. Furthermore, [169] used data from G. ruber and T. sacculifer from a single core in the Indian Ocean and the lab-based size of their species-specific CIE to deconvolve surface-ocean [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ]. Here we will use our new mono-specific ${\mathit{\delta }}^{\mathrm{13}}$ C stacks, which, due to the underlying number of records, have a much higher signal-to-noise ratio to test the robustness of their findings.

In the following, we will investigate the connection of ${\mathit{\delta }}^{\mathrm{13}}$ C in the atmosphere and ocean in closer detail in order to improve our understanding of the ${}^{\mathrm{13}}$ C cycle. The remainder of the article is structured as follows. Firstly (Sect. ), we describe the construction of our mono-specific ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ anomaly stacks, $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ (and of the accompanied ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ anomalies). Some published benthic ${\mathit{\delta }}^{\mathrm{13}}$ C data are also needed for our understanding (Sect. ). For a deeper interpretation, the global isotope-enabled carbon cycle model, the Solid Earth version of the Box model of the Isotopic Carbon cYCLE (BICYCLE-SE) [74], which has been proven to simulate glacial–interglacial (G-IG) in the carbon cycle reasonably well, is used. The model is briefly described in Sect. , including a completely revised parametrisation of the ${}^{\mathrm{13}}$ C cycle. In Sect. , we discuss what we already know from data on the ${\mathit{\delta }}^{\mathrm{13}}$ C cycle and the role the CIE might play. We then analyse the simulated ${\mathit{\delta }}^{\mathrm{13}}$ C cycle in our model results in Sect. . This enables us to evaluate (Sect. ) if our stacks, $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ , are good representations of changes in ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the wider tropical surface ocean or if corrections such as the CIE need to be applied. Finally, we briefly discuss the CIE in ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}$ and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}$ (Sect. ) before we come to our conclusions (Sect. ).

To construct time series of low-latitude ${\mathit{\delta }}^{\mathrm{13}}$ C variations through the past 160 kyr, we selected 112 and 40 ${\mathit{\delta }}^{\mathrm{13}}$ C records of the shallow-dwelling planktic foraminifera G. ruber and T. sacculifer, respectively, predominantly from the World Atlas of late Quaternary Foraminiferal Oxygen and Carbon Isotope Ratios [119]. A list of the isotope records contributing to our stacks with relevant metadata, references to the original publications and data sources is provided in Table S1. All data citations are summarized in a condensed way in Appendix A. In three sediment cores, time series from both G. ruber (white) and G. ruber (pink) contribute to our G. ruber stacks, while data from 22 cores contain mono-specific data from both G. ruber and T. sacculifer. All combined, our data selection is based on material from 127 sediment cores. The core sites cover a latitudinal range from 37.6° N to 36.7° S for G. ruber and from 32.8° N to 31.3° S for T. sacculifer in all major ocean basins (Fig. ), although the contributions from individual cores (and therefore the latitudinal range) changed over time (Fig. c). Our age models are based on either radiocarbon ages or oxygen isotope stratigraphy or a combination of both methods. To calibrate radiocarbon ages, we first subtracted a simulated local reservoir age from the nearest grid box of the modelling experiments conducted for Marine20 [[21], [54]] and then calibrated the corrected radiocarbon age with the IntCal20 calibration curve [138]. For core sections with insufficient radiocarbon coverage or outside the radiocarbon dating range, ages were added through the visual alignment with the isotope stacks by [90] and [91] using the software PaleoDataView [82]. In a few cases, age models were derived by visual alignment with the oxygen isotope records of well-dated nearby cores. The details of the age model construction are available in the netCDF files of the age models in the corresponding PaleoDataView collection [73]. A continuous age model was then constructed with the age modelling software BACON [15]. For each record, we produced an ensemble of 1000 time series by combining 1000 BACON-generated age models with 1000 down-core ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathit{\delta }}^{\mathrm{18}}$ O series by adding a random value within the typical analytical $\mathrm{1}\mathit{\sigma }$ -uncertainty of 0.05 % and 0.07 ‰ to each down-core ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathit{\delta }}^{\mathrm{18}}$ O value, respectively. The resulting 1000 ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathit{\delta }}^{\mathrm{18}}$ O time series were then interpolated to a time step of 1 kyr to calculate the mean and the standard deviation of the time series ensembles. The averaging of the individual ensemble members led to a considerable smoothing of the final time series.

Graph: Figure 1 Location of the 127 sediment cores from which data have been compiled for this study. Data from the planktic species G. ruber have been included in 87 cores, and in 18 cores data from T. sacculifer have been included, while 22 cores provide mono-specific data from both species.

Graph: Figure 2 Stacks of anomalies in (a, e)δ18O and (b, f) δ13C from the planktic species G. ruber and T. sacculifer across the last 160 kyr. Mean anomalies (±1 SE) are calculated with respect to the mean of 21–19 kyr BP (vertical blue band). Details of underlying data are found in Table S1. (c, g) Latitudinal distribution of cores contributing to the stack (mean and full range) and (d, h) stack count. Either data from all cores for each species are compiled (a–d) or (e–h) they are compiled from a reduced core selection, in which contributing cores cover both Termination 1 and 2 (T1 + T2).

Although the size class used for stable isotope measurements can vary considerably among records, it is common practice to use a fairly constant shell size down-core to minimise size-related effects on both oxygen and carbon isotope ratios [125]. To provide a common baseline, we corrected all isotope records by their individual mean values for the period from 21 to 19 kyr BP marked as Last Glacial Maximum (LGM) in various plots. To produce final isotope stacks, we averaged all corrected time series and calculated the standard error (SE) of the means at 1 kyr intervals. The final mono-specific stacks of both ${\mathit{\delta }}^{\mathrm{18}}$ O and ${\mathit{\delta }}^{\mathrm{13}}$ C anomalies based on G. ruber and T. sacculifer are plotted in Fig. a, b. The oxygen isotope stacks are also shown here to give a clear reference for G-IG changes: ${\mathit{\delta }}^{\mathrm{18}}$ O has its maxima during peak glacial times and its minima during peak interglacials. In Sect.  we will come back to these data to discuss the CIE in ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}$ and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}$ . To test the extent to which the data distribution affects the stacks, we generated two versions of the stacks: one based on all records (Fig. a–d) and an alternative based only on records which contain both Terminations (T1  $+$  T2; Fig. e–h). The stack counts (Fig. d, h) show that the two versions differ mainly in the younger half and that they are identical beyond 85 kyr BP. The latitudinal ranges in the younger half are slightly smaller for the compilations T1  $+$  T2 than when all cores are compiled, but the mean latitudes of all cores are throughout the covered time window of the last 160 kyr in all cases (for both species and for both compilations) close to the Equator (Fig. d, g). This consistency in the mean latitude suggests that the incoming light which varied in its annual mean values between $\sim$  420 W m ${}^{-\mathrm{2}}$ at the Equator and $\sim$  330 W m ${}^{-\mathrm{2}}$ around latitudes of 40° [83] should only marginally affect the isotopic fractionation [167].

The focus of this study is the ${\mathit{\delta }}^{\mathrm{13}}$ C of the surface ocean. However, for a rough comparison of ${\mathit{\delta }}^{\mathrm{13}}$ C changes in the deep ocean, we rely on the published ${\mathit{\delta }}^{\mathrm{13}}$ C stack constructed from six deep Pacific cores as contained in [89]. The six cores are all ODP cores (677, 846, 849, 1123, 1143 and 1208) from between 2700 and 3500 m water depth, located between 42° S and 36° N. The deep Pacific ${\mathit{\delta }}^{\mathrm{13}}$ C stack should cover the most-depleted end member of the marine ${\mathit{\delta }}^{\mathrm{13}}$ C cycle (Fig. d) and should give some indication how ${\mathit{\delta }}^{\mathrm{13}}$ C in the deep ocean is performing in our simulations. More details on the stack are found in [89].

Graph: Figure 3 Carbon cycle time series of the last 160 kyr, including the Penultimate Glacial Maximum (PGM) and the Last Glacial Maximum (LGM) and Terminations 1 and 2 (T1 and T2). Spline of atmospheric CO2 (a) and δ13CO2 (b) based on data from various ice cores (grey, ±1σ around the mean; [[77], [39]]) and highly resolved recent data from the "horizontal ice core" approach in Taylor Glacier (yellow; [[11], [104]]). (c) Δ(δ13Crub) and Δ(δ13Csac) averaging signals in the wider tropical surface ocean [119]. (d) Deep-ocean δ13C from benthic foraminifera stacked from six Pacific cores [89].

At the core of BICYCLE – the Box model of the Isotopic Carbon cYCLE – sits an ocean (O) with 10 boxes and a terrestrial biosphere consisting of 7 boxes (B) together with a 1-box atmosphere (A), in which the concentration of carbon (as DIC in the ocean, as $p$ CO ${}_{\mathrm{2}}$ in the atmosphere and as organic carbon in the biosphere) and both of the isotopes ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathrm{\Delta }}^{\mathrm{14}}$ C are traced [75]. Furthermore, in the ocean alkalinity, PO ${}_{\mathrm{4}}^{\mathrm{3}-}$ as a macro-nutrient and O ${}_{\mathrm{2}}$ are represented. From the two variables in the marine carbonate system (DIC and alkalinity), all other variables (CO ${}_{\mathrm{2}}$ , HCO ${}_{\mathrm{3}}^{-}$ , CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ and pH) are calculated according to [195] with updates of the dissociation constants pK ${}_{\mathrm{1}}$ and pK ${}_{\mathrm{2}}$ [112]. The 10 ocean boxes distinguish 100 m deep equatorial (or wider tropical) surface waters in the Atlantic and Indo-Pacific from 1000 m deep surface-ocean boxes in the high latitudes (North Atlantic, Southern Ocean and North Pacific). In the model, wider tropical boxes range from 40° S to 40° N in the Indo-Pacific and to 50° N in the Atlantic, rather similar to the latitudinal coverage of the sediment cores from which $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ have been constructed. Deep-ocean boxes represent all waters below 1 km in the three basins: Atlantic, Southern Ocean and Indo-Pacific. In the equatorial regions the waters between 100 and 1000 m depth are described by intermediate boxes. The terrestrial biosphere [72] distinguishes C ${}_{\mathrm{3}}$ and C ${}_{\mathrm{4}}$ photosynthesis of grasses and trees and soil carbon with different turnover times of up to 1000 years.

The model extension towards the version BICYCLE-SE used here, which can take care of solid Earth processes, is sketched in Fig. . The main improvement documented in detail in [74] is the implementation of a sediment module, which captures early diagenesis in an 8 cm deep sedimentary mixed layer (M), under which numerous historical layers are implemented. In effect, we now simulate the subsystem of the global carbon cycle consisting of atmosphere, ocean, terrestrial biosphere and sedimentary mixed layer (AOBM) within BICYCLE-SE. In each of the three ocean basins (Atlantic, Southern Ocean and Indo-Pacific) the pressure-dependent carbonate system is calculated for every 100 m water depth, and, depending on the over- or undersaturation of the carbonate ion concentration, CaCO ${}_{\mathrm{3}}$ is either accumulated or dissolved. The parametrisation and realisation of the sedimentary processes directly follow [121] and [120]. The carbon isotopes in the sedimentary mixed layer are only followed in aggregated boxes (one for each of the three ocean basins).

Graph: Figure 4 Sketch of the Box model of the isotopic carbon cycle, version solid Earth (BICYCLE-SE), modified from [74]. V is the outgassing of CO2 from volcanoes on land potentially and temporally overlain by land ice and from hot-spot island volcanoes (and mid-ocean ridges, not shown) influenced by the changing sea level. C is shallow-water carbonate deposition due to coral reef growth. Si-W is silicate weathering, and Ca-W is carbonate weathering, with different sources of C but both delivering HCO3--ions into the ocean. P is the PO43- riverine input and sedimentary burial, and S is the CaCO3 sedimentation and dissolution. A–B is the atmosphere–biosphere exchange of CO2, and A–O is the atmosphere–ocean exchange of CO2. The cyan-coloured broken circles mimic the two overturning cells in the Atlantic and Indo-Pacific oceans. The isotopic fractionation ε during exchange processes, or the prescribed δ13C of external fluxes, is given, summarising the parametrisation of the 13C cycle within the model.

Equipping BICYCLE with a process-based sediment module enables the revised model version, BICYCLE-SE, to address questions related to changes in solid Earth carbon fluxes in detail and in the long term. Roughly speaking, the following processes are considered. CO ${}_{\mathrm{2}}$ outgassing from volcanoes on land, hot-spot island volcanoes and mid-ocean ridge (MOR) hydrothermal activity is realised as partly being dependent on the changing sea level. Coral reef growth is a known shallow-water carbonate sink, which to some extent also follows sea level rise. The weathering of silicate or carbonate rocks on land consumes different amounts of atmospheric CO ${}_{\mathrm{2}}$ , with both leading to bicarbonate fluxes into the ocean. These solid Earth processes are not directly coupled with each other. Their implementation in the model might therefore lead to temporal offsets in various variables, to which the sediment module might react in a carbonate compensation feedback. Further details on the model and the time-dependent forcing are found in [74]. Part of this brief model description has been taken from [70].

The following isotopic fractionations are now considered in the BICYCLE-SE model. For this study the whole ${\mathit{\delta }}^{\mathrm{13}}$ C cycle has been revised. While isotopic fractionations are given here in the ${\mathit{\epsilon }}_{(\mathrm{A}-\mathrm{B})}$ -notation (in ‰) they are implemented after [195] in the model as fractionation factors ${\mathit{\alpha }}_{(\mathrm{A}-\mathrm{B})}$ . Both are related following

1 $${\mathit{\epsilon }}_{(\mathrm{A}-\mathrm{B})}=\mathrm{1000}\cdot ({\mathit{\alpha }}_{(\mathrm{A}-\mathrm{B})}-\mathrm{1}).$$

Furthermore, ${\mathit{\alpha }}_{(\mathrm{A}-\mathrm{B})}$ is related to ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ in reservoirs A and B:

2 $${\mathit{\alpha }}_{(\mathrm{A}-\mathrm{B})}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{A}}+\mathrm{1000}}{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{B}}+\mathrm{1000}}}.$$

There is no convention if the initial or final reservoir is given as A or B here; however, here, A is always the final reservoir and B is the initial reservoir of the fractionation process. In some cases, a specific process instead of two reservoirs is mentioned in the subscript, e.g. ${\mathit{\epsilon }}_{(\mathrm{a}\to \mathrm{o})}$ and ${\mathit{\epsilon }}_{(\mathrm{o}\to \mathrm{a})}$ for the atmosphere–ocean gas exchange, for which not only the two different reservoirs but also the direction of the flux play a role in the size of the isotopic fractionation. In that case, the quantified fractionation implies an isotopic depletion connected with the related process for $\mathit{\epsilon }<\mathrm{0}$  ‰.

• Air–sea gas exchange. page997 Using the measurements from we formulate, following in most parts , for the isotopic fractionation during gas exchange to consist of contributions from equilibrium (α eq) and kinetic (α k) fractionation (α total = α eq ⋅ α k). For the atmosphere-to-ocean CO 2 flux, a temperature-dependent equilibrium fractionation of ε a → o eq = ε (aq - g) eq = - 1.31 + 0.0049 ⋅ T C between dissolved (aq) and gaseous (g) CO 2 and a ε a → o k = - 1.08  ‰ is used. Note that ε a → o k differs by - 0.2  ‰ from ε o → a k = - 0.88  ‰ for the ocean-to-atmosphere flux, a necessary correction already given in but to our knowledge only rarely applied. For the reverse ocean-to-atmosphere flux we use the equilibrium fractionation α o → a eq = α (aq - DIC) eq = ∑ i f i α (aq - i) , with f i being the relative shares of CO 2 , HCO 3 - and CO 3 2 - in DIC in the representative ocean box. Furthermore, from the available measurements in , we derive α (aq - HCO 3 -) = α (aq - g) α (HCO 3 - - g) , α (aq - CO 3 2 -) = α (aq - g) α (CO 3 2 - - g) and α (aq - CO 2) = 1 using ε (CO 3 2 - - g) = 7.22 - 0.052 ⋅ T C and ε (HCO 3 - - g) = 10.78 - 0.114 ⋅ T C , with T C being the sea surface temperature in °C.
• Marine biology. The preindustrial marine export production of organic carbon at 100 m water depth is set to 10 PgC yr - 1 (which in the model can increase in glacial periods due to iron fertilisation in the Southern Ocean up to 13 PgC yr - 1 ; Fig. S1d) with a fixed molar rain ratio of organic C:CaCO 3 of 10 : 1. Existing data on fractionation during marine organic matter production (marine photosynthesis) are rather weak in determining if and how it depends on CO 2. Furthermore, as discussed in , some species might contain so-called carbon concentrating mechanisms and use not CO 2 but HCO 3 - as the source of their carbon, in which case a completely different isotopic fractionation during marine photosynthesis (ε (C org - DIC)) would follow. We base our initial formulation of ε (C org - DIC) on the data compilation of δ 13 C POC in , who found a dependency on latitude. Using average preindustrial δ 13 C DIC of + 2.5 ‰ as a starting value and the δ 13 C POC in of - 22 ‰, - 24 ‰ and - 28 ‰ for low, high northern and high southern latitudes, respectively, and approximating ε (C org - DIC) ≈ δ 13 C POC - δ 13 C DIC , we come up with the following isotopic fractionation ε (C org - DIC) of - 24.5  ‰, - 26.5  ‰ and - 30.5  ‰ accordingly. This approximation is motivated by the high uncertainties in δ 13 C POC as documented in. page998 The spread in δ 13 C POC in the data of is huge, ranging from - 15  ‰ to - 35  ‰. Furthermore, they confirmed the finding of earlier studies that δ 13 C POC becomes much more depleted over time than what is explainable by the 13 C Suess effect. Between 1960 and 2010 δ 13 C POC decreased by about 3 ± 4 ‰. The Suess effect shows a decrease in atmospheric δ 13 CO 2 of about 1.5 ‰ during that time , and it is known that in the ocean the Suess effect is decreasing with depth. At the same time, the global mean temperature rose by about 0.8 K. This shift in δ 13 C POC is probably caused by a shift in the composition of the phytoplankton communities. We therefore use the values derived in the previous paragraph from as our preindustrial parameter values of ε (C org - DIC) , to which we add a temperature-dependent part of - 1.5  ‰ for any K the sea surface temperature in the relevant surface-ocean box disagrees with from its preindustrial value. The assumed value fits in the range of recent temperature-dependent δ 13 C POC found in and has been obtained by tuning to simulate δ 13 CO 2 at preindustrial times to be similar to its values at the LGM, as seen in the ice core data (Fig. b). This leads to ε (C org - DIC) at the LGM of - 19.3  ‰, - 20.4 ‰ and - 24.4  ‰ for low, high northern or high southern latitudes, respectively. Data are also rather uncertain for the isotopic fractionation during the formation of CaCO 3. We assume, in agreement with , that 65 % of the CaCO 3 exported in the abyss consists of aragonite and 35 % of calcite. Calcite is either produced by coccolithophores or planktic foraminifera. Some coccolithophore species suggest an enrichment, while others suggest a depletion in δ 13 C in their shells with respect to δ 13 C DIC in the surrounding water. For planktic foraminifera, the CIE is one of various possible processes of isotopic fractionation hypothesised to occur during hard-shell formation. Isotopic fractionation factors are in comparison to ε (C org - DIC) rather small and, in the case of the CIE, species-specific. We therefore choose in the model to set the fractionation during calcite production to be neutral with respect to 13 C, thus ε (cal - DIC) = 0  ‰, but we will consider the CIE in post-processing when comparing simulations with reconstructions. For simplicity and due to missing further evidence for fractionation during aragonite production, ε (ara - DIC) was also kept at 0 ‰. More generally, we keep ε (CaCO 3 - DIC) = 0  ‰. The shallow-water sink of carbonate in corals is assumed to have a δ 13 C that follows after an isotopic fraction of ε (corals - DIC) = - 2  ‰ from the δ 13 C of the DIC in the surface waters. This value is based on a combination of recent data, paleo-data from the Great Barrier Reef and insights from simulations.
• Terrestrial biosphere. On land, isotopic fractionation is only assumed to occur during photosynthesis, with ε (C 3 - CO 2) = - 19  ‰ and ε (C 4 - CO 2) = - 5  ‰ for C 3 (all woody plants and some grasses) and C 4 (some other grasses) photosynthesis, respectively.
• External fluxes to the AOBM subsystem. The volcanic CO 2 outgassing flux is assumed to have a fixed δ 13 C signature (δ 13 C V) of - 5.0  ‰, the typical mean value for volcanic outgassing e.g. , but note that the uncertainty is ± 3  ‰. From the two weathering fluxes based on either silicate or carbonate rocks, only the latter has a contribution which brings new carbon into the system. Here, 50 % of the carbon that is entering the ocean as bicarbonate (the weathering product) has a δ 13 C signature (δ 13 C rock) of + 2 ‰, which is identical to the most likely δ 13 C values in carbonate rocks formed during the Phanerozoic. The carbon for the other half of the carbonate weathering flux and for all of the silicate weathering flux is assumed to come from CO 2 in the soil environment. We therefore assume that this CO 2 is dominated by soil respiration fluxes; therefore, a δ 13 C signature that corresponds to the mean value of the two soil carbon boxes is assumed here. To balance the inflow of 13 C via volcanism and weathering, the model has been tuned for long-term stable mean δ 13 C values in the AOBM subsystem by the following sink: about 6 % of the organic carbon that is exported from the surface boxes into the abyss is assumed to be lost in the sediment. Note that this number has been tuned with the previous version of the 13 C cycle in operation , but it has not been revised thereafter.

The BICYCLE-SE model simulates the global carbon cycle as a function of changing time-dependent physical boundary conditions (forcing), which are nearly identical to the simulations published in [74] and which are also described in that study in detail. Briefly, ocean circulation is prescribed from modern data of the WOCE experiment, while its main temporal changes are restricted to (a) the AMOC, which is reduced from modern/interglacial 16 to 10 Sv during glacial periods (Fig. S1b), and (b) Southern Ocean (SO) vertical deep mixing, which is a function of SO sea surface temperature (Fig. S1c). Ocean and land temperatures are prescribed from reconstructions (Fig. S1e), and ocean salinity is varied as a function of the prescribed sea level (Fig. S1a). Additionally, aeolian iron input in the SO is assumed to follow dust fluxes measured in Antarctic ice cores, which might change marine biology in the SO from an iron-limited to an iron-unlimited regime, increasing glacial export production of organic matter to the deep ocean (Fig. S1d). The standard scenario used here, SEi, is, apart from the revised ${\mathit{\delta }}^{\mathrm{13}}$ C cycle, nearly identical to the scenario SE in [74]. The only difference here is that in the application we revised the applied equatorial sea surface temperature (SST). It has been based in previous applications on changes in planktic ${\mathit{\delta }}^{\mathrm{18}}$ O in only one ODP record. Now we use the SST stack from [10], which is based on a compilation of SSTs from 15 non-polar sediment cores. This leads to only minor changes in atmospheric CO ${}_{\mathrm{2}}$ of less than 5 ppm, but it is important for the ${}^{\mathrm{13}}$ C cycle and its temperature dependencies (isotopic fractionation during atmosphere–ocean gas exchange and during carbon uptake by the marine biology). Simulations are started from interglacial conditions around 210 kyr BP. Scenario SEi0 is only performed to illustrate how the implementation of the temperature dependency in ${\mathit{\epsilon }}_{({\mathrm{C}}_{\mathrm{org}}-\mathrm{DIC})}$ improves the simulated ${}^{\mathrm{13}}$ C cycle, illustrated by plotting atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ against data in Fig. b.

Graph: Figure 5 Simulation results and comparison to data splines for (a) atmospheric CO2 and (b) atmospheric δ13CO2. Results for scenario SEi (standard) and SEi0 are shown. The latter differs from the standard run by a lack of temperature dependency in ε(Corg-DIC).

Simulated changes in the atmospheric record are already not too far away from the reconstructions in scenario SEi, especially in CO ${}_{\mathrm{2}}$ (Fig. a). However, to bring the carbon cycle in the atmosphere and the surface ocean as close to the reconstructions as possible, we perform additional simulations in which the atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ alone (scenario C1) or together with atmospheric CO ${}_{\mathrm{2}}$ (scenario C1CO2) is forced by the reconstructions. Here, we use the data splines as plotted in Fig. a, b [[39], [77]] and ignore the higher resolved data from Taylor Glacier, since these more abrupt changes in ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ are either covered to a large extent during the last 50 kyr in the dynamics of the spline [[11]] or are probably not recorded in our marine sediment records around 70 kyr BP [104]. This implies that internally calculated fluxes are overwritten by changes that are necessary to keep the simulated atmospheric carbon variables identical to the reconstructions. This approach is typically applied in CO ${}_{\mathrm{2}}$ -concentration-driven present-day or future ocean carbon cycle simulations [53]. It has already been used in BICYCLE-SE for ${}^{\mathrm{14}}$ C to obtain radiocarbon in the surface ocean as close to the data as possible during the construction of the most recent marine radiocarbon calibration curve Marine20 [55] and subsequent studies [79]. However, since atmospheric CO ${}_{\mathrm{2}}$ and ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ are normally prognostic variables of the model and their calculated changes should be derived from the model's differential equations followed by a proper integration scheme, this approach slightly violates the mass conservation. It nevertheless guarantees that simulated surface-ocean variables of the carbon cycle within the model realm are as consistent as possible with the atmospheric reconstructions. An overview of the applied simulation scenarios is compiled in Table .

Table 1 Overview of simulation scenarios.

Linear regression was performed with the software MATLAB [173]. The uncertainties of the fits are approximated by root-mean-square errors calculated after $s=\sqrt{\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n(y_i-f_i{)}^{\mathrm{2}}}$ , with $f_i$ being the calculated values according to the linear regression equations. In cases in which the uncertainties in both variables should be considered, we used the function "linfitxy", version 1.2.0.0 [19]. The frequency analysis was performed using R [137], including the function "coh" from the R package seewave, version 2.2.3, calculating coherence.

Reconstructed changes in the late Quaternary carbon cycle are still not completely understood. The ice cores give us a precise picture of atmospheric CO ${}_{\mathrm{2}}$ [[13], [77]] (Fig. a), which in the meantime has also been met reasonably well with various different carbon cycle models [105], [46], [66], [74]. These findings suggest that the main processes responsible for the observed changes on orbital timescales might indeed have been identified, although results are to some extent model-dependent and improvements in details are certainly necessary.

The corresponding atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ , now available over the last 155 kyr [39], is, however, still waiting for a process-based interpretation of all its features (Fig. b). Since ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ helps to pinpoint processes responsible for CO ${}_{\mathrm{2}}$ changes, any simulation that is able to explain one without the other might need to be interpreted with caution. Models suggest that especially physical and biological processes in the Southern Ocean processes robustly influence ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ , while the impact of the Atlantic meridional overturning circulation (AMOC) on ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ seems to be model-dependent [106]. Consequently, the abrupt drop in ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ at the onset of Termination 1 (T1) [[162], [150]] is nowadays understood to be caused by marine processes, while subsequent ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ changes during T1 and its recovery during the Holocene to LGM-like values were potentially related to a mixture of oceanic and terrestrial processes [[75], [11]].

Two largely unexplained features stand out in the 155 kyr ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ record. Firstly, there is a long-term trend of $+$ 0.45 ‰ from the Penultimate Glacial Maximum (PGM) and the Last Glacial Maximum (LGM). When first discovered [154], it was hypothesised that changes in the isotopic composition of solid Earth fluxes or of their intensities or long-term peat build-up might be responsible for them. Secondly, a 0.5 ‰ deep and nearly 20 kyr long minimum centred around 58 kyr BP happened, rather uncorrelated with CO ${}_{\mathrm{2}}$ changes. [39] hypothesise that the ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ minimum might have been partially caused by a change in ocean stratification between Marine Isotope Stage (MIS) 4 and MIS 3, allowing a different amount of isotopically light carbon to be stored in the deep ocean. Recently, high-resolution data of ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ from Taylor Glacier covering 74 to 59.5 kyr BP, including MIS 4 and the drop into the ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ minimum, have been published [104], showing more variability and, between 66 and 60 kyr BP with $-$ 1 ‰, a change twice as large as that previously contained in the smoothed record of [39]. [104] also first performed model simulations in order to understand which processes might be responsible for the reconstructed changes in the carbon cycle. However, to our knowledge, none of the ideas put forward in [154] for the long-term trend in ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ have so far been convincingly and successfully verified with carbon cycle model simulation. Furthermore, 400–500 kyr variability in ${\mathit{\delta }}^{\mathrm{13}}$ C related to slow eccentricity changes found throughout the Cenozoic [127], [146], [99], [184], [126] might be superimposed on faster variations, making a process-based understanding of observed changes in ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ even more challenging.

Sediment cores covering the Anthropocene clearly show that the ${\mathit{\delta }}^{\mathrm{13}}$ C of G. ruber and T. sacculifer shells ( ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}$ and ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}$ ) faithfully reflects changes in ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ caused by the ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ Suess effect [[3], [16]], albeit with a notable offset. This offset might be influenced by the CIE [168], light intensity [167] and the size of the foraminiferal shells [125]. Our new mono-specific stacks from the wider tropical surface ocean of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ (Fig. c) contain a G-IG rise of 0.25 ‰ across T1 but of only 0.15 ‰ across T2, while atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ at the same time rose by 0.1 ‰ (T2) or stayed constant (T1) (Fig. b), showing local minima during terminations in both records. Deep-ocean benthic ${\mathit{\delta }}^{\mathrm{13}}$ C (Fig. d) is approximated here by a stack from six deep Pacific cores [89], which contains a G-IG rise of 0.45 ‰ across both T1 and T2. This value is on the upper end of the 95 % confidence interval of compilations of marine ${\mathit{\delta }}^{\mathrm{13}}$ C changes across T1 [[132], [131]], which suggests a representation of global ocean-wide changes. The marine time series, both from the surface ocean and the deep ocean, also contain wide and deep minima around 60 kyr BP, similarly to the smoothed atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ data of [39] but differently to the more highly resolved Taylor Glacier ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ of [104]. Furthermore, all marine ${\mathit{\delta }}^{\mathrm{13}}$ C data, similarly to the atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ , contain a long-term rise from the PGM to the LGM (about $+$ 0.33 ‰ in the wider tropical surface ocean and $+$ 0.18 ‰ in the deep Pacific; Fig.), which might be potentially connected to the 400–500 kyr variability.

Graph: Figure 6 Scatter plot of our new stacks, (a) Δ(δ13Crub) versus Δ(δ13Csac) and (b) Δ(δ18Orub) versus Δ(δ18Osac). Data stacks without corrections for the CIE are plotted. The time series are restricted to data from the last 150 kyr to allow a comparison later on with simulation results, which were based on the only 155 kyr long atmospheric δ13CO2 record. Linear regressions using only the mean values and also when using uncertainties in both x and y are performed. The root-mean-square error is depicted by s.

Before we start with a deeper model-based interpretation of the ${}^{\mathrm{13}}$ C cycle, we have a closer look at our new isotope stacks. The size of the CIE as detected from laboratory experiments in both species differs by a factor of nearly 2: a change of $-\mathrm{0.0089}$  ‰ and $-\mathrm{0.0047}$  ‰ in ${\mathit{\delta }}^{\mathrm{13}}$ C per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ of [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] for G. ruber and T. sacculifer, respectively, and a change of $-\mathrm{0.0022}$  ‰ and $-\mathrm{0.0014}$  ‰ in ${\mathit{\delta }}^{\mathrm{18}}$ O per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ of [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] for G. ruber and T. sacculifer, respectively [169]. Therefore, if the CIE plays a role in how the isotopes of the surface ocean are recorded in the foraminifera shells on orbital timescales, then the two mono-specific time series in both ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathit{\delta }}^{\mathrm{18}}$ O should differ. At first glance (Fig. a, b) the time series are remarkable similar. A more quantitative evaluation is obtained by calculating the linear regression from scatter plots when results based on one species are plotted against those of the other. Doing so (Fig. ) reveals for ${\mathit{\delta }}^{\mathrm{13}}$ C that, on average, changes are identically recorded in both species. In other words, the linear slope of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ against $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ is 0.98 $(r^{\mathrm{2}}=\mathrm{0.95},s=\mathrm{0.04}\mathrm{‰})$ , or 0.99 $\pm$ 0.03 $(r^{\mathrm{2}}=\mathrm{0.95})$ when considering the uncertainties of our stack during regression. For ${\mathit{\delta }}^{\mathrm{18}}$ O the agreement is only slightly worse: the regression slope of ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}$ against ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}$ is 0.96 $(r^{\mathrm{2}}=\mathrm{0.96},s=\mathrm{0.09}\mathrm{‰})$ , or 0.98 $\pm$ 0.01 $(r^{\mathrm{2}}=\mathrm{0.96})$ with uncertainties. Since $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ are on average recording virtually the same changes, it is difficult to image how the species-specific CIE can play a role here. Due to the small amplitudes of the CIE in ${\mathit{\delta }}^{\mathrm{18}}$ O, it is as yet inconclusive if the CIE plays a role for $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}\right)$ versus $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}\right)$ .

General dynamics of the global carbon cycle in the BICYCLE-SE model have been analysed in detail in [74]. We focus here on the revised ${\mathit{\delta }}^{\mathrm{13}}$ C cycle but see how atmospheric CO ${}_{\mathrm{2}}$ in scenario SEi meets the ice core data in Fig. a. Note that some analysis of ${\mathit{\delta }}^{\mathrm{13}}$ C in the precursor model BICYCLE without solid Earth contributions has been described in [76], who showed that the model misses variations in ${\mathit{\delta }}^{\mathrm{13}}$ C related to periodicities longer than 100 kyr.

Atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ [39] is met by the results from scenario SEi only roughly, including some millennial-scale variations around 50–30 kyr BP and the transition from the LGM to preindustrial, showing some deficit in the second half of T1 and in the Holocene (Fig. b). The PGM-to-LGM trend of 0.45 ‰ and the minimum around 60 kyr BP are both largely unexplained in this simulation. The attribution of changes in ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ to individual processes in the ocean and land carbon cycle has been done before for the precursor model BICYCLE [[75]] and is not repeated here, since the misfit to the data indicates some fundamental shortcomings.

Graph: Figure 7 Simulated surface- and deeper-ocean δ13C time series from scenario SEi (a, b) and scenarios C1 and C1CO2 (c, d) compared with reconstructions. (a, c) Simulated δ13CDIC in the global mean surface ocean and in the wider tropical surface ocean together with simulated atmospheric δ13CO2 (right y axis) are plotted together with our new stacks from the wider tropical surface ocean, with Δ(δ13Crub) and Δ(δ13Csac) shifted by +2.65 ‰ to meet simulated surface δ13CDIC at the LGM. In panels (b) and (d), simulated δ13CDIC for the deep Indo-Pacific (I-Pac), the mean deep ocean and the mean global ocean are plotted together with δ13C from benthic foraminifera stacked from six cores in the deep Pacific (dPac) [89]. In panels (c) and (d), the scenarios C1 (closed lines) and C1CO2 (broken lines) are plotted together. Most of the time the differences between both are so small that the lines are indistinguishable.

The way in which simulated changes in atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ compare to simulated changes in various marine ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ time series is shown for scenario SEi in Fig. a, b. Both global mean surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ and wider tropical surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ show clear similarities with atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ . Here, surface values are area-weighted averages covering either the global ocean or the two equatorial ocean boxes in the case of the wider tropics, which spatially cover a similar area to the sediment cores used for our new stacks, $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ . During glacial times and the onset of deglaciations, the dynamics in the global mean surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ (cyan line in Fig. a) are in close agreement with ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ in the atmosphere (broken black line in Fig. a), while for the later part of the deglaciations and the interglacials the dynamics in the wider tropical surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ (magenta line in Fig. a) fit better to ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ in the atmosphere. This difference is probably explained by the dynamics in the polar oceans. During glacial times, the Southern Ocean is highly stratified with little vertical exchange between the surface ocean and the deep ocean. This stratification breaks down during the terminations and in the interglacials, allowing a faster exchange of tracers between the surface ocean and the deep ocean, leading to smaller surface-to-deep gradients in ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the polar oceans. In other words, the lower deep-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ values have a larger impact on polar surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ during interglacials than during glacials, leading to a divergence in ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the global mean surface and in the wider tropical surface ocean. The scatter plots between atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and either global mean surface or wider tropical surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ show that the latter has the higher correlation (Fig. S2; $r^{\mathrm{2}}=\mathrm{0.82}$ vs.  $r^{\mathrm{2}}=\mathrm{0.59}$ ). Furthermore, frequency analysis showed that the coherence between atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and wider tropical surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ is in periodicities slower than 20 kyr higher than between atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and global mean surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ (Fig. S3a). This implies that simulations which agree in atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ with reconstructions (which will be achieved later on in scenarios C1 and C1CO2) should contain a very likely realisation of ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the wider tropical surface ocean. A comparison of these simulated time series with our new mono-specific ${\mathit{\delta }}^{\mathrm{13}}$ C stacks should therefore enable us to address if and how ${\mathit{\delta }}^{\mathrm{13}}$ C has been modified during hard-shell formation. For scenario SEi the misfit in simulated wider tropical surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ and the new ${\mathit{\delta }}^{\mathrm{13}}$ C reconstructions (Fig. a) is large, but it is as yet unclear if this discrepancy can be explained by the CIE or by other processes.

To understand how representative the reconstructed ${\mathit{\delta }}^{\mathrm{13}}$ C stack from benthic foraminifera in six deep Pacific cores [89] might be, we compare it with various different simulated time series: ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the deep Indo-Pacific, in the mean deep ocean and in the mean ocean (Fig. b). Here, deep-ocean results from the model refer to ocean boxes that contain waters deeper than 1 km. As expected, the deep Indo-Pacific contains the end member of the ${\mathit{\delta }}^{\mathrm{13}}$ C cycle with the most depleted values. The mean deep-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ is offset by 0.2 ‰–0.4 ‰ towards more positive values and shows larger G-IG amplitudes than ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ does in the deep Indo-Pacific. The mean ocean is again 0.2 ‰–0.4 ‰ more positive in ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ than the mean deep ocean, again with smaller G-IG amplitudes of 0.53 ‰ across T1. This number compares ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the last 6 kyr with the mean at the LGM (23–19 kyr BP), similarly as in [132], who proposed a mean ocean rise in ${\mathit{\delta }}^{\mathrm{13}}$ C by $\mathrm{0.34}\pm \mathrm{0.19}$  ‰. However, be aware that in [132] the CIE in benthic foraminifera as deduced in [152] is not included. This suggests that the reconstructions are potentially recording a smaller G-IG change in ${\mathit{\delta }}^{\mathrm{13}}$ C than how ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the deep ocean might have changed.

When discussing the results of scenario SEi (Fig. a), we have shown that, once changes in the atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ are met by the simulations, the model should then also give a reasonable answer for what ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the wider tropical surface ocean might have looked like. Furthermore, the close agreement in simulated and reconstructed atmospheric CO ${}_{\mathrm{2}}$ (Fig. a) suggests that the assumed carbon cycle changes in our approach might be one possible realisation that is not too far away from the real-world changes. However, the misfit between simulation results from scenario SEi and reconstruction in the ${\mathit{\delta }}^{\mathrm{13}}$ C cycle, where linear regressions between simulations and reconstructions found no correlation at all ( $r^{\mathrm{2}}\le \mathrm{0.02}$ ; Fig. S4a, b), is not easily fixed. To improve our results, we force the model with the atmospheric records (scenario C1 only using ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and scenario C1CO2 using both ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and CO ${}_{\mathrm{2}}$ ) to have conditions in the surface ocean as close to reconstructions as possible. Doing so leads to even tighter correlations between simulated atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and simulated ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the surface ocean than what we obtained for scenario SEi: the $r^{\mathrm{2}}$ correlations between these variables are in scenarios C1 and C1CO2, with prescribed atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ $\ge$  0.77 and $\ge$  0.88 for global mean surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ and wider tropical surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ , respectively (Fig. S2). Again, the coherence is higher between atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and the wider tropical surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ than between atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and the global mean surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ (Fig. S3b). Furthermore, in both scenarios, the changes in simulated ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the wider tropical surface ocean agree remarkably well ( $r^{\mathrm{2}}$ between 0.76 and 0.78; Fig. S4c–f) with changes in our new stacks $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ without consideration of the CIE (Fig. c), at least on orbital timescales. This effect is also seen by the rise in coherence between simulated wider tropical surface ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ and both our stacks from less than 0.1 (scenario SEi) to higher than 0.7 (scenario C1CO2) in the 41 and 100 kyr bands (Fig. S3c, d), while in the precession bands (19 and 23 kyr) the coherence stayed below 0.6. Some more abrupt changes contained in the simulations are not recorded in the reconstructions, probably because bioturbation in the surface sediments, together with the stacking procedure, prevent our marine records from successfully resolving millennial-scale features. Thus, our forcing of atmospheric carbon records with data therefore seems to be a promising approach to obtain simulated surface ocean in agreement with reconstructions for the slow-frequency bands (41 kyr and beyond), while it seems to fail for precession and faster changes. When forcing atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ by data, the temperature-dependent isotopic fractionation during marine photosynthesis in ${\mathit{\epsilon }}_{({\mathrm{C}}_{\mathrm{org}}-\mathrm{DIC})}$ is only of minor importance for the simulated surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ . If this effect is switched off, the ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the wider tropical surface ocean differs in general by less than 0.05 ‰ from the values in scenario C1.

Furthermore, deep-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ is now, on an orbital timescale, also in better agreement with the data (Fig. d), the $r^{\mathrm{2}}$ of a linear regression between simulated deep Indo-Pacific ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ , and the reconstructed deep Pacific rises from 0.49 for scenario SEi to 0.77 and above for the scenarios forced by atmospheric carbon records (Fig. S5), although the rise in mean ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ during T1 has now been increased to 0.59 ‰. Considering a CIE of $-\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{3}}$  ‰ per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ of [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] disturbance for epibenthic foraminifera [152], simulated variations in deep-ocean [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] of $+\mathrm{20}$   $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ [74] would translate to a comparably small reduction in deep Pacific benthic ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ of up to 0.05‰. While the timing of changes in deep-ocean [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] with highest values during the deglaciation is crucial to assess how such a benthic CIE would reduce the existing data–model mismatch, a more thorough assessment of the benthic CIE would require the comprehensive compilation of benthic ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ time series in different ocean basins, which is beyond the scope of this study. Note that the approximated amplitude of this benthic CIE is close to the measurement error of benthic ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ .

Although the initial analysis of our results when forced with atmospheric records already suggests only a minor, if any, role for the CIE in the interpretation of stacked mono-specific ${\mathit{\delta }}^{\mathrm{13}}$ C on orbital timescales, in the following section we make a more quantitative assessment. The CIE has not yet been implemented in the ${}^{\mathrm{13}}$ C cycle of the model, but it is only investigated here in post-processing. The carbonate ion concentration of either global mean surface waters or wider tropical mean surface waters in our simulations is tightly anti-correlated with atmospheric CO ${}_{\mathrm{2}}$ ( $r^{\mathrm{2}}\ge \mathrm{0.93}$ ; Fig. S6), which is a consequence of the marine carbonate system [195]. Both scenarios C1 and C1CO2 lead to rather similar results here, which suggests that the CO ${}_{\mathrm{2}}$ forcing in scenario C1CO2 and its violation of mass conservation are perturbing the carbon cycle only slightly. To stay as closely as possible to the reconstructions, we nevertheless continue in the following section by using results from scenario C1CO2, but results differ only slightly when based in scenario C1; thus our conclusions are independent from this choice.

Thus, CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ in the wider tropical surface ocean in the simulation typically falls from maximum glacial values of $\sim$  320  $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ to interglacial minima of $\sim$  250  $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ across both Terminations 1 and 2 (Fig. a). This translates into a potential CIE of about 0.62 ‰ (Fig. b) for G. ruber when we use the slope of $m=-\mathrm{0.0089}$  ‰ per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ change in $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ and of 0.33 ‰ for T. sacculifer (slope of $m=-\mathrm{0.0047}$  ‰ per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ change in $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ ) [169]. The $y$ axis intercepts of the complete regressions for the CIE is determined in order to have a maximum agreement between reconstructions and simulations during the LGM. When comparing the potential CIE to the simulated LGM-to-preindustrial (PRE) amplitude of only 0.16 ‰ in wider tropical surface waters (Fig. c), the CIE-to-G-IG ratios are between a factor of 2 and 4, and CIE signals should clearly stand out in the paleo-records. If we add this CIE to our simulated mean equatorial surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ (Fig. c), we end up with time series, which should compare well with the mono-species stacks of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ (Fig. d). However, this is not the case. The $r^{\mathrm{2}}$ in the linear regressions between CIE-corrected ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in wider tropical surface waters and reconstructions is reduced to 0.54 (G. ruber) and 0.68 (T. sacculifer), while it had been $\ge$  0.76 without CIE correction (Figs. S4, S7). When plotting results as hypothetically recorded in both species against each other, we obtain a slope of 1.26 (Fig. S8a). The slope between the stacked mono-specific ${\mathit{\delta }}^{\mathrm{13}}$ C time series without further correction for a CIE was $\sim$  0.99 (Fig. a). The consideration of the CIE did not lead to time series which agree better with each other. Thus, we conclude that both species, G. ruber and T. sacculifer, are already good recorders of changes in ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in wider tropical surface-ocean waters on orbital timescales.

Graph: Figure 8 Calculating the suggested carbonate ion effects (CIE) on G. ruber and T. sacculifer. The left-hand column shows the effects on δ13C, while the right column shows the effects on δ18O. (a, e) Surface-ocean [CO32-]. (b, f) Potential CIE using slopes from [169]. (c, g) Surface-ocean conditions when ignoring the CIE or (d, h) when considering the CIE. Mean anomalies (±1 SE) of the isotope stacks are calculated with respect to the mean of 21–19 kyr BP (vertical blue band). Simulations use the results from scenario C1CO2. Different surface-ocean areas are distinguished: North Atlantic (NAtl; north of 50° N), equatorial Atlantic (EqAtl; 40° S–50° N), Southern Ocean (SO; south of 40° S), equatorial Indo-Pacific (EqIPac; 40° S–40° N) and North Pacific (NPac; north of 40° N). The mean wider tropical ocean in the model is the mean from the equatorial boxes (eq-mean).

The focus of this study is on stable carbon isotope ${\mathit{\delta }}^{\mathrm{13}}$ C. However, during the construction of our mono-specific wider tropical stacks of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sac}}\right)$ , the corresponding stacks of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}\right)$ are easily generated by-products initially used to cross-check the applied age models. However, these ${\mathit{\delta }}^{\mathrm{18}}$ O data give us the possibility to also have a closer look at the role of the CIE in the recording of oxygen isotopes in foraminiferal shells. For that endeavour, we need a background time series of ${\mathit{\delta }}^{\mathrm{18}}$ O which represents the signals when not modified by the CIE. Such a mean ${\mathit{\delta }}^{\mathrm{18}}$ O in the wider tropical surface ocean should record the same sea-level-related variations as the average global ocean, but it might differ in the recorded temperature effect if the change in the average wider tropical sea surface temperature differed from the mean ocean temperature (MOT) change. [134] showed that the LGM-to-PRE change in the MOT derived from the model-based interpretation of noble gas reconstructions in ice cores is $\mathrm{2.1}\pm \mathrm{0.7}$  K. The reconstructed rise in the MOT is slightly higher when ignoring the effect of past saturation changes on noble gases [156]. The data assimilation effort in LGM temperature changes by [176] is broadly in agreement with the MOT change of [134] and proposes that the tropical (30° S to 30° N) sea surface was around $\mathrm{2.6}$  K colder at the LGM than at the PRE, agreeing within the uncertainties with the MOT change. To a first order, we therefore assume that the planktic foraminifera should record the same temperature effect in ${\mathit{\delta }}^{\mathrm{18}}$ O as contained in the mean ocean. Thus, the global ocean ${\mathit{\delta }}^{\mathrm{18}}$ O calculated from stacking benthic time series [91] represents the CIE-free background against which we compare our new $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}\right)$ stacks.

From the simulated LGM-to-PRE change in mean wider tropical surface-ocean CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ of about $-\mathrm{70}$   $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ (Fig. e) and the laboratory-based amplitudes of the CIE ( $-$ 0.0022 ‰ and $-$ 0.0014 ‰ change in ${\mathit{\delta }}^{\mathrm{18}}$ O per $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ for G. ruber and T. sacculifer, respectively [169]), we determined that $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{rub}}\right)$ and $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sac}}\right)$ should record the changes since the LGM by $+$ 0.15 ‰ and $+$ 0.10 ‰ differently than how ${\mathit{\delta }}^{\mathrm{18}}$ O in the surface waters truly changed (Fig. f). Compared to the G-IG amplitude in mean ocean ${\mathit{\delta }}^{\mathrm{18}}$ O of $-$ 1.65 ‰ (Fig. g), these potential CIEs represent corrections of $-$ 9 % and $-$ 6 %, a difference of 3 % which might be difficult to detect in the paleo-records. A linear regression through a scatter plot of ${\mathit{\delta }}^{\mathrm{18}}$ O + CIE ${}_{\mathrm{rub}}$ versus ${\mathit{\delta }}^{\mathrm{18}}$ O  $+$  CIE ${}_{\mathrm{sac}}$ has a slope of 0.97 ( $r^{\mathrm{2}}=\mathrm{1.00}$ ; Fig. S8b), which is indistinguishable from the slope obtained from regression through the data stacks (Fig. b), while the slope when considering the CIE should move to unity (indicating that both species were recording the same signal underneath the CIE) if the effect plays an important role during data interpretation. The evidence for or against the CIE in ${\mathit{\delta }}^{\mathrm{18}}$ O from both data and models is therefore inconclusive.

The CIE for ${\mathit{\delta }}^{\mathrm{13}}$ C and ${\mathit{\delta }}^{\mathrm{18}}$ O recorded in planktic foraminifera was first identified in laboratory experiments [[168]], and it was, based on theory, suggested for both isotopes that the underlying processes are directly related to the pH in the surrounding seawater during hard-shell formation [[196], [194]]. However, these theoretical studies were already unable to confirm the full range of the CIE as contained in the experiments. Furthermore, according to [14], it is impossible to determine if pH or [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] is responsible for the observed fractionation effects. If this theoretical understanding is correct, we would expect to see the CIE in neither or both isotopes in the our mono-specific stacks. Thus, although the interpretation of ${\mathit{\delta }}^{\mathrm{18}}$ O with respect to the CIE is uncertain due to the signal-to-noise ratio, we argue, based on the clear evidence of a lack of the CIE in the recording of ${\mathit{\delta }}^{\mathrm{13}}$ C in G. ruber and T. sacculifer, that there is probably also no significant CIE contained in the ${\mathit{\delta }}^{\mathrm{18}}$ O time series of both species. This finding argues against the suggestion of [169] that the CIE and ${\mathit{\delta }}^{\mathrm{13}}$ C time series from G. ruber and T. sacculifer might be used to calculate a record of surface-ocean [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ]. Furthermore, we suggest using our new stack of $\mathrm{\Delta }\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rub}}\right)$ as representative of ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in the wider tropical surface ocean.

Various possible explanations for the lack of a CIE on orbital timescales exist. Firstly, it might be that the isotopic fractionation during hard-shell formation in G. ruber and T. sacculifer is rather insensitive to [CO ${}_{\mathrm{3}}^{\mathrm{2}-}$ ] in the range of interest (250–320  $\mathrm{\mu }$ mol kg ${}^{-\mathrm{1}}$ ). Such an insensitivity has been suggested for other species [14], but, due to a lack of published data (the slopes of the CIE in G. ruber and T. sacculifer were only summarised in [169], while underlying experiments have never been published in the peer-reviewed literature), it cannot be properly checked for the two species investigated here. Secondly, not the CIE, but alternatively the incorporation of respired CO ${}_{\mathrm{2}}$ (depleted in ${\mathit{\delta }}^{\mathrm{13}}$ C) during shell formation, might be responsible for the observed isotope data in laboratory experiments performed with Orbulina universa and Globigerina bulloides[14]. This process might also play a role in G. ruber and T. sacculifer, but it would only explain observed effects in ${\mathit{\delta }}^{\mathrm{13}}$ C and not in ${\mathit{\delta }}^{\mathrm{18}}$ O. However, since our stacks are inconclusive with respect to the CIE and ${\mathit{\delta }}^{\mathrm{18}}$ O, they might be of relevance here. A third explanation might be related to homeostasis. In symbiont-bearing planktic foraminifera, such as G. ruber and T. sacculifer, the pH at the shell surface critically depends on photosynthesis and hence light levels and symbiont density [62]. In order to facilitate calcification, G. ruber and T. sacculifer may actively influence the pH at the shell surface by seeking specific (optimum) light levels through vertical migration, thereby keeping the CIE constant over time. Planktic foraminifera are known to move vertically in the water column [67]. Vertical migration to optimise both nutrient uptake and light has been proposed to play an important role in phytoplankton by modelling [189], an effect which has recently been supported by field data [199]. We speculate that similar behaviour could occur in the two planktic foraminifera species. Indeed, [61] and [27] found that ${\mathit{\delta }}^{\mathrm{18}}$ O in various planktic foraminifera (including G. ruber and T. sacculifer) is best explained by also considering calcification in waters deeper than their expected living depth.

It is too early to be able to generalise our finding that on orbital timescales the CIE plays no role in the interpretation of signals in planktic foraminifera in paleo-records. For that endeavour, more mono-specific stacks are necessary, preferably from conceptually different foraminifera species without symbionts or spines, as these might potentially show a different behaviour with respect to light (and pH) optimisation. However, our findings might suggest that previous studies on planktic ${\mathit{\delta }}^{\mathrm{13}}$ C, which ignored the CIE [98], [96], might not be biased.

Our carbon cycle simulations confirm that atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and mean surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ are tightly related to each other, highlighting the importance of air–sea gas exchange for carbon isotopes. This is not entirely new, and it has already been discussed before [98], [157], [133]. However, the ${}^{\mathrm{13}}$ C cycle is more complex than stated previously [[98], [58], [133]]. These studies suggested that one might calculate a mean surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ as a function of atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ and a temperature-dependent fractionation during gas exchange. We assumed here, based on modern data from [180], that species composition, and therefore isotopic fractionation during marine photosynthesis, might also be temperature-dependent, having an important impact on surface-ocean ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ . Furthermore, our simulation results show that ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}$ in polar oceans and in the wider tropical surface ocean have a different and time-dependent relation to atmospheric ${\mathit{\delta }}^{\mathrm{13}}$ CO ${}_{\mathrm{2}}$ .

Finally, since our simulations were forced by atmospheric carbon records, we are unable to identify specific processes being responsible for the simulated changes in the ${}^{\mathrm{13}}$ C cycle. Recent climate simulations [192] emphasise the importance of the 405 kyr eccentricity cycle in tropical hydroclimate. It therefore seems reasonable that the missing long-term variability in ${\mathit{\delta }}^{\mathrm{13}}$ C in our setup might indeed be connected to weathering fluxes, as proposed before [154], [184], something which needs to be tested in more detail in future carbon cycle simulation studies.

The following citations contain the original publications for the data used to produced our stacks. More details are found in Table S1: [1]; [2]; [4]; [5]; [6]; [7]; [8]; [17]; [22]; [23]; [25]; [26]; [28]; [31]; [32]; [33]; [34]; [35]; [36]; [37]; [38]; [45]; [47]; [48]; [49]; [50]; [51]; [52]; [56]; [57]; [59]; [60]; [64]; [65]; [68]; [69]; [80]; [81]; [84]; [85]; [86]; [87]; [94]; [97]; [101]; [107]; [108]; [109]; [110]; [111]; [113]; [114]; [115]; [116]; [117]; [119]; [122]; [128]; [129]; [130]; [135]; [136]; [139]; [141]; [142]; [145]; [147]; [148]; [149]; [153]; [155]; [158]; [159]; [160]; [161]; [170]; [171]; [172]; [174]; [175]; [177]; [178]; [179]; [182]; [183]; [185]; [186]; [187]; [188]; [193]; [200].

The supplement related to this article is available online at: https://doi.org/10.5194/cp-20-991-2024-supplement.

This work is supported by the German Federal Ministry of Education and Research (BMBF) as a Research for Sustainability initiative (FONA; , last access: 18 April 2024) through the PalMod project. We thank Peter U. Clark for providing the SST data connected to [10] and Jelle Bijma for discussions.