Isotope re-equilibration experiments

A synthetic pure 18O seawater solution was prepared by adding NaCl and NaHCO 3 to pure H 2 18O until reaching 0.55 and 0.003 mol L−1, respectively. Sub-modern planktonic foraminifera (G. bulloides) from modern sediments of the Gulf of Lion, France, with bulk δ 18O ~ 1.35 ± 0.05‰ (VPDB) were rinsed three times in ethanol using an ultrasonic bath. Gold capsules, each containing 160 µg of cleaned foraminifera tests (∼12 specimens) and 100 µL of this synthetic seawater solution, were sealed using a Lampert PUK 4 welding machine and placed in Parr autoclaves for 82 days at T = 300 °C and P = 200 bars. Note that with a 18O/(16O+18O) ratio of the foraminifera test close to 0 (because of the rarity of 18O) when the ratio of the solution is close to 1, the mass balance is such that the water δ 18O can be considered constant during the experiments, regardless of the magnitude of the isotope exchange during the experiments.

Speciation calculations were performed by running the geochemical code CHESS34 using the thermodynamic database of the EQ3/6 code35. Activity coefficients for aqueous species were calculated using the Davies equation36; the electrical balance was achieved using the H+ concentration. At ambient temperatures (293 K), the solution is slightly undersaturated with respect to calcite (within uncertainty bounds). At this temperature, assuming that the dissolution of calcite is thermodynamically controlled (i.e. neglecting kinetic barriers), equilibrium would be reached after the dissolution of <2% of the tests. At the experimental temperature (573 K), because of the retrograde solubility of calcite, saturation with respect to calcite is reached for a dissolution progress corresponding to only 0.1% of the tests. Therefore, we can assume that the NanoSIMS images reflect an isotope re-equilibration very close to chemical equilibrium, with little or no contributions from secondary crystallisation from the bulk solution.

Because the carbonate ‘building blocks’ of planktonic and benthic foraminifera are the same at microscales29,30, we confidently use the results of experiments conducted on planktonic foraminifera as representing the results for both planktonic and benthic foraminifera.

Characterisation techniques

SEM observations were performed on foraminifera tests deposited on aluminium stubs and with 5-nm thick coatings of gold using a SEM-FEG Ultra 55 Zeiss (IMPMC, Paris, France) microscope operating at a 2-kV accelerating voltage and a working distance of 2 mm for imaging with secondary electrons and at a 15-kV accelerating voltage and a working distance of 7.5 mm for imaging with backscattered electrons and EDXS mapping.

Isotope maps were produced on triplicate samples with the Cameca NanoSIMS 50 (IMPMC, Paris, France). A thorough technical explanation of the NanoSIMS instrument has been provided in a recent review23. Briefly, the NanoSIMS ion microprobe is a secondary ion mass spectrometry instrument characterised by an extremely high spatial resolution, a high sensitivity, a high mass-resolving power and a multicollection capability34. Using a focused primary beam of 133Cs+ ions, secondary ions were sputtered from the sample surface, typically to a depth of ∼100 nm. 18O– and 16O– ions from the sample were simultaneously detected (multicollection mode) by electron multipliers at a mass-resolving power of ∼9000 (M/∆M). At this mass-resolving power, the measured secondary ions were resolved from the potential interferences. Images were obtained from a presputtered surface area (the same surface area previously imaged using SEM) by rastering the primary beam across the sample surface. The primary beam was focused to a spot size of ~150 nm, and the pixel size was adjusted so that it was smaller than the size of the primary beam. An electron gun supplied the electrons to the sputtered surface during the analysis to compensate for positive charge deposition from the primary beam and to minimise the surface charging effects. Imaging data were processed using custom-made software (LIMAGE, L. Nittler, Carnegie Institution of Washington).

Modelling the experiments

In the present experiments, the driving force for the isotope re-equilibration through solid-state diffusion is the large 18O excess of the fluid. The most generic expression describing the temperature dependence of the diffusion coefficient can be written as follows:

$$D = {D_0}{\rm{exp}}\left( { - \frac{{{\rm{Ea}}}}{{RT}}} \right),$$ (1)

where D 0 is the intrinsic diffusion constant, Ea is the activation energy of oxygen diffusion in calcite, R is the ideal gas constant and T is the absolute temperature. Although the bulk diffusion of oxygen in calcite is the sum of the contributions of both the grain boundary and volume diffusion, it is commonly assumed that the volume diffusion is slow and can be ignored37,38.

