On the validity of foraminifera-based ENSO reconstructions

On the validity of foraminifera-based ENSO reconstructions Brett Metcalfe1,2, Bryan C. Lougheed1,3, Claire Waelbroeck1 and Didier M. Roche1,2 1Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, F91191 Gif-sur-Yvette, France 2Earth and Climate Cluster, Department of Earth Sciences, Faculty of Sciences, VU University Amsterdam, de Boelelaan 5 1085, 1081 HV, Amsterdam, The Netherlands 3Department of Earth Sciences, Uppsala University, Villavägen 16, 75236 Uppsala, Sweden Correspondence to: Brett Metcalfe (brett.metcalfe@lsce.ipsl.fr or b.metcalfe@vu.nl)

A chief parameter that has been employed in previous ENSO proxy work using individual foraminifera analysis (IFA) is the σ(δ 18 Oc) parameter.This parameter σ(δ 18 Oc) is used to reconstruct the downcore individual foraminifera variance in δ 18 Oc, and increased σ(δ 18 Oc) is considered to correlate to increased variation in SST and, in turn, increased ENSO incidence and/or magnitude (Leduc et al., 2009;Zhu et al., 2017a) or even increased inter-annual variance in temperature.However, foraminifera are not passive recorders of environmental conditions such as SST, in that the very the ambient environment that researchers wish to reconstruct also modifies the foraminiferal population as well (Mix, 1987).Here, we use the recently developed Foraminifera as Modelled Entities (FAME) model (Roche et al., 2017) to take into account potential modulation of δ 18 Oc and the Temperature recorded in the calcite, herein Tc, by foraminifera growth, an estimate of the proxy Mg/Ca.A number of models exist to determine the foraminiferal responses to present (Fraile et al., 2008(Fraile et al., , 2009;;Kageyama et al., 2013;Kretschmer et al., 2017;Lombard et al., 2009Lombard et al., , 2011;;Roy et al., 2015;Waterson et al., 2016;Žarić et al., 2005, 2006), past (Fraile et al., 2009;Kretschmer et al., 2016) and future (Roy et al., 2015) climate scenarios, FAME uses the associated temperature and δ 18 Oeq at each grid cell to compute a time averaged δ 18 Oc and Tc for a given species.This model utilises the temperature component of the Lombard et al. (2009) equations to simulate temperature-derived growth (Kageyama et al., 2013;Lombard et al., 2009Lombard et al., , 2011) ) and subsequently growth rate-weighted δ 18 Oc (Roche et al., 2017) and Tc for Globigerinoides ruber, Globigerinoides sacculifer and Neogloboquadrina dutertrei, forced by fifty-eight years of monthly ORAS S4 ocean reanalysis temperature and salinity data (Balmaseda et al., 2013).The original Lombard et al. (2009Lombard et al. ( , 2011) ) equations are based upon a synthesis of culture studies, pooled together irrespective of experimental design or rationale, therefore they can be considered to conceptually represent the fundamental niche of a given foraminiferal species, i.e. the range in environment that the species can survive.Using the MARGO core top δ 18 Oc database (MARGO Project Members*, 2009), Roche et al. (2018) computed the optimum depth habitat (the depth habitat that exhibits the strongest correlation when comparing FAME δ 18 Oc and MARGO δ 18 Oc) for each species in the MARGO database (MARGO Project Members*, 2009).The MARGO database does not include N. dutertrei, meaning that we concentrate our efforts mainly on G. ruber and G. sacculifer.By identifying the optimum depth habitat, Roche et al. (2018) established the realised niche, i.e. the range in environment that the species can be found, for these species for the late Holocene.Unlike some foraminiferal models, FAME does not include limiting factors such as competition, respiration or predation variables as no reliable proxy exists for such parameterisation in the geological record, therefore aspects such as interspecific competition that may limit the niche width of a species are not computed.Consequently, we allow all the species of foraminifera to grow down to ~ 400 m (depending if optimal temperature conditions are met) to capture the total theoretical niche width.As the optimised depths of Roche et al. (2018) are shallower, and upper ocean water is more prone to temperature variability, our approach likely dampens both the modelled δ 18 Oc and Tc.Therefore, the sensitivity of the model was tested by applying the same procedure but with the limitation of the depth set to 60; 100 and 200 m.

