FAME (v1.0): a simple module to simulate the effect of planktonic foraminifer species-specific habitat on their oxygen isotopic content

The oxygen-18 to oxygen-16 ratio recorded in fossil planktonic foraminifer shells has been used for over 50 years in many geoscience applications. However, different planktonic foraminifer species generally yield distinct signals, as a consequence of their specific living habitats in the water column and along the year. This complexity is usually not taken into account in model–data integration studies. To overcome this shortcoming, we developed the Foraminifers As Modeled Entities (FAME) module. The module predicts the presence or absence of commonly used planktonic foraminifers and their oxygen-18 values. It is only forced by hydrographic data and uses a very limited number of parameters, almost all derived from culture experiments. FAME performance is evaluated using the Multiproxy Approach for the Reconstruction of the Glacial Ocean surface (MARGO) Late Holocene planktonic foraminifer calcite oxygen-18 and abundance datasets. The application of FAME to a simple cooling scenario demonstrates its utility to predict changes in planktonic foraminifer oxygen-18 to oxygen-16 ratio in response to changing climatic conditions.

1 Introduction

Since the early work of Emiliani (1955), oxygen-18 isotopic abundance in calcite fossil foraminifer tests recovered from oceanic sediments has been widely used to reconstruct the past variations in oxygen-18 content of seawater as well as its temperature, the two main variables that affect the ratios of oxygen-18 to oxygen-16 in calcite (hereafter R¹⁸O_c). The recognition that different species of foraminifers from the same sediment core yielded different R¹⁸O_c was made early on (e.g., Duplessy et al., 1970; Lidz et al., 1968; Berger, 1969; Fairbanks and Wiebe, 1980; Deuser, 1987), though it was an attempt by Emiliani (1954) to relate depth habitat of foraminifers to the density of seawater that led to the revelation that the R¹⁸O_c recorded by fossil foraminifers likely reflected the average depth habitat of individual species. Through in situ water column sampling via opening–closing plankton nets, Jones (1967) through faunal abundance counts corroborated the depth habitats that Emiliani (1954) inferred through isotopic analysis. However, increased plankton sampling (Bé and Tolderlund, 1971) and the advent of the sediment trap have shown that different species have different living habitats in the water column and throughout the year and that in some cases the foraminiferal R¹⁸O_c presents an offset with respect to the equilibrium calcite oxygen-18 to oxygen-18 ratio (Mix, 1987; Bijma and Hemleben, 1994; Ortiz et al., 1995; Kohfeld et al., 1996; Bauch et al., 1997; Schiebel et al., 2002; Simstich et al., 2003; Mortyn and Charles, 2003; Rebotim et al., 2017; Jonkers and Kucera, 2015). This complexity is usually not accounted for in paleoceanographic studies. Instead, the approximation is often made that each planktonic foraminifer species has an apparent living depth – defined as the water depth where equilibrium calcite formation would approximate their measured calcite R¹⁸O_c value in the water column – that can vary by hundreds of meters from one region to another. To correctly interpret the wealth of information coming from the calcite oxygen-18 to oxygen-16 ratio record, especially when multiple species are measured at the same geographical location, there is a need to take into account the impact of depth habitat and growth season on each species' calcite oxygen-18. Minor contributors to the resultant calcite R¹⁸O_c, such as carbonate ion concentration (Spero et al., 1997; Pearson, 2012) and symbionts (e.g., Ezard et al., 2015; Spero, 1998), may modulate the absolute values, making species-specific comparisons problematic; however, their overall contribution may also covary and/or autocorrelate with temperature and latitudinal gradients; therefore, this paper focuses on the major components only.

In recent years, the development of water-isotope-enabled ocean models has allowed the simulation of the two variables at the root of the calcite R¹⁸O_c record: seawater temperature and oxygen-18 abundance. So far, water-isotope-based model–data studies have generally compared planktonic oxygen-18 ratios to equilibrium calcite R¹⁸O_c values. The equilibrium calcite ratio in that case is computed from annual averaged seawater temperatures and oxygen-18 abundance taken either at surface or averaged over the upper 50 to 100 m of the water column (e.g., Caley et al., 2014; Werner et al., 2016). To go one step further, it is necessary to account for species-specific habitat when computing calcite R¹⁸O_c. The Foraminifers As Modeled Entities (FAME) approach is underpinned by two arguably simple principles: (1) the weighting due to species habitat is reflected in the calcite oxygen-18 ratio record, and (2) the model derived should be kept simple to allow its offline application to the output of climatic models without the need of rerunning the entire climate model simulations.

After having developed the FAME methodology, we found out that the idea was already present in a theoretical framework in Mix (1987) and in one following study (Mulitza et al., 1997), in the latter referred to as Mix's model. The most notable difference between the early study of Mix (1987) and the present one is the actual definition of the weighting functions. Mix (1987) assumed them to be simple Gaussians, whereas we build ours on the laboratory culture-based equations of planktonic foraminifer growth rates as a function of temperature given in Lombard et al. (2009).

