Modulation of late Pleistocene ENSO strength by the tropical Pacific thermocline

The El Niño Southern Oscillation (ENSO) is highly dependent on coupled atmosphere-ocean interactions and feedbacks, suggesting a tight relationship between ENSO strength and background climate conditions. However, the extent to which background climate state determines ENSO behavior remains in question. Here we present reconstructions of total variability and El Niño amplitude from individual foraminifera distributions at discrete time intervals over the past ~285,000 years across varying atmospheric CO2 levels, global ice volume and sea level, and orbital insolation forcing. Our results show a strong correlation between eastern tropical Pacific Ocean mixed-layer thickness and both El Niño amplitude and central Pacific variability. This ENSO-thermocline relationship implicates upwelling feedbacks as the major factor controlling ENSO strength on millennial time scales. The primacy of the upwelling feedback in shaping ENSO behavior across many different background states suggests accurate quantification and modeling of this feedback is essential for predicting ENSO’s behavior under future climate conditions.


Individual foraminifera selection.
To simulate the selection of individual foraminifera, we randomly selected 80 monthly temperatures from synthetic time series based on the SODA 2.1.6 mixed-layer temperature reanalysis data set. We altered the time series for both length (by appending copies of the time series) and for changes in ENSO amplitude and/or frequency to determine the parameters that ENSO (and more specifically, El Niño) change is most apparent. For each change in ENSO parameters, we performed the random selection of monthly temperature data 1000 times and applied random temperature uncertainty from a normal distribution with standard deviation equal to our analytical uncertainty (±0.47 °C) to simulate selection and analysis of individual foraminifera. We then calculated quantiles for each realization, then the mean value for every quantile from all realizations. We calculated a 90% confidence interval about that mean, and plotted the mean quantiles vs. the quantiles of the original, unmodified temperature series and the calculated confidence interval to test whether ENSO change is identifiable. We find that simulated sampling of individual foraminifera from an altered SST distribution is capable of capturing change in ENSO amplitude relative to the reference interval, and such changes appear in the tails of the distribution (Supplementary Figure 1). Changes in ENSO frequency does alter the distributions in the same fashion, with the same response, or with the same consistency, demonstrating that change in the tails of the distribution is likely the result of amplitude change (Supplementary Figure 2). We find that the number of individuals (60 vs 80) and the length of the synthetic time series the foraminifera are chosen from (50 vs 500 years) has a minor effect on the ability of the method to identify changes in the synthetic time series, and that denser sampling leads to higher chances of identification at lower ENSO change. Sampling rates comparable to our sample rates (80 specimens over 500 years) demonstrate that sampling analogous to individual foraminifera analysis will indicate decreases in ENSO amplitude of 19% and increases in ENSO amplitude of 21% in over 90% of cases using the 98 th quantile as a metric. Using the 96 th quantile as a metric for ENSO change, ENSO amplitude changes of 21% and 25% for decrease and increase, respectively, are captured in 90% of simulations. For the 94 th quantile, change in ENSO amplitude of 27% is captured in 90% of simulations.

Diagnostic features of quantiles
In these synthetic SST distributions, we find that the 94 th , 96 th and 98 th quantiles respond strongly and linearly to changes in ENSO amplitude (Supplementary Figure 3). Significant changes in these quantiles may be identified with ENSO changes of ±10% in our SST simulations. In our simulations, the 98 th quantile is the most sensitive to changes in ENSO amplitude, while changes in ENSO frequency or the seasonal cycle are not as sensitive at the 98 th quantile. However, this quantile is subject to the highest frequency of false positive results, and thus analysis of additional quantiles and the total number of quantiles is necessary. Intervals with multiple quantiles showing deviation from the reference population have increasingly lowered frequencies of false positives, and we find that the mean of the 94 th -98 th quantile has a lower frequency of false positives than the 98 th quantile while retaining diagnostic capability (Supplementary Figure 3).