Input variables (temperature; salinity and δ 18 Osw)
The temperature and salinity of the ocean reanalysis data product (Universiteit Hamburg, DE) ORA S4 (Balmaseda et al., 2013) were extracted at one-degree resolution for the tropical Pacific (-20°S to 20°N and 120°E to -70°W), with each single grid cell comprised of data for 42 depth intervals (5 -5300 m water depth) and 696 months (January 1957-December 2015).A global 1-degree grid was generated, and each grid cell was classified as belonging to one of 27 distinct ocean regions, as defined by either societal and scientific agencies, for identifying regional δ 18 Osw -salinity relationships (LeGrande and Schmidt, 2006).Using the δ 18 Osw database of LeGrande and Schmidt (2006) a regional δ 18 Osw -salinity relationship was defined, of which the salinity is the salinity measured directly at the isotope sample collection point (included within the database).Two matrices were computed; one giving values of the slope (m) and the other of intercept (c) of the resultant linear regression equations, these were used as look-up tables to define the monthly δ 18 Osw from the monthly salinity Ocean reanalysis product ORAS S4 (Balmaseda et al., 2013), which was used for the calculation of δ 18 Oeq, i.e. the expected δ 18 O for foraminiferal calcite formed at a certain temperature (Kim and O'Neil, 1997).The δ 18 Oeq is calculated from a rearranged form of the following temperature equation: Specifically, we used the quadratic approximation (Bemis et al., 1998) of Kim and O'Neil (1997), where  0 = 16.1, a = 0.09, b = -4.64 and converted from V-SMOW to V-PDB using a constant of -0.27 ‰ (Hut, 1987;Roche et al., 2017):

Foraminifera as modelled entities (FAME)
For each latitude and longitude grid, a monthly growth rate was calculated using the ORA S4 temperature and a Michaelis-Menton kinetics to predict growth calculated from culture experiments, described in Lombard et al. (2009).Growth rate was artificially constrained to arbitrary values of the upper 60; 100; 200 and 400 m for all species to reflect the symbiotic nature of various foraminiferal species.FAME is based upon FORAMCLIM which can compute eight foraminiferal species (Kageyama et al., 2013;Lombard et al., 2009Lombard et al., , 2011;;Roche et al., 2017), however comparison with a core top database has been limited to five foraminiferal species (Roche et al., 2017).Here the output of FAME is further restricted to three species that have been the main focus of foraminifera-based studies that have been used to infer ENSO variability, namely the upper ocean dwelling Globigerinoides sacculifer and Globigerinoides ruber, as well as the thermocline dwelling Neogloboquadrina dutertrei (Ford et al., 2015;Koutavas et al., 2006;Koutavas and Joanides, 2012;Koutavas and Lynch-Stieglitz, 2003;Leduc et al., 2009;Sadekov et al., 2013).

Monthly growth rate (GR) weighted δ 18 O and TC
The FAME growth rate output was used to compute the monthly depth-weighted oxygen isotope distribution (∑ , where z is the depth) for each species, using the computed δ 18 Oeq for a given latitudinal and longitudinal grid point, where () is either 60; 100; 200 or 400 m.No correction for species specific disequilibria, such as vital effect, was applied to the data.At each grid point the total growth rate (∑ ∑ , where nt is the number of time steps) was calculated and each individual monthly growth rate (∑  () =0 ) normalised to this value.A weighted histogram was constructed from the corresponding oxygen isotope values (Figure S1), which will sum to 1 (or 100 %).The rationale for constructing a (flux/GR-)weighted histogram is related to the probability of an oxygen isotope value occurring, essentially for an unweighted distribution this is , with each month contributing 1  and therefore, having a maximum value for   of t (= 696).However, in reality the monthly contribution, if remineralisation of the settling flux and benthic seafloor processes (Lougheed et al., 2018) are for now largely ignored, is modulated by various hydrographic processes and related to the temperature sensitivity (Jonkers and Kučera, 2017;Mix, 1987).Although other processes may also impact species such as mixed layer depth and nutrients these for now can be ignored.Therefore, both the resultant δ 18 Oc and flux are sufficiently perturbed/altered by T°C that for a weighted distribution each monthly contribution is . The same procedure was executed for the Tc, where the computed δ 18 Oeq is replaced with the ORA S4 temperature profile for a given latitudinal and longitudinal grid point.

Statistical analysis
The tropical Pacific Ocean is divided into four Niño regions based on historical ship tracks, from east to west: Niño 1 and 2 (0° to -10°S, 90°W to 80°W), Niño 3 (5°N to -5°S, 150°W to 90°W), Niño 3.4 (5°N to -5°S, 170°W to 120°W) and Niño 4 (5°N to -5°S, 160°E to 150°W).One index for ENSO, the Oceanic Niño Index (ONI), based upon the Niño 3.4 region (because of the region's importance for interactions between ocean and atmosphere) is a 3-month running mean of SST anomalies in ERSST.v5(Huang et al., 2017).However, Pan-Pacific meteorological agencies differ in their definition (An andBong, 2016, 2018) of an El Niño, with each country's definition reflecting socio-economic factors, therefore, for simplicity we utilise a threshold of χ ≥ +0.5°C as a proxy for El Niño, -0.5°C ≤ χ ≥ +0.5°C for neutral climate conditions and -0.5°C ≤ χ for a La Niña in the Oceanic Niño Index.Many meteorological agencies consider that five consecutive months of χ ≥ +0.5°C must occur for the classification of an El Niño event.However, here it is considered that any single month falling within our threshold values as representative of El Niño, neutral or La Niña conditions (grey bars in Figure 1).By using this threshold, three weighted histograms for each δ 18 Oc and Tc and their resultant distributions (El Niño; Neutral; and La Niña) were computed for every month and for every latitude and longitude grid-point for the 1958-2015 period.
Monthly output representing different climate conditions were compared against each another for δ 18 Oc and Tc separately.An Epanechnikov-kernel distribution was first fitted to the monthly output, after which any two desired distributions could be compared for similarity using an Anderson-Darling test.The areas where the two distributions are found to agree, at the 5% significance level, are shown as, e.g., the grey/white area in Figure's 3 to 6. Comparison, between the weighted ) and unweighted (�δ 18 O   � = As the weighted distributions are effectively probability distributions, in order to fit a distribution, we multiplied the bin counts by 1000, effectively converting probability into a hypothesised distribution.Using the repeat matrix function (MatLab function: repmat), a matrix of δ 18 Oc was produced using each bin's mid-point (δ 18 Omid-point) there is a threefold error combined with this methodology which may account for minor variation between discrete runs of the model: first the counts values were rounded to whole integers so an exact number of cells could be added to a matrix; secondly the δ 18 Omid-point was used which gives an error associated with the bin size (±0.05‰) that is symmetrical close to the distributions measures of central tendency but asymmetrical at the sides; and finally, the associated rounding error at the bin edges within a histogram (±0.005 ‰).

FAME Output
Using a basin-wide statistical test, we examine whether the mean δ 18 Oc values of a given El Niño foraminifera population (FPEN) and a given non-El Niño ('Neutral conditions') foraminifera population (FPNEU) can be expected to be significantly different at any given specific location.Where FPEN and FPNEU exhibit significantly different mean δ 18 Oc values, ENSO events can potentially be detected by paleoceanographers and unmixed using, for example, a simple mixing algorithm with individual foraminiferal analysis (Wit et al., 2013).In cases where FPEN and FPNEU do not exhibit significantly different means, then the chosen species and/or location represent a poor choice to study ENSO dynamics.We compute growthweighted δ 18 Oc (Figure 3 and 5) and temperature (Figure 6) distributions for each 1° by 1° grid cell in the fifty-eight year simulation using FAME (Roche et al., 2018), constraining the calculation to the Tropical Pacific Ocean (between -20°S to 20°N and 120°E to -70°W).
Our model produces 696 individual monthly maps for all three species (Figure 2), a series of species-specific histograms for different climatic states for each grid-point (Figure S1) and their underlying cumulative distribution functions (Figure S2).
While two of the three species (G.ruber and G. sacculifer) have similar ecologies, they show differences in their resultant δ 18 Oc for the same ocean conditions.These results can subsequently be compared to the recorded climate conditions associated with each grid cell and simulation time step: identification of timesteps that represent El Niño (EN), Neutral (NEU), and La Niña climate conditions was done using the Oceanic Nino Index (ONI) derivative (Huang et al., 2017) (Figure 1).Comparison, for each species, FAME's predicted growth-weighted δ 18 Oc distributions associated with each climate event was done using an Anderson-Darling (AD) test.This statistical test can be used to determine whether or not two distributions can be said to come from the same population.Our results show that much of the Pacific Ocean can be considered to have statistically different population between FPEN and FPNEU.We consider that the likely cause for such a remarkable result is due to FAME computing a weighted average and, therefore, the lack of a signal found exclusively within the regions demarked in Figure 1 as El Niño regions could represent how the temperature signal is integrated via an extension of the growth rate; growing season and depth habitat of distinct foraminiferal populations.To determine how robust these results are we use the 1σ values of the observed (FAME) minus expected (MARGO), as computed by Roche et al. (2018) with the MARGO core top δ 18 Oc database, as the potential error associated with the FAME model.We have computed regions in which the difference in oxygen isotopes between the two populations (Δδ 18 Oc) compared with the ADtest is smaller than the aforementioned error (Hatching in Figures 3 and 4), i.e.where the mean difference between FPEN and FPNEU is within the error.The hatched regions in Figure 3 considerably reduce the areal extent of significant difference between FPEN and FPNEU, with the remaining regions aligning with the El Niño 3.4 region (Figure 1).No such test was performed on the N. dutertrei dataset, because of its absence from the MARGO dataset.However, N. dutertrei areas may be unsuitable for the extraction of a palaeoclimate signal for sedimentological reasons, as will be outlined later.The model-driven results were assessed with the underlying observational dataset, to check if they consistent the input data (temperature and δ 18 Oeq) underwent the same statistical test (Figure 4 and Figure S3).Instead of a variable depth, we opted for fixed depths at 5, 149 and 235 m, giving a Eulerian view (Zhu et al., 2017a) in which to observe the implications of a dynamic depth habitat.By using a fixed depth, these results show that the shallowest depths produce populations that are significantly different both in terms of their mean values and their distributions.In the upper panel of Figure 4, the canonical El Niño 3.4 region is clearly visible at 5 m depth.Whilst differences exist between the temperature (Figure 4) and the FAME AD results (Figure 3), for instance close to the Panama isthmus, there are significant similarities between the plots.These plots also show that our FAME data (Figure 3), in which we allow foraminiferal growth down deeper than the depths in Roche et al. (2018), are a conservative estimate and thus are on the low-end (Figure 4), to account for potential discrepancies with depth habitats.In the original paper on depth habitats based upon temperatures derived from δ 18 Oc, Emiliani (1954) cautioned that the depth habitats obtained would represent a weighted average of the total population, and while foraminiferal depth habitats are likely to vary spatiotemporally, the average depth habitat is skewed toward the dominant signal (Mix, 1987).
To further test the model-driven results and to assess if they are still consistent when the depth limitation (400 m) is shoaled (to 60, 100, 200 m), the analysis was rerun with the aforementioned range of depths (Figure 5 and 6).Whilst it is possible to discern differences between the depths, it is important to note that a large percentage of the tropical Pacific remains accessible to palaeoclimate studies.A shallower depth limitation in the model increases the area for the 'warm' species, suggesting that the influence of a reduced variability in temperature or δ 18 Oeq with a deeper depth limit causes the differences between FPEN and FPNEU to be reduced.

Water depth and SAR
The resolution of the ocean reanalysis data for the time period 1958-2015 would essentially be analogous with a sediment core representing 50 yr -1 cm -1 (or 20 cm -1 kyr -1 ).Based on our analysis, such a hypothetical core with a rapid sediment accumulation rate (SAR) could allow for the possible disentanglement of El Niño related signals from the climatic signal using IFA, but only in a best-case scenario involving minimal bioturbation, which is unlikely in the case of oxygenated waters.Indeed, one should view discrete sediment intervals, and the foraminifera contained within them, as representative of an integrated multi-decadal or even multi-centennial signal (Lougheed et al., 2018;Peng et al., 1979), as opposed to other proxies such as corals (Cole and Tudhope, 2017), speleotherms (Chen et al., 2016), or varves in which distinct layers correspond to discrete time intervals (true 'time-series' proxies).Therefore, in order to reliably extract short-term environmental information from foraminiferal-based proxies, the signal that one is testing or aiming to recover must exhibit a large enough amplitude in order to perturb the population by a significant degree from the background signal, otherwise it will be lost due to the smoothing effect of bioturbation (Hutson, 1980a;Löwemark, 2007;Löwemark et al., 2005Löwemark et al., , 2008;;Löwemark and Grootes, 2004) upon the downcore, discrete-depth signal (Cole and Tudhope, 2017;Mix, 1987).Similarly, the individual characters of El Niño events, which are very short in duration, become lost in the bioturbated sediment record.New geochronological tools, such as dual 14 C-δ 18 O measurements on single foraminifera (Lougheed et al., 2018), show that low sedimentation rate cores can have large variances in age between individual foraminifera present within a discrete 1 cm depth interval (Berger and Heath, 1968;Lougheed et al., 2018).In the case of a high SAR core (20 cm -1 kyr -1 ), assuming a sediment mixing layer of 10 cm, benthic seafloor processes would result in a minimum 700-year 1σ confidence interval for a 1 cm discrete depth age (Lougheed et al., 2018).A consequence of bioturbation in sediment cores is that a series of high magnitude, but low frequency El Niño events could be smoothed out of the downcore, discrete-depth record.To investigate which areas of the sea floor can potentially preserve a foraminiferal downcore signal, we overlay our results with a map of time-averaged deep-sea SAR (Olson et al., 2016), adapted by Lougheed et al. (2018) to also show seafloor areas under the CCD depth, where carbonate material is not preserved (Berger, 1967(Berger, , 1970b)).Of the total area where FPEN is significantly different from FPNEU (i.e.those areas where planktonic foraminiferal flux is suitable for reconstructing past ENSO dynamics), only a small proportion corresponds to areas where the sea floor is both above the CCD (< 3500 mbsl) and SAR is at least 5 cm/ka (Figure 3D).Also of note is that fact that much of the canonical El Niño 3.4 region (Wang et al., 2017) used in oceanography (Figure 1) is also excluded from these suitable areas in Figure 3D.

Palaeoceanographic Implications
Palaeoceanographic (Ford et al., 2015;Garidel-Thoron et al., 2007;Koutavas et al., 2006;Koutavas and Joanides, 2012;Koutavas and Lynch-Stieglitz, 2003;Leduc et al., 2009;Sadekov et al., 2013;White et al., 2018) and palaeoclimatological archives (Moy et al., 2002) have been used to indirectly and directly study past ENSO.Several authors have focussed on individual foraminifera analysis (IFA) or pooled foraminiferal analysis in the Pacific region, either for trace metal or stable isotope geochemistry.Given the complexity in reconstructions of trace metal geochemistry (Elderfield and Ganssen, 2000;Nürnberg et al., 1996), and the potential error associated with determining which carbonate phase is first used when foraminifera biomineralise (Jacob et al., 2017), the focus here has been on the δ 18 Oc.The resultant data of such studies have been used to infer a relatively weaker Walker circulation, a displaced ITCZ and equatorial cooling (Koutavas and Lynch-Stieglitz, 2003); both a reduction (Koutavas and Lynch-Stieglitz, 2003) and intensification (Dubois et al., 2009) in eastern equatorial Pacific upwelling; and both weakened (Leduc et al., 2009) and strengthened ENSO variability (Koutavas and Joanides, 2012;Sadekov et al., 2013) during the LGM.A number of these results are contentious, for instance the reduction in upwelling in this region (Koutavas and Lynch-Stieglitz, 2003) is contradicted by Dubois et al. (2009), who used alkenones (i.e.,  37 ′ ratios) to suggest an upwelling intensification.Whilst the  37 ′ proxy has problems within coastal upwelling sites (Kienast et al., 2012) it does not discount their claim, especially considering that δ 18 O records can themselves be influenced by salinity upon the δ 18 Osw component (Rincón-Martínez et al., 2011) and the potential influence of [CO3 2-] upon foraminiferal δ 18 Oc ( de Nooijer et al., 2009;Spero et al., 1997;Spero and Lea, 1996).Furthermore, the discrepancies between marine cores' is worth noting, especially considering that terrestrial records suggest the number of El Niño events per century in the early Holocene (8-6 ka BP) was minimal (Moy et al., 2002), with between 0 and 10 events occurring per century.This dampened ENSO is observed within lake core colour intensity and records driven primarily by precipitation (Cole and Tudhope, 2017;White et al., 2018).If the number and magnitude of ENSO events were reduced, the relatively low downcore resolution of marine records may not accurately capture the dynamics of such lower amplitude ENSO events using existing methods.The possibility of a marine sediment archive being able to reconstruct ENSO dynamics comes down to several fundamentals: the time-period captured by the sediment intervals (a combination of SAR and bioturbation), the frequency and intensity of ENSO events, as well as the foraminiferal abundance during ENSO and non-ENSO conditions.
The results presented here imply that much of the Pacific Ocean is not suitable for reconstructing ENSO studies using palaeoceanography, yet several studies have exposed shifts within σ(δ 18 Oc) of surface and thermocline dwelling foraminifera.One can, therefore, question what is being reconstructed in such studies.One artefact of sampling is the potential occurrence of aliasing (Pisias and Mix, 1988;Wunsch, 2000;Wunsch and Gunn, 2003): a fundamental problem with proxy records is that they can be confounded by local regional climate, and/or ENSO's teleconnections, that mimic ENSO changes albeit at a different temporal frequency (Dolman and Laepple, 2018).Our own analysis using our FAME δ 18 Oc and Tc output mimics foraminiferal sedimentary archives, pooling several decades worth of data in which the resolution is coarse enough to obscure and prevent individual El Niño events being visible but allowing for some kind of long-term mean state of ENSO activity to be reconstructed (Cole and Tudhope, 2017).The results of our Anderson-Darling test may be unduly influenced by the Pacific decadal variability (PDV), also referred to as the Pacific Decadal Oscillation (PDO) (Pena et al., 2008).In much of the tropical Pacific the ratio of decadal to interannual σSST suggests that they are comparable in magnitude, therefore fluctuations in SST are more obviously apparent outside of the purely canonical regions of ENSO (Wang et al., 2017).It could be that the areas outside of these canonical ENSO regions (Figure 1) reflect the PDO (Pena et al., 2008;Wang et al., 2017).The use of σ(δ 18 Oc) as an indicator of ENSO is dependent on whether the variance is dominated by interannual variance (Ford et al., 2015;Thirumalai et al., 2013;Zhu et al., 2017a).Zhu et al. (2017) computed the total variance change with and without the annual cycle suggesting that, for some cores the increased assumed ENSO variability at the LGM as deduced by proxy records (Koutavas et al., 2006;Koutavas and Joanides, 2012;Koutavas and Lynch-Stieglitz, 2003) may be purely a by-product of the annual cycle.

Limitations of the methods applied and assessment of model uncertainties
For simplicity we have assumed that our model is 'perfect', of course that is inaccurate, there are four potential sources of error: the input variables (temperature, salinity and their conversion into δ 18 Osw and δ 18 Oeq); the model's error with respect to real world values (Roche et al., 2018); the statistical test's errors (associated Type I -in which attribution of significance is given to an insignificant random event, a false 'positive' -and Type II -in which a significant event is attributed to be insignificant, a false 'negative' -errors); and reducing the complexities of foraminiferal biology via parameterization.The input variables can have errors associated with both the absolute values of temperature and salinity used here; and the  (1997) equation and the foraminiferal microenvironment (de Nooijer et al., 2008(de Nooijer et al., , 2009;;Spero et al., 1997;Spero and DeNiro, 1987;Spero and Lea, 1996) which has further implications due to the upwelling of cool, low pH, waters in the eastern Tropical Pacific (Cole and Tudhope, 2017;Raven et al., 2005).The spatial variability in salinity, particularly within regions underlying the intertropical convergence zone (ITCZ) and the moisture transport from the Caribbean into the eastern Pacific along the topographic low that represents Panama Isthmus, the resultant conversion of salinity to δ 18 Osw and then δ 18 Oeq may contain further error.If such errors are independent of the absolute value of the variable, i.e. the error on cold temperature is the same and not larger than warm temperatures, then the error terms effectively cancel one another out.A point of note, is that the δ 18 O to °C conversion of Kim and O'Neil (1997) is considered to be marginally larger at the cold end then at the warm end (0.2 ‰ per 1°C to 0.22 ‰ per 1°C) than that originally discerned (O'Neil et al., 1969).The comparison of the temperature signal produced here (Tc) to a value corresponding to that reconstructed from measurements of Mg/Ca should be done with caution.As stated previously, several other parameters can alter this technique, this includes abiotic effects such as salinity (Allen et al., 2016;Gray et al., 2018;Groeneveld et al., 2008;Kısakürek et al., 2008) or carbonate ion concentration (Allen et al., 2016;Evans et al., 2018;Zeebe and Sanyal, 2002); biotic effects such as diurnal calcification (Eggins et al., 2003;Hori et al., 2018;Sadekov et al., 2008Sadekov et al., , 2009;;Vetter et al., 2013); or additional factors such as sediment (Fallet et al., 2009;Feldmeijer et al., 2013) or specimen (Barker et al., 2003;Greaves et al., 2005) 'cleaning' techniques.Given the role of Mg in inhibiting calcium carbonate formation, the manipulation of seawater similar to the modification of the cell's pH (de Nooijer et al., 2008(de Nooijer et al., , 2009) ) may aid calcification and explain the formation of low-Mg by certain foraminifera (Zeebe and Sanyal, 2002), however scaling these processes up to a basin-wide model is beyond the remit of this current paper.
Our modelling results also depend upon the species symbiotic nature and potential genotypes.For instance, mixotrophs, those organisms that utilise a mixture of sources for energy and carbon (planktonic foraminifera such as G. ruber; and/or G. sacculifer) can outcompete heterotropic (or photoheterotropic) organisms (planktonic foraminifera such as N. pachyderma; N. incompta) especially in stratified-oligotrophic waters.FORAMCLIM has placed both G. bulloides and N. dutertrei as symbiont barren for the purposes of ecophysiological parameters used for computing expected potential habitats for both future and past scenarios.However, recent culture studies have shown that G. bulloides Type IId instead of being symbiont barren actually plays host to endobiotic cyanobacteria (Bird et al., 2017), whilst N. dutertrei Type Ic contains what are reputed to be algal symbionts (Bird et al., 2018).They do however, confirm that N. incompta is symbiont barren (Bird et al., 2018).The authors have suggested that ecophysiological models, such as FORAMCLIM, should be updated to include both the photosynthetic changes and resultant respiration parameters (Bird et al., 2018).Whilst FAME uses only the temperature dutertrei) are non-symbiotic.A species that hosts symbionts will likely have a restricted temperature that is associated with the temperature tolerance of their symbionts, given that the next geneartion of a species of planktonic foraminifera must be re-infected with their symbionts.Variation in cryptic species of G. bulloides do exist, given that genotypes of morphotypes exhibit distinct environmental preferences (Darling et al., 2000(Darling et al., , 1999;;Huber et al., 1997;Morard et al., 2013;de Vargas et al., 1999de Vargas et al., , 2002)), especially between genotypes Type Ia, Ib, IIa, IIb, and IIc and the genotype Type IId that has a restricted geographic distribution (Darling et al., 2004).Morard et al. (2013) simulated the impact of genotypes upon palaeoceanographic reconstructions (in particular transfer functions) using a theoretical abundance, calculated with a best-fit gaussian response model, depending upon SST later using a similar approach (Morard et al., 2016) to deduce the impact upon δ 18 O.Incorporation of both a theoretical genotype abundance (Morard et al., 2013) and ecophysiological tolerances of different genotypes (Bird et al., 2018) within an ecophysiological model could further reduce error within modelling of planktonic foraminiferal habitats, and thus reduce data-model comparison error.

Conclusion
Previous work has focussed on comparing either the difference between time slices or between downcore and modern IFA variance, using a reduction in σ(δ 18 Oc) as an indicator of reduced ENSO (Zhu et al., 2017a).Concentrating on the period spanning the instrumental record by using the FAME module, we present a species-specific δ 18 O and Temperature output (Tc) for ocean reanalysis data to validate the use of foraminiferal proxies in reconstructing ENSO dynamics in the Pacific Ocean.Overall, our results suggest that foraminiferal δ 18 O for a large part of the Pacific Ocean can be used to reconstruct ENSO, especially if an individual foraminiferal analysis (Lougheed et al., 2018;Wit et al., 2013) approach is used (Ford et al., 2015;Koutavas et al., 2006;Koutavas and Joanides, 2012;Koutavas and Lynch-Stieglitz, 2003;Sadekov et al., 2013;White et al., 2018), contrary to previous analysis (Thirumalai et al., 2013).However, the sedimentation rate of ocean sediments in the region is notoriously slow (Olson et al., 2016) and much of the ocean floor is under the CCD.These factors reduce the size of the area available for reconstructions considerably (Lougheed et al., 2018), thus precluding the extraction of a temporally valid palaeoclimate signal using long-standing methods.We further highlight that the conclusions drawn from foraminiferal reconstructions should consider both the frequency and magnitude of El Niño events during the corresponding sediment time interval (with full error) to fully understand whether or not a strengthening or dampening occurred.The use of ecophysiological models (Kageyama et al., 2013;Lombard et al., 2009Lombard et al., , 2011) ) are not limited to foraminifera and provide an important way to test whether proxies used for palaeoclimate reconstructions are suitable for the given research question.
made for El Niño climate; combined Neutral and La Niña climates; only Neutral climate; and only La Niña climate.