Since the early work of Mix (1987), other methods were developed to approach the species-specific complexity of planktonic foraminifers. Schmidt (1999) developed a simple module to compute planktonic foraminifer R¹⁸O_c in a water-isotope-enabled global ocean model. However, in his approach, water depths at which planktonic foraminifers calcify and their seasonal growth patterns are fixed for each species. Therefore, such a module cannot properly account for the impact of climatic changes on foraminifer living conditions. Fraile et al. (2008) and Lombard et al. (2011) developed models predicting the abundance of common planktonic foraminifer species in response to hydrographic data and food concentration. Both these models predict the relative abundances of the different simulated foraminifer species, an information which is not needed to assess individual species' oxygen-18 but entails a large number of empirical parameters, i.e., 21 and 15 parameters per planktonic foraminifer species in Fraile et al. (2008) and Lombard et al. (2011), respectively. Moreover, Fraile et al. (2008) derived the sensitivity of each species with respect to temperature from sediment-trap data, so that their model can only account for changes in seasonality and not in depth habitat. In contrast, the FORAMCLIM model (Lombard et al., 2011) predicts both season and water depth of each species' potential maximum abundance. In fact, FAME can be viewed as a simplified version of FORAMCLIM (only retaining FORAMCLIM computation of growth rates as a function of temperature), expanded by a mechanistic calculation of species-dependent calcite oxygen-18. FAME is only forced by hydrographic data and only uses six parameters per planktonic foraminifer species that are all derived from culture experiments, plus one parameter accounting for the effect of the accretion of a calcite crust by N. pachyderma. Taken together, these characteristics make FAME a uniquely simple and robust model designed to predict changes in the oxygen-18 to oxygen-16 ratio of commonly used planktonic foraminifers in response to changing climatic conditions.

2 Methodology

The calcite oxygen-18 ∕ oxygen-16 ratio (reported in δ¹⁸O_c, in permil versus Vienna Pee Dee belemnite – V-PDB – in what follows) of planktonic foraminifers is intrinsically a four-dimensional signal, acquired at a specific season (time dimension), over a specific depth range and area in the ocean (space dimensions). The mean δ¹⁸O_c signal measured on a sample composed of a number of individual foraminifer shells of one species is thus the integration of many different single δ¹⁸O_c paths in this four-dimensional space. If we suppose that the sampled population is representative of the living conditions of the species, it is thus likely that there is an oversampling of the areas and time representing favorable growth conditions and an undersampling of area and time with unfavorable growth conditions. Hence, a reasonable way to predict the mean δ¹⁸O_c of a foraminifer sample constituted of a number of individuals is to weight the oceanic conditions by the growth rate of each individual. The predicted δ¹⁸O_c is then a weighted sum of these conditions in space and time.

2.1 Basic equations

To define the effect of the habitats of the different foraminifer species, we first consider a subset of the growth functions derived by Lombard et al. (2009) from culture experiments (Fig. 1) following the original formulation of Kooijman (2000). For each foraminifer species k considered, the growth function is written as

$$\begin{array}{}\text{(1)}& \mathit{\mu }\left(T,k\right)={\displaystyle \frac{\mathit{\mu }\left(T_{\mathrm{1}},k\right)\cdot \exp \left(\frac{T_{\mathrm{A}}}{T_{\mathrm{1}}}-\frac{T_{\mathrm{A}}}{T}\right)}{\mathrm{1}+\exp \left(\frac{T_{\mathrm{AL}}\left(k\right)}{T}-\frac{T_{\mathrm{AL}}\left(k\right)}{T_{\mathrm{L}}\left(k\right)}\right)+\exp \left(\frac{T_{\mathrm{AH}}\left(k\right)}{T_{\mathrm{H}}\left(k\right)}-\frac{T_{\mathrm{AH}}\left(k\right)}{T}\right)}},\end{array}$$

where μ(T,k) is the growth rate at temperature T for the species k, μ(T₁,k) is the growth rate for a chosen reference temperature T₁ (20 ^∘C or 293 K), T_A is the Arrhenius temperature, T_H(k) and T_L(k) define the upper and lower boundaries of the growth tolerance range for the species k, and T_AH(k) and T_AL(k) are the Arrhenius temperatures for the decrease in growth rate, respectively, above and below these boundaries for species k (Lombard et al., 2009). In the present study, we use the nominal values of Eq. (1) parameters given in Lombard et al. (2009) with the exception of T_L for G. bulloides. Indeed, comparing the output of FAME with sediment-trap data from the subpolar North Atlantic (Jonkers et al., 2013) showed that the nominal value of T_L=281.1 K was likely too high, causing an absence of growth outside of the 3 summer months. In contrast, subpolar North Atlantic sediment-trap data indicate that, on average over the 4 years of observations, significant G. bulloides fluxes prevailed from the end of June to the middle of November. We hence chose a value of T_L closest to the nominal value of Lombard et al. (2009) that would allow the extension of the growing season into the fall in agreement with the data pattern. Hence, a value of T_L=280 K was used for G. bulloides within FAME.