Applied to the present case, Eq. (1) becomes:

$${D_{{\rm{foram}}}} = {D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right),$$ (2)

where T xp is the absolute temperature of the experiment and Ea foram is the activation energy of oxygen solid-state diffusion in foraminifera calcite. No consensus exists in the literature on the value of this last parameter (see below).

In the present experiments, the initial and boundary conditions can be expressed as follows:

$$C\left( {x,t = 0} \right) = 0,\,\forall \,x >0$$ (3)

$$C\left( {x = 0,t} \right) = {C_0} = 1$$ (4)

with C standing for the 18O concentration (i.e. the molar 18O/(16O+18O) ratio) in foraminifera calcite, t is time and x is the length (positive distance from the fluid/solid interface into the solid).

The diffusion equation for C reads as follows (Fick’s second law):

$$\frac{{\partial C}}{{\partial t}} = {D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right)\frac{{{\partial ^2}C}}{{\partial {x^2}}},$$ (5)

According to Crank39, Because the calcite fraction quantitatively affected by the diffusion process is generally small compared to the total calcite size/volume, the analytical solution for the diffusion into an infinite one-dimensional medium can be applied to calculate the 18O concentration along the profile after the experimental duration t xp :

$$C\left( x \right) = {C_0}\left( {1 - {\rm{erf}}\left( {x{\rm{/}}\left( {2\sqrt {{D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right){t_{{\rm{xp}}}}} } \right)} \right)} \right).$$ (6)

Because the error function is defined as follows:

$${\rm{erf}}\left( y \right) = \frac{2}{{\sqrt \pi }}\mathop {\int }

olimits_0^y \exp \left( { - {\tau ^2}} \right){\rm{d}}\tau ,$$ (7a)

with

$$y = \frac{x}{{2\sqrt {{D_{0,{\rm{foram}}}}\exp \left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right)\,{t_{{\rm{xp}}}}} }},$$ (7b)

its derivative can be written as follows:

$$\frac{\partial }{{\partial {\rm{y}}}}\left( {{\rm{erf}}\left( y \right)} \right) = \frac{2}{{\sqrt \pi }}\exp \left( { - {y^2}} \right).$$ (8)

The flux of matter F entering the calcite domain at the interface x = 0 can be expressed as follows:

$$F\left( {0,t} \right) = - {D_{{\rm{foram}}}}\frac{{\partial C}}{{\partial x}},$$ (9)

with (according to Eqs. (6) and (8)):

$$\frac{{\partial C}}{{\partial x}} = \frac{{\partial C}}{{\partial y}}\frac{{\partial y}}{{\partial x}} = \frac{{ - {C_0}{\rm{exp}}( - {{\rm{y}}^2})}}{{\sqrt {\pi {D_{{\rm{foram}}}}t} }},$$ (10)

leading to the following expression for the flux F at the interface x = 0:

$$F\left( {0,t} \right) = \frac{{{C_0}\sqrt {{D_{{\rm{foram}}}}} }}{{\sqrt {\pi t} }},$$ (11)

The cumulative flux (F c ), which represents the amount of 18O that has entered the calcite domain during the experimental duration t xp , can be expressed as follows:

$${F_c}\left( {0,{t_{{\rm{xp}}}}} \right) = \mathop {\int }

olimits_0^{{t_{{\rm{xp}}}}} F\left( {0,\tau } \right){\rm{d}}\tau = {C_0}\frac{{\sqrt {{D_{{\rm{foram}}}}} }}{{\sqrt \pi }}\mathop {\int }

olimits_0^{{t_{{\rm{xp}}}}} \frac{1}{{\sqrt \tau }}{\rm{d}}\tau = 2{C_0}\frac{{\sqrt {{D_{{\rm{foram}}}}{t_{{\rm{xp}}}}} }}{{\sqrt \pi }}.$$ (12)

From a physical standpoint, this result illustrates that the integral of any diffusion profile that satisfies the abovementioned boundary conditions is equivalent to that of a step function of the height C 0 with a length d of:

$$d = \frac{2}{{\sqrt \pi }}\sqrt {{D_{{\rm{foram}}}}{t_{{\rm{xp}}}}} .$$ (13)