Clim.
Past Discuss., https://doi.org/10.5194/cp-2019-9Manuscript under review for journal Clim.Past Discussion started: 1 March 2019 c Author(s) 2019.CC BY 4.0 License.limitation of input values to a single value per month.Whilst it is possible to interpolate to a daily resolution, this is problematic for two reasons: (1) daily temperature records have much more high frequency oscillations than the data here and (2) the lifecycle of a single foraminifera is approximately monthly, therefore by using monthly data it provides an estimate of the average population signal.Conversion of salinity and temperature into δ 18 Osw and δ 18 Oeq uses a quadratic approximation, one source of error is the unknown influence of carbonate ion concentration on both the Kim and O'Neil

Clim.
Past Discuss., https://doi.org/10.5194/cp-2019-9Manuscript under review for journal Clim.Past Discussion started: 1 March 2019 c Author(s) 2019.CC BY 4.0 License.component of FORAMCLIM (Roche et al., 2018) and we have only modelled N. dutertrei, it is important to note that there are distinctions between the fundamental niche that FAME computes, i.e. the conditions that an organism can survive, and the realised niche, i.e. what an organism actually occupies given limiting factors within the environment.Likewise, FAME and FORAMCLIM are based upon the original culture experiments that assumed that both species (G.bulloides and N.