Figure 1 Growth functions corresponding to Eq. (1) over the full temperature range considered, replotted from Lombard et al. (2009).

We compute the μ(T,k) coefficient for all values of $T(x,y,z,t)$ in the World Ocean, T being a 4-D variable of space and time. This, in turn, gives us the growth rate of the different foraminifer species considered in a four-dimensional space as

$$\begin{array}{}\text{(2)}& \mathit{\mu }\left(T,k\right)=\mathit{\mu }\left(T\left(x,y,z,t\right),k\right).\end{array}$$

To avoid numerical issues in the code, we limit the value of μ(T,k) on the low end as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{llll}{\mathit{\mu }}^{\prime }\left(T\left(x,y,z,t\right),k\right)& =& \mathit{\mu }\left(T,k\right)& \text{if }\mathit{\mu }\left(T,k\right)\ge \mathrm{0.1}\\ & & & \cdot \underset{T}{max}\mathit{\mu }\left(T,k\right)\\ & =& \mathrm{0}& \text{otherwise.}\end{array}\end{array}$$

Given a four-dimensional input field for oceanic temperatures and δ¹⁸O of seawater, the equilibrium inorganic calcite δ¹⁸O value can be computed from the temperature equation of Kim and O'Neil (1997). Here, we use the quadratic approximation of that equation given in Bemis et al. (1998):

$$\begin{array}{}\text{(4)}& T=T_{\mathrm{0}}-b\cdot \left({{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sw}}\right)+a\cdot {\left({{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sw}}\right)}^{\mathrm{2}},\end{array}$$

where T₀=16.1 ^∘C, $b=-\mathrm{4.64}$ and a=0.09, δ¹⁸O_sw is the seawater δ¹⁸O. Since the seawater temperature (^∘C) and δ¹⁸O_sw (permil) are inputs, we can solve this equation to determine the value of δ¹⁸O_c. With the discriminant of the second-degree equation being

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }=b^{\mathrm{2}}-\mathrm{4}a\cdot \left(T_{\mathrm{0}}-T_{\mathrm{sw}}\right),\end{array}$$

it becomes

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c},\mathrm{eq}}={\displaystyle \frac{-b-\sqrt{\mathrm{\Delta }}}{\mathrm{2}a}}+{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{s}}\mathrm{w}-\mathrm{0.27},\end{array}$$

where the constant, 0.27, correction (Hut, 1987) accounts for the difference in the reference scales of seawater (permil versus Vienna mean standard ocean water – V-SMOW) and calcite (permil versus V-PDB).

In previous studies, we and others (Caley et al., 2014; Werner et al., 2016) computed the above δ¹⁸O equilibrium value, averaged over time and the surface layer (typically the first 50 m) to compare model results and measured δ¹⁸O_c from planktonic foraminifers. In the following, we will refer to this method as the “old method”, written formally as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c},\mathrm{om}}\left(x,y\right)={\displaystyle \frac{\mathrm{1}}{n_t}}\sum _{t=\mathrm{1}}^{n_t}{\displaystyle \frac{\mathrm{1}}{n_z}}\sum _{z=\mathrm{0}}^{z_b}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c},\mathrm{eq}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sw}}\left(x,y,z,t\right),T_{\mathrm{sw}}\left(x,y,z,t\right)\right),\end{array}$$

where n_t is the number of time steps, n_z the number of vertical levels and z_b the maximum depth.

The formalism used clearly expresses the fact that the old method is not species-specific nor season-specific since all time steps and vertical levels are averaged with the same weight. In contrast, the FAME method weighs the δ¹⁸O_c both in time and in the horizontal and vertical space according to the population abundances using the foraminifer growth formula (1). We thus write

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c},\mathrm{fm}}\left(x,y,k\right)={\displaystyle \frac{\mathrm{1}}{n_t}}\sum _{t=\mathrm{1}}^{n_t}{\displaystyle \frac{\mathrm{1}}{n_z}}\sum _{z=\mathrm{0}}^{z_b\left(k\right)}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c},\mathrm{eq}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sw}}\left(x,y,z,t\right),T_{\mathrm{sw}}\left(x,y,z,t\right)\right){\mathit{\mu }}^{\prime }\left(T(x,y,z,t),k\right),\end{array}$$

where z_b(k) is dependent on the species and constrained by core-top data (see below).

Using this set of equations, for any given seawater temperature and δ¹⁸O provided as a four-dimensional field and a given species k, we compute this species' δ¹⁸O_c over x,y (latitude, longitude) coordinates.

It should be clearly understood that this approach is not able and does not attempt to determine the relative abundances of the different species. Instead FAME provides a simplified approach to compute the δ¹⁸O_c of a generic population of foraminifers if environmental conditions permit its growth. From a model–data perspective, this approach enables one to compute the calcite δ¹⁸O for a given species, were it to exist in the sedimentary record. Due to the limitations set by Eq. (3), no calcite isotopic content is computed if μ^′ is zero, and hence these areas will be masked out in the following.

2.2 Growth function uncertainties