Because the initial 18O/(16O+18O) value of the foraminifera is close to 0, the 18O/(16O+18O) value of the foraminifera calcite crystals at the end of the experiment (Q) is equivalent to the ratio of the re-equilibrated volume to the initial volume. Assuming that the foraminifera calcite crystals have a simple spherical geometry with a radius r 0 and that d << r 0 (as suggested by NanoSIMS imaging), Q can be written as follows:

$$Q = \frac{{\frac{4}{3}\pi r_0^3 - \frac{4}{3}\pi {{\left( {{r_0} - d} \right)}^3}}}{{\frac{4}{3}\pi r_0^3}}.$$ (14)

Combining Eqs. (2), (13) and (14), the mean radius of the calcite grains (r 0 ) can thus be expressed as follows:

$${r_0} = \frac{d}{{1 - {{\left( {1 - Q} \right)}^{1{\rm{/}}3}}}} = \frac{2}{{\sqrt \pi }}\frac{{\sqrt {{D_{{\rm{foram}}}}{t_{{\rm{xp}}}}} }}{{1 - {{\left( {1 - Q} \right)}^{1/3}}}} = \frac{2}{{\sqrt \pi }}\frac{{\sqrt {{D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right){t_{{\rm{xp}}}}} }}{{1 - {{\left( {1 - Q} \right)}^{1/3}}}}.$$ (15)

Given the actual calcite grain size of foraminifera (50–250 nm)29,30, i.e. r 0 values of 25–125 nm, the quantity of 18O that diffused within the foraminifera tests during the experiments can be used to estimate the activation energy and the intrinsic diffusion coefficient of oxygen solid-state diffusion in foraminifera calcite (Ea foram and D 0,foram ). A large range of several tens of kJ mol−1 exists in the literature for the activation energy of oxygen solid-state diffusion in calcite aggregates. Farver and Yund27 estimated Ea values of 127 ± 17 kJ mol−1 (D 0,Farver = 7.6*10−9 m2 s−1), while Anderson26 suggested that diffusion within calcites can even proceed with an activation energy as low as 70 kJ mol−1 (D 0,Anderson = 4.6*10−20 m2 s−1) at low temperatures. Applied to the present experiments, Eq. (15) yields r 0 values of 7 nm (i.e. a calcite grain size of ∼14 nm, assuming spherical grains) with the values reported by Anderson26, while an r 0 value of 8.25 µm (i.e. a calcite grain size of ∼16.5 µm) is obtained when using the values reported by Farver and Yund27. This discrepancy from the actual grain size of foraminifera is easily explained by (1) the large uncertainty (typically several tens of kJ mol−1) in the experimental determination of the Ea values and (2) the large dispersion of foraminifera grain sizes and effective grain boundaries, which require intermediate values of Ea foram and D 0,foram computed as weighted average defined by:

$$\left\{ {\begin{array}{*{20}{c}} {{\rm{E}}{{\rm{a}}_{{\rm{foram}}}} = x{\rm{E}}{{\rm{a}}_{{\rm{Farver}}}} + \left( {1 - x} \right){\rm{E}}{{\rm{a}}_{{\rm{Anderson}}}}} \\ {{\rm{log}}\left( {{D_{0,{\rm{foram}}}}} \right) = x{\rm{log}}\left( {{D_{0,{\rm{Farver}}}}} \right) + \left( {1 - x} \right){\rm{log}}\left( {{D_{0,{\rm{Anderson}}}}} \right)} \end{array}} \right.,$$ (16)

with x varying between 0 and 1.

Figure 5 shows that x has to be between 0.18 and 0.41 to yield foraminifera calcite grain size values between 50 and 250 nm (i.e. r 0 values of 25–125 nm). This range corresponds to Ea foram values ranging from 81 to 94 kJ mol−1 and D 0,foram values ranging from 4.8 × 10−18 to 1.83 × 10−15 m2 s−1.

Fig. 5 Determination of the Ea foram and D 0,foram values for the oxygen solid-state diffusion in foraminifera calcites. a Foraminifera calcite grain sizes (2*r 0 values) obtained using Eq. (15) for the different x values (Eq. (16)). b D 0,foram and Ea foram values obtained using Eq. (16) for the different x values Full size image

Modelling natural settings

We developed a numerical model to estimate the present-day δ 18O of fossil tests of benthic foraminifera assumed to have formed during the Cretaceous and the Cenozoic in equilibrium with a seawater/sediment interface at 3.5 °C and that underwent partial isotope re-equilibration via solid-state diffusion during sediment burial. This model was also used to estimate the present-day δ 18O values of fossil tests of planktonic foraminifera assumed to have formed during the Cretaceous and the Cenozoic at different latitudes, i.e. at equilibrium with seawater at different temperatures, and that underwent a similar burial history.

Because the volumes impacted by re-equilibration may not be negligible compared to the size of the foraminifera calcite building blocks, it is preferable to solve the diffusion equation with a spherical geometry rather than with a simpler 1D semi-infinite medium. Because the spherical diffusion equation expressed below does not have a simple analytical solution, the diffusion equation was numerically solved using the finite difference with an implicit time discretisation.

The driving force for isotope re-equilibration through solid-state diffusion over geological times within sediments is the temperature-dependent O isotope fractionation between foraminifera test calcite and seawater at thermodynamic equilibrium (\({\Delta _{{\rm{cc}} - {\rm{sw}}}}(\theta )\)), where θ stands for the temperature in °C. A number of formulations of \({{\rm{\Delta }}_{{\rm{cc}} - {\rm{sw}}}}(\theta )\) can be found in the literature2. The most widely used equation is the one reported by Anderson and Arthur 3), which relates \({\Delta _{{\rm{cc}} - {\rm{sw}}}}\) and θ as follows:

$$\theta = 16.0 - 4.14\left( {{\Delta _{{\rm{cc}} - {\rm{sw}}}}} \right) + 0.13{\left( {{\Delta _{{\rm{cc}} - {\rm{sw}}}}} \right)^2},$$ (17)

with

$${\Delta _{{\rm{cc}} - {\rm{sw}}}}\left( \theta \right) = {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}} - {\delta ^{18}}{{\rm{O}}_{{\rm{sw}}}},$$ (18)

where \({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\) is the O isotope composition of the test calcite (in ‰, relative to VPDB) and \({\delta ^{18}}{{\rm{O}}_{{\rm{sw}}}}\) is the O isotope composition of seawater (in ‰, relative to VSMOW). This equation allows the \({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\) to be written as a function of \({\delta ^{18}}{{\rm{O}}_{{\rm{sw}}}}\) and T (in K):

$${\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}} = \frac{{4.14 - \sqrt {{{\left( {4.14} \right)}^2} - 4{\rm{*}}0.13{\rm{*}}\left( {16 - T + 273.15} \right)} }}{{2{\rm{*}}0.13}} + {\delta ^{18}}{{\rm{O}}_{{\rm{sw}}}}.$$ (19)

Assuming spherical calcite crystals, the mathematical model describing diffusion in spherical coordinates is:

$$\frac{{\partial C}}{{\partial t}} = {D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{RT}}} \right)\frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\left( {{r^2}\frac{\partial }{{\partial r}}} \right),$$ (20)

with the following initial and boundary conditions in natural settings:

$${\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,t = 0} \right) = {{\rm{\Delta }}_{{\rm{cc}} - {\rm{sw}}}}\left( {T(t = 0)} \right) + {\delta ^{18}}{{\rm{O}}_{{\rm{sw}}}}\left( {t = 0} \right),\,\forall \,r \in [0;{r_0}],$$ (21)

$${\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r = {r_0},t} \right) = {{\rm{\Delta }}_{{\rm{cc}} - {\rm{sw}}}}\left( {T(t)} \right) + {\delta ^{18}}{{\rm{O}}_{{\rm{sw}}}}\left( t \right),$$ (22)

where r represents the distance to the solid centre of the spherical foraminifera calcite crystals and \({r_0}\) represents their radii. Equation (20) was numerically solved using the implicit finite differences method40 with M cells numbered from the solid/fluid interface to the solid centre. The numerical discretisation leads to the following mass balance equation:

$$\begin{array}{ccccc} \frac{{ - {S_{i,i - 1}}D\left( z \right)}}{{{\rm{d}}r}}{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_{i - 1}^{n + 1} + \left( {\frac{{{V_i}}}{{\Delta t}} + \frac{{{S_{i,i - 1}}D\left( z \right)}}{{{\rm{d}}r}} + \frac{{{S_{i,i + 1}}D\left( z \right)}}{{{\rm{d}}r}}} \right){{\rm{\delta }}^{18}}{O_{{\rm{cc}}}}_i^{n + 1} \\ - \frac{{{S_{i,i + 1}}D\left( z \right)}}{{{\rm{d}}r}}{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_{i + 1}^{n + 1} = \frac{{{V_i}}}{{\Delta t}}{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_i^n,\end{array}$$ (23)

where \({S_{i,i - 1}} = 4\pi {\left( {r - {\rm{d}}r{\rm{/}}2} \right)^2}\) (respectively \({S_{i,i + 1}}\)) is the surface area between the ith and (i−1)th (respectively the ith and (i+1)th) spherical shells of calcite with thicknesses of dr, located at a distance r–dr from the centre. \({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_{i - 1}^{n + 1}\), \({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_i^{n + 1}\) and \({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_{i + 1}^{n + 1}\) represent the 18O concentrations within the (i−1)th, ith and (i+1)th shells at time step n+1; and \({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_i^n\) is the 18O concentration within the ith shell at the time step n. \({V_i}\) is the volume of the ith shell defined as \(\frac{4}{3}\pi \left( {{{\left( {r + {\rm{d}}r{\rm{/}}2} \right)}^3} - {{\left( {r - {\rm{d}}r{\rm{/}}2} \right)}^3}} \right)\), Δt is the time interval and D(z) is the diffusion coefficient at a given burial depth z.

Applying the appropriate boundary conditions permits the calculation of the mass balance for the cell of the mesh in contact with the pore water:

$$\begin{array}{ccccc}{V_1}\frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_1^{n + 1} - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_1^n}}{{\Delta t}} = {S_0}D\left( z \right)\frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_{in}^{n + 1} - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_1^{n + 1}}}{{{\rm{d}}r}} \hfill \\ + {S_{1,2}}D\left( z \right)\frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_2^{n + 1} - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_1^{n + 1}}}{{{\rm{d}}r}}, \end{array}$$ (24)

with S 0 representing the surface area of the calcite crystals of radius r 0 .

Similarly, the mass balance can be calculated for the Mth cell corresponding to the centre of a calcite sphere by supposing that the inner interface is impermeable, which sets a no flux boundary condition:

$${V_M}\frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_M^{n + 1} - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_M^n}}{{\Delta t}} = {S_{M,M - 1}}D\left( z \right)\frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_{M - 1}^{n + 1} - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_M^{n + 1}}}{{{\rm{d}}r}}.$$ (25)

Finally, the overall present-day isotope composition of the foraminifera (\({\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right)\)) is calculated via the cumulative isotope composition of calcite on each cell in the domain (t age being the time needed to reach present day):

$${\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right) = \frac{{\mathop {\sum}\limits_{i = 1}^M {{V_i}{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}_{_i}({t_{{\rm{age}}}})} }}{{\left( {4{\rm{/}}3} \right)\pi r_0^3}}.$$ (26)

Input parameters of the simulations

The main parameters of the present simulations were the seawater temperature, the O isotope composition of seawater, the sediment burial rate, the geothermal gradient and the activation of oxygen diffusion in foraminifera calcite (as demonstrated below, the present simulations do not depend on the intrinsic diffusion coefficient).

Here, benthic foraminifera were assumed to have grown in equilibrium with the contemporary deep seawater at a constant temperature of 3.5 °C, and the planktonic foraminifera were assumed to have grown in equilibrium with the surface seawaters at −2, 2, 7, 13, 20, 25, 27 and 27.5 °C, roughly corresponding to the surface waters of the modern ocean at 70, 60, 50, 40, 30, 20, 10 and 0° latitude10. In other words, we assumed that the ocean temperatures remained constant during the late Cretaceous and the entire Cenozoic.

The O isotope composition of seawater has varied through geologic times, primarily as a consequence of ice sheet growth and decay6,7,8,9. The calculation of the influence of the global ice volume on the average seawater δ 18O depends on how much continental ice was on the planet at a given time. Zachos et al.7 suggested that the growth of the ice sheets on Antarctica and in the Northern Hemisphere have been responsible for a total decline of ∼2.3‰ of the seawater δ 18O values (1.2‰ for Antarctic ice sheets and 1.1‰ for Northern Hemisphere ice sheets). However, a value of −1.2‰ for an ice-free ocean is commonly encountered in the literature2. Models of present-day ice sheet growth have indicated that an ice-free ocean would have an average δ 18O value between −0.89 and −1.1‰ (VSMOW)41,42. Assuming a more or less ice-free ocean over the late Cretaceous and the entire Cenozoic, the δ 18O of the seawater was fixed to −1‰ (VSMOW) for the present simulations. The sediment pore water in which the foraminifera have been fossilised was assumed to be isotopically similar to the seawater in which they lived32 (as supported by the δ 18O close to −1‰ (VSMOW) of the pore water in the old sediments).