Centennial-scale variability
We tested the impact of centennial-scale variability on the distribution of mixed-layer temperatures by generating multi-century time series beginning with 50-year modern monthly reanalysis temperature data and repeatedly appending the modern reanalysis data set. 2000 (Supplementary Figure 4) and 400 (Supplementary Figure 5) year time series were generated to bracket the length of time a given sediment interval in our cores may represent. In all simulations, random temperature anomalies from a normal distribution of mean of zero and standard deviation equal to initial reanalysis data set standard deviation were added. Centennial scale variability was simulated by changing the properties of the resulting temperature data. We tested warming and cooling trends, sudden changes in average temperature, and random changes temperature in the appended data. To test trends, a trend of either + or -0.05°C was added to each appended temperature set, resulting in an overall warming or cooling trend of ±0.1°C/ century. Changes in average mixed layer temperature were simulated by altering the mean temperature of one half of the appended intervals by ±0.5°C. The final test was the random addition of data with either a mean of +0.5°C, -0.5°C, or unchanged data. 80 months of the resulting temperature series were then randomly selected to simulate selection of individual foraminifera. Quantiles 0-100 at 2% intervals were calculated for these distributions. This process was repeated 1000 times and the median values for each quantile was calculated. These quantiles were then used for Q-Q analysis against the modern mixed-layer reanalysis data. We found that no simulations showed significant change in the distributions. The warming/cooling trends showed deviation, although not significant, from the 1:1 line. In the case of the 2000-year trends, an increase in mixed-layer temperature variability was observed as an increase in median absolute deviation (MAD). However, as this is unaccompanied by significant changes in the distributions, it is unlikely that our results represent this sort of trend. We conclude from these tests that non-ENSO centennial-scale variability is unlikely to manifest itself as a change in the extreme ends of the distributions that are used as indicators of ENSO strength.