In the original work of Lombard et al. (2009), 95 % confidence intervals were shown per species' growth functions, but no equation was given for them. In order to nonetheless estimate the bias introduced in FAME by using the given functions, we combined the different possible values for the parameters of the growth functions to obtain functions that are close to the 95 % confidence intervals mentioned for most species and larger than the 95 % confidence intervals for others. This results in an upper and lower growth function for each species, as shown in Fig. A1, where the original data points of Lombard et al. (2009) are also given.

2.3 Reference datasets

In an attempt to validate the FAME approach, we apply its methodology to reference datasets, close to present-day observations. The first necessary step is the computation of a reference δ¹⁸O_c field as obtained when forced by climatological data.

For seawater temperature, we use the World Ocean Atlas 2013 (WOA13) (Locarnini et al., 2013) data at a monthly resolution. Considering that there is no equivalent seawater oxygen-18 gridded dataset available in the World Ocean Atlas fields and that the existing NASA Goddard Institute for Space Studies (GISS) gridded dataset (LeGrande and Schmidt, 2006) presents large deviations with respect to the seawater oxygen-18 (δ¹⁸O_sw) raw data in numerous locations, we derived a δ¹⁸O_sw dataset based on seawater salinity to δ¹⁸O_sw relationships. This dataset is built in two steps: (a) derivation of regional δ¹⁸O_sw – salinity relationships from GISS δ¹⁸O_sw and salinity (Schmidt et al., 1999) clustered by oceanic regions, and (b) computation of a δ¹⁸O_sw field based on the World Ocean Atlas 2013 (Zweng et al., 2013) salinity fields. The resulting field is at the World Ocean Atlas spatial resolution and is used as reference seawater oxygen-18 in the following. Details on the derivation of the δ¹⁸O_sw dataset are given in Appendix A.

Ganssen and Kroon (2000); Anand et al. (2003); Kuroyanagi and Kawahata (2004); Farmer et al. (2007)Pak and Kennett (2002); Pak et al. (2004); Kuroyanagi and Kawahata (2004); von Langen et al. (2005)Ganssen and Kroon (2000); Anand et al. (2003); Kuroyanagi and Kawahata (2004); Farmer et al. (2007)Bauch et al. (2002); Ganssen and Kroon (2000); Pak and Kennett (2002); Pak et al. (2004); Kuroyanagi and Kawahata (2004)Kohfeld et al. (1996); Bauch et al. (1997, 2002); Simstich et al. (2003); Kuroyanagi and Kawahata (2004)

Table 1 Maximum depth per species as computed from the optimization procedure. z_b is the depth yielding the smallest difference to the MARGO Late Holocene data (Waelbroeck et al., 2005). We computed a confidence interval $\left[{\mathit{\sigma }}_{z_d}^{\uparrow },{\mathit{\sigma }}_{z_d}^{\downarrow }\right]$ corresponding to a change of ±0.1 permil in the mean error. The ∞ sign indicates that no value of z_b within the range $\left[\mathrm{0},-\mathrm{1500}\right]$ yields the desired ±0.1 permil change.

¹ An encrustation term of +0.1 permil is taken into account in the case of N. pachyderma (see text).

As an independent test of the FAME results, we use the planktonic δ¹⁸O_c measurements from the Multiproxy Approach for the Reconstruction of the Glacial Ocean surface (MARGO) Late Holocene dataset (Waelbroeck et al., 2005) restricted to high chronozone quality levels (i.e., levels 1 to 4). A few errors have been corrected in the published dataset: these concern the suppression of 10 Neogloboquadrina incompta (or N. pachyderma right) data points from the Nordic Seas where only Neogloboquadrina pachyderma should have been listed, and one outlier N. pachyderma value with no age control that was erroneously listed as having a level-4 chronozone quality. The corrected version of MARGO Late Holocene planktonic oxygen-18 dataset is available in the Supplement.

As a result, the dataset used in the present study contains 248 values for Neogloboquadrina pachyderma, 128 values for Globigerina bulloides, 59 values for Neogloboquadrina incompta, 135 values for Globigerinoides ruber and 51 values for Globigerinoides sacculifer. In the remainder of the paper and following the genus reassignment of Spezzaferri et al. (2015), we will refer to the latter as Trilobatus sacculifer.

2.4 Calculation of the best-fitting maximum depth per foraminifer species

In Eq. (8), the maximum depth of integration per foraminifer species, z_b(k), is a free parameter and needs to be determined. We have chosen to calculate it as the depth where the δ¹⁸O_c simulated by FAME driven by the World Ocean Atlas 2013 temperature and derived seawater δ¹⁸O datasets is closest on average to MARGO Late Holocene δ¹⁸O_c data. The rationale behind this choice is to specifically design FAME to enable model–data comparison with isotopic records from marine sediment cores.