To properly estimate the burial rate and consider the inevitable sediment compaction, the burial depth of the sediments from the DSDP/ODP/IODP sites used in the Friedrich et al.9 compilation of the deep sea benthic δ 18O values have been plotted against their ages, excluding sites that have undergone significant erosion. A sliding five-point average of these data is reported in Fig. 3a with a reasonable fit with z max = 500 m, α = 10 and β = 1.5:

$$z = {z_{{\rm{max}}}}{\left( {1 - \exp \left( { - \frac{t}{\alpha }} \right)} \right)^\beta }.$$ (27)

The geothermal gradients in oceanic sediments have been estimated to range between 30 and 70 °C km−1 32. Because the temperature experienced by the fossil foraminifera directly controls the extent of solid-state diffusion, the present model was computed using different geothermal gradients, ranging from 40 to 60 °C km−1. The temperature experienced by the sediments θ (in °C) is the product of the burial depth z (in metres) and the geothermal gradient \(

abla \theta \) (in °C m−1):

$$\theta =

abla \theta * z{\left( {1 - \exp \left( { - \frac{t}{\alpha }} \right)} \right)^\beta }.$$ (28)

As detailed above, Ea foram values ranging from 81 to 94 kJ mol−1 had to be considered. The present numerical simulations were thus conducted using different activation energies ranging from 85 to 95 kJ mol−1.

Non-dependence on the intrinsic diffusion coefficient

The present numerical simulations do not depend on the intrinsic diffusion coefficient (D 0,foram ), as demonstrated here. If calcite crystals are approximated by simple rods, the Q ratio corresponds to:

$$Q = \frac{d}{{{r_0}}} = \frac{2}{{{r_0}\sqrt \pi }}\sqrt {{D_{{\rm{foram}}}}{t_{{\rm{xp}}}}} ,$$ (29)

where r 0 is the radius of the calcite domain. Combining Eq. (29) with Eqs. (2) and (13) yields the following:

$${r_0} = \frac{2}{{Q\sqrt \pi }}\sqrt {{D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right){t_{{\rm{xp}}}}} .$$ (30)

As emphasised here, the integral of any diffusion profile can be approximated via a step function with length d and height C 0 . Therefore, at any time, the instantaneous isotope composition of the calcites (\({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( t \right)\)) can be expressed using the following mass balance:

$${\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( t \right) = \frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r = {r_0},t} \right)d + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,t = 0} \right)\left( {{r_0} - d} \right)}}{{{r_0}}},$$ (31)

which yields, after rearrangement,

$${\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( t \right) = \frac{d}{{{r_0}}}\left( {{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},t} \right) - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right)} \right) + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right).$$ (32)

Therefore, after a given duration t age , the isotope composition of the foram (\({\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right)\) ) as previously defined is as follows:

$${\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right) = \frac{1}{{{t_{{\rm{age}}}}}}\frac{1}{{{r_0}}}\mathop {\int }

olimits_0^{{t_{{\rm{age}}}}} d(\tau )\left( {{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},\tau } \right) - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right)} \right){\rm{d}}\tau + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right).$$ (33)

Replacing r 0 and d(τ) with the values taken from Eqs (2), (13) and Eq. (30) yields the following:

$$\begin{array}{l}{\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right) = \frac{1}{{{t_{{\rm{age}}}}}}\frac{1}{{\frac{2}{{Q\sqrt \pi }}\sqrt {{D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{xp}}}}} \right){t_{{\rm{xp}}}}} }}\\ \mathop {\int }

olimits_0^{{t_{{\rm{age}}}}} \frac{2}{{\sqrt \pi }}\sqrt {{D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{xp}}}}} \right)} \left( {{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},\tau } \right) - } \right.\\ \left. {{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right)} \right){\rm{d}}\tau + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right),\end{array}$$ (34)

which yields, after rearrangement,

$$\begin{array}{ccccc}{\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right) = \frac{1}{{{t_{{\rm{age}}}}}}\frac{1}{{\frac{1}{Q}\sqrt {{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right){t_{{\rm{xp}}}}} }}\hfill \\ \mathop {\int }

olimits_0^{{t_{{\rm{age}}}}} \sqrt {{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{RT(\tau )}}} \right)} \left( {{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},\tau } \right) - {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right)} \right){\rm{d}}\tau \\ + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right).\end{array}$$ (35)