Figure 1 .
Figure 1.Oceanic Niño Index and the temperature anomaly for a single El Niño event.(Top) Oceanic Niño Index (ONI), sourced from a 3-month running mean of SST anomalies in ERSST.v5 of the Niño 3.4 region (Huang et al., 2017).Grey vertical bars represent the periods in which El Niño-like conditions exist using a simple one-month threshold.(Bottom) The sea surface temperature difference between week beginning 1st December 1997 minus the long-term climatic mean (1971 -2000) for December.The 1997 -1998 El Niño represents an EP-ENSO.The long term monthly climatology, the NOAA optimum interpolation (OI) SST V2, based upon the methodology of Reynolds and Smith (Reynolds and Smith, 1995) using two distinct climatologies for 1971 -2000 and 1982 -2000 (Reynolds et al., 2002).Boxes represent the Nino region: Niño 1 and 2 (0° --10°S, 10

Figure 2 .
Figure 2. A snapshot of the output of FAME.Each panel represents an individual species (Top Panel G. sacculifer; Middle Panel G. ruber and Bottom Panel N. dutertrei) δ 18 O for a single time step (t = 696).The δ 18 O for each grid-point is based upon FAME, a computed integrated signal based upon the growth rate and the δ 18 Oeq values.Values in per mil (‰ V-PDB).

Figure 3 .
Figure 3. Maps of those locations that can reconstruct El Niño events for G. ruber; G. sacculifer; N. dutertrei and all species.Upper (three) panels represent the Anderson-Darling test results for the species G. ruber; G. sacculifer; and N. dutertrei where the black areas reflect latitudinal and longitudinal grid points that failed to reject the null hypothesis (H0) and therefore the foraminiferal population (FP) of the El Niño is similar to the Non-El Niño, and therefore the distribution between the neutral climate and El Nino cannot be said to be different (FPEl Niño = FPNon-El Niño).White/grey areas reflect regions in which the H1 hypothesis is accepted, therefore the distributions can be said to be unique (FPEl Niño ≠ FPNon-El Niño).Hatching represents the FAME species specific values of 0.32 ‰ and 0.20 ‰.The final, bottom, panel represents a synthesis of the upper panels areas in which reconstructions are possible, represented by both the blue and white area.The black area represents those areas where the null hypothesis for all three species cannot be rejected.Blue represents sedimentological limitations that preclude the extraction of a 10

Figure 4 .
Figure 4. Results of an Anderson-Darling test between El Nino and Neutral climate conditions based upon the Temperature input data: Fixed depth.White represents where the H1 hypothesis is accepted, therefore the distributions of the foraminiferal population (FP) for El Niño and Non-El Niño can be said to be unique (FPEl Niño ≠ FPNon-El Niño).Black that the null hypothesis (H0), 5

Figure 5 .
Figure 5. Oxygen isotope (δ 18 Oc): Maps of those locations that can reconstruct El Niño events for (Left Panel) G. ruber; (Middle panel) G. sacculifer; and (Right panel) N. dutertrei.Panels represent the Anderson-Darling test results for the species G. ruber; G. sacculifer; and N. dutertrei where the black areas reflect latitudinal and longitudinal grid points that failed to reject the null hypothesis (H0) and therefore the foraminiferal population (FP) of the El Niño is similar to the Non-El Niño, and therefore the distribution between the neutral climate and El Nino cannot be said to be different (FPEl Niño = FPNon-El Niño).White/grey areas reflect regions in which the H1 hypothesis is accepted, therefore the distributions can be said to be unique (FPEl Niño ≠ FPNon-El Niño).Hatching represents the FAME species specific values of 0.32 ‰ and 0.20 ‰.Each row represents a particular depth limit for the model, from top to bottom: 60 m; 100 m; 200 m and 400 m.10

Figure 6 .
Figure 6.Temperature signal (TC): Maps of those locations that can reconstruct El Niño events for (Left Panel) G. ruber; (Middle 5