To determine the optimal value of z_b(k), we repeated successive runs of FAME with values of z_b ranging from 1500 m to the surface along the standard World Ocean Atlas vertical grid. The only difference between the different species at this stage are the species-specific terms in the equations presented and the data of each species from the MARGO Late Holocene set. The results obtained through this optimization procedure are given in Table 1. The maximum depths of calcification derived this way are remarkably close to what is known from the ecology of G. ruber, N. incompta, T. sacculifer and G. bulloides (Berger, 1969; Bijma and Hemleben, 1994; Ortiz et al., 1995; Schiebel et al., 2002; Mortyn and Charles, 2003; Rebotim et al., 2017). Only in the case of N. pachyderma, the computed value of z_b was much too deep (900 m) with respect to what studies based on opening–closing plankton nets show. Also, plankton haul studies have revealed that, whereas N. pachyderma seems to grow at relatively shallow depth, i.e., where the chlorophyll maximum is found, a calcite crust is added between ∼50 and 250 m, which greatly increases its mass (Kohfeld et al., 1996; Simstich et al., 2003). As a consequence, the δ¹⁸O of N. pachyderma collected in deep sediment traps and in surface sediment is systematically heavier than that of living non-encrusted N. pachyderma. To account for this effect, we have added a 0.1 permil “encrustation term” to our calculation of N. pachyderma calcite δ¹⁸O weighted by that species' culture-based growth rates. The encrustation term value has been chosen in order to simulate maximum depths in agreement with the literature. N. pachyderma simulated depth of maximum growth (Fig. 4e and Table 1) does indeed match very well the available observations. For instance, Fig. 4e shows a deepening of N. pachyderma depth of maximum growth from 0 to 30 m in the Greenland Sea to 100–350 m in the Norwegian Sea, in agreement with the apparent calcification depths reconstructed by Simstich et al. (2003).

Concerning T. sacculifer, although this species bears symbionts, and thus lives in the photic zone like G. ruber, it is known to produce calcite with higher δ¹⁸O values than G. ruber. These heavier δ¹⁸O values are thought to result from the accretion of gametogenetic calcite (for a certain unknown fraction of the shell mass) or from the precipitation of its final sac-like chamber deeper in the water column (Duplessy et al., 1981). This characteristic explains the deeper habitat depth computed for T. sacculifer versus G. ruber (maximum calcification depth estimates range from −200 to −75 m, best estimate of −100 m). Note that a deeper habitat for T. sacculifer than G. ruber is in agreement with observations (Table 1).

2.5 Evaluation of the model performance

2.5.1 Error distribution

Since the depth parameter was constrained using the MARGO Late Holocene dataset by error minimization, it is not surprising that the errors obtained with FAME are very small on average for each species considered (Fig. 2). The error distribution obtained with FAME is very similar to the one obtained with the simple surface equilibrium assumption for the two species closest to the surface (G. ruber and N. incompta). For deeper dwellers (T. sacculifer, G. bulloides & N. pachyderma), FAME results are better than those obtained with the old method, as expected, since deeper layers in the ocean are accounted for.

Figure 2 Error distribution for the “old method” (grey) and the “FAME method” (orange) using climatological datasets as compared to MARGO Late Holocene dataset (Waelbroeck et al., 2005). Best-fitting distributions are calculated and plotted as a solid line for the FAME method and as a dashed line for the old method, except for T. sacculifer, for which the small number of available data points yields a poor fit both for FAME and the old method. The mean and deviation are given for FAME and the old method at the top of each panels.

2.5.2 Robustness of results

To test the robustness of our calibration in depth or error distribution, we performed a full set of additional analyses using the lower and upper growth functions as presented in Fig. A1, introduced here above.

Regarding depth calibration, we find our results to be largely insensitive to the use of these upper-bound and lower-bound values for the growth functions. Specifically, the uncertainty in the maximum growth depth is largest on N. pachyderma (range of 475 to 600 m) and G. bulloides (400 to 450 m). It is somewhat smaller for T. sacculifer (100 to 125 m) and N. incompta (60 to 65 m). There is no impact for G. ruber.

Another method to check the impact of the uncertainty in the growth functions on our δ¹⁸O_c results is by keeping maximum computed depths constant and looking at the impact of the growth function on the mean difference between simulated and MARGO δ¹⁸O_c values shown in Fig. 2. When doing so (not shown), the resulting change is lower than 0.1 permil for all individual species. It therefore shows that our results are very robust and largely insensitive to the errors arising in the growth functions used (Lombard et al., 2009).

Figure 3 Model–data comparison of species' abundances. Ocean regions where FAME predicts that the species is present at some time of the year (${\mathit{\mu }}^{\prime }>\mathrm{0}$) are plotted in blue, with shades of blue indicating the number of months of presence. Overlaid are the MARGO Late Holocene data (quality levels 1–4) species' abundance data, plotted using the yellow–white to dark red color bar and given in percent. A qualitative correspondence between simulated FAME presence/absence and the occurrence of 10 % level in the difference species is noted.

2.5.3 Geographical distribution

To further check our methodology against the MARGO Late Holocene dataset, we compare the zone of presence of each species predicted by FAME (grossly determined by μ^′) with the observed reported abundances in the MARGO dataset (restricted to chronology quality values 1–4). As noted above, we cannot predict the relative abundance of each species. However, the method determines the species' absence or presence.