Thus, \({\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right)\) does not depend on D 0,foram .

This result holds true in the case of spherical calcites, at least when the characteristic length of diffusion is negligible compared to the calcite grain size, which is the case for most of the present simulations. In fact, rewriting Eq. (13) and Eq. (15) yields the following:

$$d = \frac{2}{{\sqrt \pi }}\sqrt {{D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{RT(t)}}} \right)t} = \sqrt {{D_{0,{\rm{foram}}}}} f(t)$$ (36)

$${\rm{and}}\quad \quad \quad {r_0} = \frac{2}{{\sqrt \pi }}\frac{{\sqrt {{D_{0,{\rm{foram}}}}{\rm{exp}}\left( { - \frac{{{\rm{E}}{{\rm{a}}_{{\rm{foram}}}}}}{{R{T_{{\rm{xp}}}}}}} \right){t_{{\rm{xp}}}}} }}{{1 - {{\left( {1 - Q} \right)}^{1/3}}}} = \sqrt {{D_{0,{\rm{foram}}}}} {\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}},$$ (37)

in which all the constant parameters of Eq. (15) but \(\sqrt {{D_{0,{\rm{foram}}}}} \) are included in Cte xp .

If the characteristic length of diffusion is negligible compared to the calcite grain size, the instantaneous isotope composition of calcite (\({\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( t \right)\)) can be written as follows:

$$\begin{array}{ccccc}{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( t \right) = \frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},t} \right){V_{{\rm{shell}}}} + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right){V_{{\rm{int}}}}}}{{{V_0}}} \hfill\\ = \frac{{{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},t} \right)4\pi {{({r_0} - d)}^2}d + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right)\frac{{4\pi }}{3}{{({r_0} - d)}^3}}}{{\frac{{4\pi }}{3}{r_0}^3}}, \end{array}$$ (38)

where V shell stands for a fully re-equilibrated external shell of calcite crystals, V int is the internal volume unaffected by the re-equilibration and V 0 is the initial volume.

Replacing d and r 0 by their values yields the following:

$${{\delta ^{18}}{{{\rm {O}}} _{{\rm{cc}}}}\left( t \right) = \frac{\begin{array}{ccccc} {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},t} \right)4\pi {\left( {\sqrt {{D_{0,{\rm{foram}}}}} {\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}} - \sqrt {{D_{0,{\rm{foram}}}}} f(t)} \right)^2} \sqrt {{D_{0,{\rm{foram}}}}} f(t) + {\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right)\frac{{4\pi }}{3}{\left( {\sqrt {{D_{0,{\rm{foram}}}}} {\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}} - \sqrt {{D_{0,{\rm{foram}}}}} f(t)} \right)^3}\\ \end{array}}{{\frac{{4\pi }}{3}{{\left( {\sqrt {{D_{0,{\rm{foram}}}}} {\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}}} \right)}^3}}},}$$ (39)

which yields, after rearrangement,

$${\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( t \right) = \frac{{\left( {{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},t} \right){{\left( {{\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}} - f(t)} \right)}^2}f(t) + \frac{1}{3}{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right){{\left( {{\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}} - f(t)} \right)}^3}} \right)}}{{\frac{1}{3}{\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}}^3}}.$$ (40)

After a given duration t age , the isotope composition of a foraminifera (\({\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right)\)) is written as follows:

$$\begin{array}{ccccc} {\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right) = \frac{3}{{{t_{{\rm{age}}}}{\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}}^3}}\mathop {\int }

olimits_0^{{t_{{\rm{age}}}}} \left( {{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {{r_0},\tau } \right){{\left( {{\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}} - f(t)} \right)}^2}f(t)} \right. \\ \left. { + \frac{1}{3}{\delta ^{18}}{{\rm{O}}_{{\rm{cc}}}}\left( {r,0} \right){{\left( {{\rm{Ct}}{{\rm{e}}_{{\rm{xp}}}} - f(t)} \right)}^3}} \right){\rm{d}}\tau .\end{array}$$ (41)

Thus, \({\delta ^{18}}{\rm{O}}_{{\rm{cc}}}^{{\rm{overall}}}\left( {{t_{{\rm{age}}}}} \right)\) does not depend on D 0,foram in the case of spherical calcites.

Data availability

Original data and numerical codes are available from the corresponding author upon request.