Supplementary Note 2: Age Model
The age model for core 14MC1 is based on previous published radiocarbon dates ( 2  boundaries, as identified visually in the records (Supplementary Figure 6). We find that the records show no appreciable or consistent offset between benthic and planktic records. Further, creation of new records using Analyseries at selected MIS boundary transition points resulted in planktic and benthic depth-age models that were nearly identical.
Existing age models for core 17PC have been based on an orbitally-tuned Globigerinoides ruber δ 18 O stratigraphy, core-top radiocarbon dating, and Monte-Carlo based cross-correlation on sediments <150 ka ( 7 , 13 ). For this study, we realigned the oxygen isotope stratigraphy to the  Sample ages for core 17PC were linearly interpolated between averaged age control points.
The Analyseries generated age model was compared to the additional age models, including those based on the LR04 benthic stack alignment (JLS)( 7 ), the Monte-Carlo cross correlation methods ( 13 ), and the Hidden Markov Method Probstack generated δ 18 O stratigraphy ( 14 ) (Supplementary Figure 7). Our ages, generated using Analyseries, are within 5ky of LR04 and Probstack predicted ages with the exception of MIS 8 where there is significant disagreement between the JLS and Probstack models. Our model falls between these and we believe more accurately aligns with state transitions (e.g., glacial-interglacial) during this MIS 8 interval, and with state changes in MIS 5 and 7. chamber showed that variability is similar between the two methods, but absolute Mg/Ca values are consistently offset by +0.2 mmol/mol (outside-in minus inside-out). This offset is consistent across measured Mg/Ca value (~2-5 mmol/mol), and may be the result of surface effects from the concave curved inner surface of the foraminifera. Our multiple-chamber method, however, allows us to calculate Mg/Ca ratios that include earlier ontogenetic signals from the f1 chamber and shows close correspondence with Mg/Ca ratios derived from traditional cleaning and dissolution methods.

Analytical and calibration uncertainty
Average intra-chamber standard deviation (1σ) of the three analyses of the f0 chamber is ±0.112 mmol/mol, and of the two analyses of the f1 chamber is ±0.107 mmol/mol. This translates to ±0.24 °C for f0 and ± 0.24 °C for f1 in a foraminifer with an average Mg/Ca ratio of 5 mmol/mol. Anomalies in 27 Al, 55 Mn, and 66 Zn were identified to indicate possible contamination from clays, metal hydroxides, or other sources. As such contamination was uncommon in these sediments, suspect chambers/individuals were removed from analysis. Trace element peaks associated with the inside and outside layers of foraminifera shells were automatically removed from the analysis via an automated function. NIST glass standard 610 was repeatedly analyzed during each analytical run to compute absolute response to, and drift in, elemental intensities during and across sample runs. NIST glass was analyzed at 4Hz and fluence between 0.59 and 0.76 J/cm 2 sufficient to obtain adequate signal strength and ablation duration for analysis. The

Impact of analytical uncertainty on ENSO change
To assess the impact of this analytical uncertainty of our determination of whether an interval shows changed El Niño amplitude, we applied bootstrap resampling to randomly altered temperature distributions. We generated random temperature uncertainties from a normal distribution with mean zero and standard deviation of 0.47°C, and applied these random anomalies to each individual foraminifera temperature. We repeated this process 10 4 times and performed Q-Q analysis on each realization for each sample interval (Supplementary Table 3).
For each realization, we determined whether each sample interval displayed changed El Niño amplitude, the number of significantly changed quantiles, and the effect on median absolute deviation. Intervals with increased El Niño amplitude in multiple quantiles are unlikely to change, intervals of only slightly enhanced El Niño amplitude (e.g., 1 quantile) may be influenced by uncertainties. Intervals with decreased El Niño amplitude are more sensitive to analytical uncertainty. While Interval MIS 7e is likely to display reduced El Niño amplitude with analytical uncertainty applied, but MIS 7a is highly sensitive to these results, and, as noted, the finding of reduced El Niño amplitude compared to the modern mixed layer is not robust. In most cases, analytical uncertainty increased the total variability of the population, although the median MAD of the Monte Carlo populations is within 1 standard deviation of the actual sample MAD.

Age and analytical uncertainty and correlations
We incorporate age and analytical uncertainty into our analysis of the relationship between MAD, the mean 94 th -98 th quantile temperature anomalies of our sample intervals and climate boundary conditions. We used weighted bivariate linear regression (WLR) routines ( 17 ) to determine correlations and statistical significance accounting for these uncertainties, and incorporated Monte Carlo techniques in order to estimate our uncertainties. First, we interpolate the value of each climate parameter at each of our interval ages. We estimate the uncertainty about this value by generating age uncertainties for each record, taking in to account the uncertainty of the specific climate record and the uncertainty of our age model in quadrature. We perform Monte Carlo simulations to generate multiple age models for each record by random alteration of the ages by an amount drawn from a normal distribution with the combined age uncertainty. We then randomly applied analytical uncertainty for each record from a normal population with standard deviation of the analytical uncertainty, and then re-interpolated the data at the resulting age using the resulting analytical value. Total uncertainty about that value for use in the WLR algorithm is then calculated as one half of the ~2-sigma range. We estimate the 2sigma range from our empirical Monte Carlo data using the 2.5 and 97.5 quantiles, which utilizes most of the generated data. Ages uncertainties for each record were taken from the original references when possible. Age and analytical uncertainty for insolation parameters was considered zero, as these parameters are mathematically solutions ( 18 ). Age uncertainty for the E-W SST gradient was estimated at 2ky, which is two times the time steps of the record. Analytical uncertainty in the original record was estimated at 0.5 °C. Age uncertainty for the latest portion of the ITCZ Ti MAR record was determined from published 14C dates ( 19 ). Ages for the remainder of the record used an alignment method similar to ours using Analyseries. We conservatively assigned this record a value of 2x our age errors. No analytical uncertainty was reported for this record, but we applied a 0.5% uncertainty factor based on the mean value of the Ti MAR data ( 20 ). The mixed layer δ 18 O data from site 851 reported an analytical uncertainty of 0.07‰, which we added in quadrature to the mixed-layer calculation ( 21 ). We applied an age uncertainty of 3ky, which is the estimated uncertainty of the core used for alignment and age control ( 22 ). For the interpolated data points at 152ka and 162ka, we incorporate additional uncertainty about our interpolated data. We estimate this additional uncertainty from the amplitude of the dominant period of the Site 851 data (obliquity, ~41ky) during the early portion of the record ( 21 ). The maximum amplitude of such change is ~1.6‰ from 4-300ky. The maximum offset from our intervals to Site 851 data is less than 20% of one obliquity cycle. We therefore estimate that the data would likely be within 20% of this maximum amplitude, or 0.33‰, and apply this additional uncertainty to these data points. We estimate the uncertainty in our variability data from our foraminifera populations and apply this to our weighted bivariate regression as well. MAD uncertainty is calculated as the standard error of the MAD. For the 94-98 th quantile mean, we our uncertainty as one half of the 2.5 and 97.5 quantile range, which approximates one standard deviation from our empirical data as above. We combine this with our calculated analytical uncertainty, in quadrature, to determine our total uncertainty. Climate mixed-layer temperature data set from 1958-2008. Amplitude reduction of 25% will be identified in over 90% of cases. Amplitude increase of 50% will likewise be identified in over 90% of cases. reference interval for all intervals is modern mixed-layer temperatures (1958-2009, 15-46m) from the 0.5°x0.5° cell surrounding the site of 17PC in the Central Equatorial Pacific ( 23 ).

Supplementary Tables
Supplementary Table 1