Figure 4 Depth of maximum growth for the species considered for the month of July. The color scale shows the depth in meters. Oceanic areas left in white correspond to areas where growth rates are below the threshold defined in Eq. (3).

The results presented in Fig. 3 show that, despite the exceptional simplicity of our approach, FAME predicts relatively well the spatial limits of the area occupied by each species. The two species whose presence distribution is best predicted are again G. bulloides and N. pachyderma, both showing a quite remarkable model–data match of the transition zones from presence to absence. N. incompta and G. ruber also show quite satisfactory results, with only a few outliers in specific areas: FAME computes overextended coverage of N. incompta in the Gulf of Guinea and of G. ruber along the coast of Namibia.

The computed spatial coverage of T. sacculifer is slightly too extended towards high northern and possibly southern latitudes. The very low number of high-quality dated data points in the latter area prevents a definitive conclusion. Also, specific zones, consistent for several species, may be noted, such as the Benguela upwelling regions where FAME fails to predict the absence of T. sacculifer and G. ruber.

Figure 5 Comparison of the depth of maximum growth for N. pachyderma and G. bulloides for January and July. Color scale is as in Fig. 3.

One possible explanation for this mismatch could be the impact of increased nutrient availability on observed abundances as a consequence of the upwelling systems, whereas nutrients are at present ignored in the FAME approach. Another possibility could be the quality of the vertical oceanic structure obtained from the World Ocean Atlas in those upwelling regions. Finally, it should be noted that our comparison ignores the natural interannual variability since we are using climatologies. The interannual variability involves changes in the location of the fronts and currents, and thus bears the potential of shifting the spatial boundaries between the different foraminifer species.

Further discussion of the abundance comparison including all data points from the MARGO Late Holocene dataset regardless of the dating quality is given in Appendix B.

To further investigate the functioning of the FAME model, it is useful to consider the spatial distribution of the depth at which each species' growth is maximum. An example is given for the month of July in Fig. 4. It clearly shows that even though the maximum depth allowed for each species is fixed through the z_b(k) parameter, the predicted/computed calcification depth varies according to the location in the World Ocean. Except for G. ruber which always calcifies in the topmost ocean layers, the depth of maximum growth exhibits large spatial variations, notably at the edge of the species' domains; in July, this is particularly marked in the case of G. bulloides and N. pachyderma (Fig. 4d and e).

Likewise, it is useful to consider the seasonal variations in the depth of maximum growth for a given species. We propose to highlight this aspect for the two species that show the largest variations: N. pachyderma and G. bulloides at two extreme months (January and July) (Fig. 5). For both species, the area of computed non-zero contribution varies along the year, with an expansion (reduction) of the area occupied by N. pachyderma in the Northern Hemisphere in January (July), while the regions occupied by G. bulloides shift towards higher (lower) latitudes in the Northern Hemisphere in July (January). These seasonal changes are a direct response of these species' preferred habitat to temperature. FAME thus mechanistically predicts the adaptation of planktonic foraminifer depth habitat to maintain optimal living conditions. For instance, Fig. 5b and d clearly show that G. bulloides is predicted to dwell deeper at low latitudes when surface temperature rises above its preferred temperature range. Similarly, Fig. 5b and d show that G. bulloides is present at higher northern latitudes in July than in January, so that the growing season actually tracks the species' preferred living conditions, as observed (Jonkers and Kucera, 2015).

2.5.4 Effect of a large climatic change on the computed oxygen-18 content of the calcite

Though FAME gives realistic results when forced by atlas data, it is mostly designed to retrieve the species-specific effect of climate change on the recorded δ¹⁸O_c. To highlight the effect of seasonal and vertical weighting of the δ¹⁸O_c signal computed by FAME, we have performed a simplified experiment showing the effect of a change in the foraminifers' living conditions on their δ¹⁸O_c signal.

To simulate a change in climatic conditions, we apply a homogeneous 4 ^∘C decrease to the WOA13 sea temperature dataset and compute the anomaly in δ¹⁸O_c between that new cold state and the original one for each species, as well as for the surface equilibrium approach (Fig. 6). This anomaly is noted Δ¹⁸O_c in what follows.

Applying a spatially homogeneous temperature change should result is a quasi-homogeneous temperature change in the equilibrium calcite Δ¹⁸O_c, following Eq. (6). It is indeed what is obtained in Fig. 6e, with Δ¹⁸O_c values between 0.8 in the tropics and 1 permil at high latitudes. When applying the FAME equations, we obtain large spatial variations in Δ¹⁸O_c with values down to −0.75 and up to 1 permil. All species share a common pattern of lower Δ¹⁸O_c at the border of their living domain and close to equilibrium values at the center of their living domain. More specifically, the smallest differences to the equilibrium are recorded by N. pachyderma and the largest, negative, differences are computed for G. bulloides. The species with the smallest vertical living range, G. ruber, has the most homogeneous distribution. The range of values (minimum to maximum) is always close to 1 permil with the exception of G. bulloides that presents a total range of 1.6 permil. This large range of Δ¹⁸O_c for G. bulloides is a consequence of its growth over a large range of temperatures (Eq. 1). In general, the maximum simulated Δ¹⁸O_c values are systematically 0.1 to 0.2 permil lower than the equilibrium value.

