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 data – model integration studies. To overcome this shortcoming, we developed the FAME (Foraminifers As Modeled Entities) module. The module predicts the presence or absence of commonly used planktonic foraminifers, and their oxygen-18 5 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 MARGO Late Holocene planktonic foraminifer calcite oxygen-18 and abundances data sets. 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.


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 ra-tios of oxygen-18 to oxygen-16 in calcite (hereafter R 18 O c ).The recognition that different species of foraminifers from the same sediment core yielded different R 18 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 18 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 18 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 18 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 18 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 18 O c record: seawater temperature and oxygen-18 abundance.So far, water-isotopebased model-data studies have generally compared planktonic oxygen-18 ratios to equilibrium calcite R 18 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 18 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 culturebased 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 18 O c in a waterisotope-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.

Methodology
The calcite oxygen-18 / oxygen-16 ratio (reported in δ 18 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 δ 18 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 δ 18 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 δ 18 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 δ 18 O c is then a weighted sum of these conditions in space and time.

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 where µ(T , k) is the growth rate at temperature T for the species k, µ(T 1 , k) is the growth rate for a chosen reference temperature T 1 (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.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 fourdimensional space as µ (T , k) = µ (T (x, y, z, t) , k) . (2) To avoid numerical issues in the code, we limit the value of µ(T , k) on the low end as follows: otherwise. (3) Given a four-dimensional input field for oceanic temperatures and δ 18 O of seawater, the equilibrium inorganic calcite δ 18 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): it becomes 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 δ 18 O equilibrium value, averaged over time and the surface layer (typically the first 50 m) to compare model results and measured δ 18 O c from planktonic foraminifers.In the following, we will refer to this method as the "old method", written formally as 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 δ 18 O c both in time and in the horizontal and vertical space according to the population abundances using the foraminifer growth for-mula (1).We thus write 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 δ 18 O provided as a four-dimensional field and a given species k, we compute this species' δ 18 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 δ 18 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 δ 18 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.

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.

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 δ 18 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 (δ 18 O sw ) raw data in numerous locations, we derived a δ 18 O sw dataset based on seawater salinity to δ 18 O sw relationships.This dataset is built in two steps: (a) derivation of regional δ 18 O sw -salinity relationships from GISS δ 18 O sw and salinity (Schmidt et al., 1999) clustered by oceanic regions, and (b) computation of a δ 18 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 δ 18 O sw dataset are given in Appendix A.
As an independent test of the FAME results, we use the planktonic δ 18 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.

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 δ 18 O c simulated by FAME driven by the World Ocean Atlas 2013 temperature and derived seawater δ 18 O datasets is closest on average to MARGO Late Holocene δ 18 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;Geosci. Model Dev., 11, 3587-3603, 2018 www.geosci-model-dev.net/11/3587/2018/(Waelbroeck et al., 2005).We computed a confidence interval σ  (1996); Bauch et al. (1997Bauch et al. ( , 2002)); Simstich et al. (2003); Kuroyanagi and Kawahata (2004) 1 An encrustation term of +0.1 permil is taken into account in the case of N. pachyderma (see text).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 δ 18 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 δ 18 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 δ 18 O values than G. ruber.These heavier δ 18 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).

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.

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 lowerbound 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 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.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 δ 18 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 δ 18 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).

Geographical distribution
To further check our methodology against the MARGO Late Holocene dataset, we compare the zone of presence Ocean regions where FAME predicts that the species is present at some time of the year (µ > 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. 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.
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 cov-erage 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.
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 com- parison 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 tempera- ture.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).

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 δ 18 O c .To highlight the effect of seasonal and vertical weighting of the δ 18 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 δ 18 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 δ 18 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 18 O c in what follows.
Applying a spatially homogeneous temperature change should result is a quasi-homogeneous temperature change in the equilibrium calcite 18 O c , following Eq.( 6).It is indeed what is obtained in Fig. 6e, with 18 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 18 O c with values down to −0.75 and up to 1 permil.All species share a common pattern of lower 18 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 18 O c for G. bulloides is a consequence of its growth over a large range of temperatures (Eq.1).In general, the maximum simulated 18 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.

Summary and conclusions
We developed the FAME (Foraminifers As Modeled Entities) module to account for planktonic foraminifer speciesspecific 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 δ 18 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 δ 18 O c are much more spatially variable than the computed change in equilibrium surface calcite.Coupling the FAME module to isotopeenabled 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  1) over the full temperature range considered, with added lower-range and upperrange 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).

A1 Construction of the basin masks
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 δ 18 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  A1.
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).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.

A2 Computation of the δ 18 O sw -salinity relationships
The basins defined in the previous section are then used to cluster the raw data, δ 18 O sw and salinity, of the GISS database (Schmidt et al., 1999)   The "*" symbol highlights the regions where FAME and WOA regions do not cover the same area (see text).
ity 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 δ 18 O sw in the WOA geographical grid from the WOA salinity fields.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.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.

Figure 2 .
Figure2.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.

Figure 3 .
Figure3.Model-data comparison of species' abundances.Ocean regions where FAME predicts that the species is present at some time of the year (µ > 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.

Figure 4 .
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).

Figure 5 .
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.

Figure 6 .
Figure 6.δ 18 O c response to a horizontally and vertically homogeneous 4 • C cooling applied to the WOA13 dataset, 18 O c .Results are expressed in permil for each species (a-d) and for the equilibrium surface calcite approach.

Figure A1 .
Figure A1.Growth functions corresponding to Eq. (1) over the full temperature range considered, with added lower-range and upperrange curves (respectively, short and long dashed curves) chosen to mimic the 95 % confidence interval ofLombard et al. (2009).Data points are the original data obtained byLombard et al. (2009).
Appendix A: Derivation of the a reference δ 18 O sw dataset We constructed our reference δ 18 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 δ 18 O sw -salinity relationships for each of these basins and (c) use of these relationships to obtain a δ 18 O sw field at WOA spatial resolution.

Figure A2 .
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 TableA1.

Figure B1 .
FigureB1.Model-data comparison of species' abundances.Ocean regions where FAME predicts that the species is present at some time of the year (µ > 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.

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 in the respective basins.Furthermore, only data locations where both δ 18 O sw and salin- Geosci.Model Dev., 11, 3587-3603, 2018www.geosci-model-dev.net/11/3587/2018/

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.