This simple scenario, though unrealistic with respect to actual climatic applications, shows the potential of FAME to unravel the climatic signal embedded in multispecies isotopic records and thus opens the door to transient climate–data intercomparison where the species' specific behavior is taken into account.

Figure 6 δ¹⁸O_c response to a horizontally and vertically homogeneous 4 ^∘C cooling applied to the WOA13 dataset, Δ¹⁸O_c. Results are expressed in permil for each species (a–d) and for the equilibrium surface calcite approach.

3 Summary and conclusions

We developed the FAME (Foraminifers As Modeled Entities) module to account for planktonic foraminifer species-specific habitat when computing their calcite oxygen-18. In contrast to models predicting the abundance of planktonic foraminifers, FAME only aims at predicting the presence or absence of a given species and its oxygen-18 value. FAME is only forced by hydrographic data and uses a very limited number of parameters, almost all derived from culture experiments. Taken together, these characteristics make FAME a uniquely simple and robust model for predicting changes in the oxygen-18 of commonly used planktonic foraminifer species in response to changing climatic conditions. FAME performance is evaluated using MARGO Late Holocene planktonic foraminifer δ¹⁸O_c and abundance datasets. We show that FAME predicts remarkably well the presence/absence of G. ruber, N. incompta, N. pachyderma and G. bulloides over most of the World Ocean, while yielding a slightly less good prediction of T. sacculifer presence/absence. Investigating the simulated seasonal pattern, we show that the predicted growing season and habitat depth track the species' preferred living conditions, as observed in plankton hauls and sediment-trap data. Finally, the application of FAME to a simple cooling scenario demonstrates that computed changes in species-specific δ¹⁸O_c are much more spatially variable than the computed change in equilibrium surface calcite. Coupling the FAME module to isotope-enabled climate models makes it possible for the first time to extract the climatic information contained in isotopic time series measured on different planktonic species at the same location. This opens the possibility to better reconstruct the evolution of the upper water column structure than ever before, notably over climate transitions.

Derivation of the a reference δ¹⁸O_sw dataset

We constructed our reference δ¹⁸O_sw dataset at World Ocean Atlas standard resolution (1^∘ grid) through a three-step methodology: (a) construction of an appropriate basin mask to allow clustering the GISS global dataset regionally, (b) derivation of δ¹⁸O_sw – salinity relationships for each of these basins and (c) use of these relationships to obtain a δ¹⁸O_sw field at WOA spatial resolution.

A1 Construction of the basin masks

Figure A1 Growth functions corresponding to Eq. (1) over the full temperature range considered, with added lower-range and upper-range curves (respectively, short and long dashed curves) chosen to mimic the 95 % confidence interval of Lombard et al. (2009). Data points are the original data obtained by Lombard et al. (2009).

Figure A2 Basin masks as defined in the WOA standard mask file (a) and in FAME (b) on the same 1^∘ resolution grid. Values correspond to the basins defined in Table A1.

Table A1 Comparative list of basin masks in WOA and FAME. The “value” field provides the integer value used in the netCDF file to specify the respective basin on the WOA grid.

The “*” symbol highlights the regions where FAME and WOA regions do not cover the same area (see text).

Our native resolution being the 1^∘ regular grid of the World Ocean Atlas, we first retrieved the available basin mask file on that grid from the National Oceanic and Atmospheric Administration (NOAA) website https://www.nodc.noaa.gov/OC5/woa13/masks13.html (last access: 5 October 2016) and converted it to a netCDF format file (http://www.unidata.ucar.edu/software/netcdf, last access: 30 August 2018). The basins defined in the WOA base mask did not perfectly fit our purpose; we hence modified the masks to isolate some particular regions where the δ¹⁸O_sw and salinities are specific (e.g., Sea of Okhotsk) or merge some regions of the WOA mask into larger ensembles (e.g., Hudson Bay). A summary of the basins in the original file and in ours is given in Table A1. The Pacific and Atlantic oceans were split into south, north and tropical parts, based on boundaries at 30^∘ north and south, respectively. The Indian Ocean has been only split in two: north and south using the 30^∘ south boundary. The Bay of Bengal has been kept a separated basin as in the original file. The Greenland–Iceland–Norwegian (GIN) seas were made a separated basin from the Arctic Ocean, using the boundaries at 80^∘ north and at 20^∘ east. We also extended the Hudson Bay mask area to include the Hudson Strait and Ungava Bay, since these do not represent the same water mass properties as the North Atlantic Ocean. The limit used is −64.5^∘ west, corresponding to the southern tip of Resolution Island on the grid given. Finally, the same procedure was applied to define the Sea of Okhotsk, using the official definition of the International Hydrographic Organization (IHO SP-23() IHO SP-23: Limits of Oceans and Seas, International Hydrographic Organization, Imp. Monégasque, Monte Carlo, special Publication no. 23, 1953). The results of this whole procedure are shown in Fig. A2. In Table A1, some values are annotated with a “*” to highlight basins having the same value in FAME as in the standard WOA but covering a different area: the Arctic Ocean from which the GIN seas have been taken out in FAME, Hudson Bay which covers a part of the former Atlantic basin of the WOA given its aforementioned expansion to the Hudson Strait and Ungava Bay.

The netCDF data file resulting from this procedure is provided in the Supplement to the paper.

Computation of the δ¹⁸O_sw – salinity relationships

The basins defined in the previous section are then used to cluster the raw data, δ¹⁸O_sw and salinity, of the GISS database (Schmidt et al., 1999) in the respective basins. Furthermore, only data locations where both δ¹⁸O_sw and salinity are given in the original database for depths less than 200 m are retained under the additional constraint that depth of the ocean should be more than 175 m. The latter are to ensure that the values are representative of high sea values and not coastal areas, possibly under fluvial influence. Additionally, all values below 5 permil in salinity are ignored for all basins. Lastly, for two basins (North Atlantic and Bay of Bengal), the existence of two different slopes where only one corresponds to open-ocean conditions renders necessary the addition of one additional condition to keep only the latter. We thus added a limit at 27 permil in salinity for those two basins.

The resulting slopes, intercept and correlation coefficients are given in the Table A2. Using those relationships, we further compute the δ¹⁸O_sw in the WOA geographical grid from the WOA salinity fields.

Table A2 Values obtained for the δ¹⁸O_sw – salinity relationships.

Appendix B: Further discussion of predicted and observed planktonic foraminifer abundances

Figure B1 Model–data comparison of species' abundances. Ocean regions where FAME predicts that the species is present at some time of the year (${\mathit{\mu }}^{\prime }>\mathrm{0}$) are plotted in blue, with shades of blue indicating the number of months of presence. Overlaid are the MARGO Late Holocene data (all quality levels) species' abundance data, plotted using the yellow–white to dark red color bar and given in percent. A qualitative correspondence between simulated FAME presence/absence and the occurrence of 10 % level in the difference species is noted.

In the main text, we have only compared the results of FAME to the data points in the MARGO Late Holocene database that were characterized by high chronological control quality. Since this drastically restricts the geographical extent covered by MARGO data, and in the interest of completeness, we propose here a short discussion based on all points of the MARGO Late Holocene database, regardless of their chronological control quality. The interest of Fig. B1 is to provide some information in the Southern, Pacific and Indian Ocean regions that are largely void of points in the previous figure. While the bulk of the conclusions given in the main text is unchanged by this new comparison, we may highlight the following.

The unsorted distribution for G. ruber is not very different from the one described above, albeit with a good definition of the Southern Ocean abundance limit where FAME results are in good accordance with MARGO. Also, one may note a series of points without the presence of G. ruber in the equatorial Pacific in MARGO, an aspect which is not predicted by FAME. However, these points are mingled with points with G. ruber presence in the MARGO database, indicating they could be an artifact resulting from the presence of older sedimentary material in the unsorted MARGO database; it is thus difficult to draw a firm conclusion.

Regarding N. incompta, the picture is pretty much the same as described in the main text to the exception of a number of mismatching sites in the tropical and midlatitudes in all Southern Ocean basins (Pacific, Indian and Atlantic).

The distribution for T. sacculifer shows a clear latitudinal mismatch of the limit of presence/absence when comparing the FAME results to the unsorted MARGO dataset. It seems obvious that the latitudinal spread of G. ruber in FAME should be considered as too extended in the middle to high latitudes in both hemispheres.

The joint comparison of unsorted G. ruber and T. sacculifer distributions points to the existence of consistent zones where FAME does not predict the absence of those two species. This was noted earlier for the Benguela upwelling region. It is also visible here for the Peru–Chile upwelling and the eastern equatorial Pacific. All these zones correspond to upwelling regions (e.g., Mackas et al., 2006) and are characterized by strong contrasts in surface water properties with respect to the surrounding regions, large interannual and intra-seasonal variability, and high phytoplankton production. The existence of this consistent pattern in upwelling regions in the unsorted database confirms that G. ruber and T. sacculifer distributions are not well simulated in upwelling regions, either because nutrients are presently not accounted for in FAME or because the increased nutrient availability and/or the vertical structure of oceanic physical properties is not faithfully depicted in the 1^∘ resolution WOA13 dataset we used in input.

The unsorted distribution for G. bulloides still presents an excellent match for the limits but with some discrepancies in the equatorial and tropical latitudes; albeit MARGO unsorted data do not present a large regional consistency outside the northern coast of Brazil (where FAME also predicts the absence of G. bulloides).

Aside from some minor mismatches in the southern Indian Ocean, the conclusions for N. pachyderma are also largely unaffected by the use of all the points from the MARGO database.

To conclude, the use of all the data points regardless of the quality of the chronological control in the MARGO Late Holocene database does not add much new information, especially since the data points should be considered with caution as they could correspond to a different climate regime than the Late Holocene.