Unraveling forced responses of extreme El Niño variability over the Holocene.

Uncertainty surrounding the future response of El Niño-Southern Oscillation (ENSO) variability to anthropogenic warming necessitates the study of past ENSO sensitivity to substantial climate forcings over geological history. Here, we focus on the Holocene epoch and show that ENSO amplitude and frequency intensified over this period, driven by an increase in extreme El Niño events. Our study combines new climate model simulations, advances in coral proxy system modeling, and coral proxy data from the central tropical Pacific. Although the model diverges from the observed coral data regarding the exact magnitude of change, both indicate that modern ENSO variance eclipsed paleo-estimates over the Holocene, albeit against the backdrop of wide-ranging natural variability. Toward further constraining paleo-ENSO, our work underscores the need for multimodel investigations of additional Holocene intervals alongside more coral data from periods with larger climate forcing. Our findings implicate extreme El Niño events as an important rectifier of mean ENSO intensity.

INTRODUCTION

The response of the El Niño–Southern Oscillation (ENSO) to external climate forcing remains uncertain. Reducing this uncertainty is critical for producing reliable projections of future ENSO behavior and its impacts in a warming climate (–). The Holocene [11.65 thousand years ago to the present ()] provides an opportunity to explore the sensitivity of ENSO to changes in orbital forcing. The precession, or “wobble,” of the Earth’s rotational axis on ~21,000-year-long time scales influences the amount of incoming solar radiation (insolation) at a given latitude during a particular season (, ). ENSO events peak during boreal winter, and insolation-driven changes in seasonality (figs. S1 and S2) may influence the strength of ENSO (–). To investigate this hypothesis, paleoclimate records of past ENSO variability provide out-of-sample tests of climate model simulations of ENSO under changing Holocene background conditions (, ).

Isolating the forced ENSO response amidst a large range of natural variability remains a challenge, and thus, the underlying connection between ENSO variability and orbital forcing remains unclear. Many climate models that participate in the Paleoclimate Modeling Intercomparison Project [PMIP; ()] simulate a weak to modest reduction in ENSO variability (multimodel mean, ~8%; N = 30) 6000 years before present [kiloannum (ka)] during the mid-Holocene (~6 ka) compared to preindustrial conditions (, , ). Previous studies have attempted to develop a causal link between changes in seasonality and the amplitude of ENSO variability [see (, , , –) and references therein]. Nevertheless, whereas some paleoclimate proxy reconstructions suggest a reduction in ENSO variability during the mid-Holocene (–), other studies suggest a more heterogeneous tropical Pacific response (–). Previous work has applied model-data comparisons for the Holocene period (, ), although the combination of spatiotemporal variability in ENSO teleconnections () and the fact that orbital forcing is a slowly evolving, time-transient climatic driver necessitates different approaches to fully address the scope of externally forced Holocene ENSO changes. In addition, PMIP3-era mid-Holocene simulations are short [~100 years ()] and preclude the calculation of robust ENSO statistics to separate forced and unforced changes in ENSO (, ). Together, these challenges hinder our ability to confidently (i) simulate future changes in ENSO under a warming boundary condition and (ii) provide accurate comparisons with paleoclimate archives spanning the Holocene.

To clarify the response of ENSO to external forcing, we performed a systematic model-data comparison of the range of simulated and observed ENSO variability throughout the Holocene using paleoclimate proxy records from corals. We synthesize published Holocene coral oxygen isotope (δ18O) records from the Northern Line Islands in the central equatorial Pacific that span the most recent 7000 years of the Holocene (, , –) (Kiritimati and Fanning Atolls; Fig. 1). We focus on the central tropical Pacific due to the large ENSO signal, and because Kiritimati and Fanning together have the highest amount of publicly archived Holocene coral data. The oxygen isotopic composition of a coral’s carbonate skeleton is an established proxy that jointly reflects sea-surface temperature (SST) and the oxygen isotopic composition of seawater (, ), which is linearly related to sea-surface salinity (SSS) (). On interannual (year-to-year) time scales, El Niño events bring warmer and wetter conditions to Kiritimati and Fanning, which collectively yield negative coral δ18O anomalies (, , ). Opposite effects occur during La Niña, in which cooler and drier conditions yield positive δ18O anomalies.

Fig. 1.

CESM1 tropical Pacific SST anomalies.

(A) SD of tropical Pacific monthly SST anomalies (SSTAs) for 0 ka as simulated by CESM1 (Materials and Methods). (B to D) Simulated difference in the SD of SSTAs for (B) 3 ka, (C) 6 ka, and (D) 9 ka relative to 0 ka. Black box outlines the Niño 3.4 region (5°S to 5°N, 120°W to 170°W). Stars indicate the selected sites in the central equatorial Pacific for coral proxy data. Kiritimati (blue star; 2°N, 157°W) and Fanning (purple star; 4°N, 160°W).

We compare the coral-inferred estimates of past ENSO variability to new simulations for the intervals 0, 3, 6, and 9 ka performed with version 1 of the Community Earth System Model (CESM1; Materials and Methods and table S1). The long, multicentury simulations allow for the estimation of changes in ENSO driven by changes in external forcings between these intervals. The use of CESM1 is motivated by the model’s (i) realistic ENSO dynamics (–) and (ii) newly generated output for multiple time slices during the Holocene (figs. S2 and S3). We use a coral proxy system model (PSM) () that advances previous work (–) to simulate coral δ18O anomalies at Kiritimati and Fanning Atolls. The simulated estimates of coral δ18O variability are then directly compared to ~1200 years of published coral proxy data (tables S2 to S4) from these key sites (, , , , ). Our study enhances previous model-data comparison efforts for Holocene ENSO variability [see, e.g., (, , )] by using multicentury time slice simulations from a model with improved ENSO dynamics, and a process-based coral PSM to appropriately contextualize the coral data alongside simulated changes. The coral PSM includes various uncertainties and assumptions inherent to coral-based paleoclimatology, including variations in coral growth rates, analytical errors, and age modeling assumptions associated with developing a monthly resolved coral chronology (Materials and Methods). Output from the CESM1 Holocene simulations is processed with the coral PSM to provide a more accurate comparison with the coral proxy data and to assess model-data agreement for changes in ENSO variability. This is the first study to apply an intermediate-complexity coral PSM () to time slice simulations of extended Holocene paleoclimate output to quantify agreement on paleo-ENSO changes with published central Pacific coral geochemical data. Furthermore, we apply quantile-quantile (Q-Q) analysis to simulated and reconstructed paleo-ENSO changes to understand how El Niño extremes can influence mean ENSO intensity.

RESULTS

Simulated changes in Holocene ENSO variability

CESM1 simulates a modest yet consistent increase in ENSO variability (~2.1%/ka) from the early to late Holocene (Fig. 2). Using the 30-year running standard deviation (SD) of Niño 3.4 monthly SST anomalies (SSTAs) as a metric for capturing temporal changes in interannual variability (Fig. 2A and Materials and Methods), ENSO variability evolves from lower values during the early Holocene to higher values during the late Holocene. The median running SD values at 3, 6, and 9 ka respectively have an 8, 13, and 22% reduction in ENSO variability relative to 0 ka (Fig. 2A and fig. S4). Linear regression of the running SD values results in a statistically significant forced change in ENSO variability of ~2.1%/ka (P < 0.01), but with a large range of 30-year estimates.

Fig. 2.

CESM1 simulated changes in Holocene ENSO variability.

(A) The 30-year running SD of Niño 3.4 monthly SSTA for 9, 6, 3, and 0 ka. Lower and upper bounds of the boxes respectively correspond to the 25th and 75th percentiles, and the center line indicates the median (50th percentile). Whiskers represent the 1.5 × interquartile range, and outliers are indicated with a black dot. The range of the box and whiskers (Fig. 2A and fig. S4) captures intervals with higher and lower ENSO variability that arise purely from internal variability within the simulated climate system of each time interval. (B) Probability density functions (PDFs) of Niño 3.4 monthly SSTA for the Holocene time slices (see legend for labels). (C and D) Scatter plot of SD versus the number of extreme El Niño events in (C) 100-year and (D) 30-year windows. Extreme El Niño events are defined when the November-December-January (NDJ) average SSTA exceeds the 95th percentile (p95) of monthly SSTA for that time interval (Materials and Methods and fig. S5). Solid symbols in (C) indicate the average number of events/century for the full-length time slice simulations. Open symbols in (C) and (D) represent individual nonoverlapping windows. In this and all subsequent figures, the reported monthly anomalies are 9-year high-pass filtered, climatology-removed, and 5-month running mean anomalies to isolate interannual variability and facilitate comparison with coral proxy data.

The probability density functions (PDFs) of monthly SSTA for the Niño 3.4 region become progressively wider from 9 to 0 ka (Fig. 2B), further corroborating the increase in ENSO variability in Fig. 2A. A Kolmogorov-Smirnov (K-S) test () at the 5% significance level confirms that the SSTA distributions for 9, 6, and 3 ka are statistically different from those for 0 ka. In addition, the monthly Niño 3.4 SSTA distributions for all time slices are positively skewed with heavy warm tails (Fig. 2B), indicating a realistic asymmetry between the magnitude of warm SSTAs during El Niño and cool SSTAs during La Niña events (, ).

A change in the magnitude and frequency of extreme El Niño events (Fig. 2C and Materials and Methods) drives the increase in ENSO variability from 9 to 0 ka inferred from the change in SD (Fig. 2A). The SDs of SSTAs for each time slice are linearly correlated (R2 = 0.97) to the frequency of extreme El Niño events per century, illustrating how changes in extreme events lead to changes in overall ENSO variability (Fig. 2C). The linear relationship between SD and the frequency of extreme El Niño is consistent whether the overall SD (Fig. 2C, solid symbols) or the SD of nonoverlapping 100-year (Fig. 2C, open symbols) or 30-year windows (Fig. 2D) is used. Thus, the linear increase in the magnitude and the frequency of extreme El Niño events leads to the CESM1 simulated increase in ENSO variability over the past 9000 years.

The simulated changes in ENSO variability occur in conjunction with changes in the mean climate of the tropical Pacific compared to 0-ka conditions (fig. S2). Precessional forcing weakens the trade winds during February–May and deepens the thermocline in the western Pacific. This signal is communicated across the Pacific basin to the east, deepening the thermocline and warming SSTs during August–November in the eastern equatorial Pacific. Together, these changes lead to reduced seasonality in the eastern equatorial Pacific.

Coral comparison over the past 7000 years

The simulated coral δ18O anomalies (Materials and Methods and fig. S6) at both Kiritimati (Fig. 3, A and B) and Fanning (Fig. 3, C and D) show a modest linear increase in SD over the Holocene of similar magnitude to the simulated Niño 3.4 region (Fig. 2A). CESM1 has a more active ENSO compared to instrumental observations (, ) as well as some other documented model biases (), so the 30-year running SD of simulated coral δ18O anomalies in Fig. 3 is reported as a percent change relative to the median SD value for 0 ka to facilitate comparison with the coral proxy data.

Fig. 3.

Comparing simulated coral and coral proxy-inferred changes in Holocene ENSO variability.

(A and C) Percent change in SD for simulated coral (n = 100 realizations) and measured coral δ18O anomalies at (A) Kiritimati and (C) Fanning. Simulated coral δ18O anomalies generated using coral sensor model coupled with a process-based PSM that includes analytical errors, and uncertainties related to coral growth rates and age modeling assumptions (Materials and Methods) (). Percent change relative to the median value for 0 ka (simulated coral output) and 1 to 0 ka for the coral data, which include the modern coral composite record for each site (, , , ). Colored boxes show the 25th and 75th percentiles, and colored vertical lines indicate the 2.5th to 97.5th percentile range. Gray circles indicate the full range of simulated coral estimates. Colored circles indicate coral proxy records for each site, with full 30-year intervals outlined in black. Coral records <30 years are not outlined. (B and D) PDFs of monthly coral δ18O anomalies for (C) Kiritimati and (D) Fanning in which all the individual δ18O anomaly records (tables S2 to S4) are grouped into bins younger than 1 ka (blue, purple) and older than 1 ka (gray).

Whereas the trend in SD of simulated coral δ18O anomalies at the sites preserves the simulated ENSO signal over the Holocene, there is a wide range of unforced changes in ENSO variability on 30-year intervals according to CESM1—the typical duration of a coral record. For example, the SD of a 30-year interval at Kiritimati for 0 ka can range between −45 and +55% relative to the median value (Fig. 3A). This range is larger than the externally forced ~2%/ka Holocene trend that explains 20% of the variance in the SD of simulated coral δ18O anomalies at Kiritimati.

Many of the central Pacific coral δ18O observations show a greater reduction in SD compared to the simulated corals (Fig. 3, A and C), but large scatter in the coral data and large internal variability within the simulated corals make it difficult to definitively quantify the percent reduction in the coral data. Kiritimati and Fanning respectively have ~900 and ~330 years of modern and fossil coral data available for comparison with the CESM1 simulated corals (Materials and Methods and tables S2 to S4). Using Kiritimati as an example, a long 177-year coral record at ~4.3 ka () shows a reduction in interannual variability of ~55% compared to the 1- to 0-ka reference interval, whereas the difference in medians between the 0- and 6-ka time slices in the simulated corals is 13%. This ~55% reduction falls below the 2.5th percentile of the range of simulated coral estimates at 6 ka (Fig. 3A, vertical blue line versus gray circles). However, the scatter in the coral data close to 4.3 ka (5 to 3 ka) spans a range of ~50% (−17 to −70%) for changes in SD, and the total range of the change in SD at Kiritimati is ~90% (+19 to −70%). For context, the full range of internal variability for the 0-ka Kiritimati simulated corals can range from +53 to −45% relative to the median value, with 95% of the data falling within −28% to +25% (Fig. 3A).

A comparison of the distribution of monthly coral δ18O anomalies from 1 to 0 ka shows a wider distribution versus 7 to 1 ka at both sites (Fig. 3, B and D), supporting the simulated and proxy-inferred increase in SD during the late Holocene. As many of the coral records are short (<30 years) and contain a high degree of scatter between them, we group the records into those that are older than 1 ka, and those that are younger than 1 ka, including the modern coral record for each site (Fig. 3, B and D). A K-S test () at the 5% significance level confirms that the δ18O anomaly distributions for 1 to 0 ka versus 7 to 1 ka are statistically different.

Constraining Holocene ENSO behavior using Q-Q analysis

To further explore changes in ENSO variability, we use Q-Q analysis to highlight changes in the tails of our climatic distributions, i.e., extreme ENSO events (Materials and Methods and figs. S7 and S8). Warm El Niño events are more extreme than cool La Niña events such that temperature variations at ENSO-sensitive locations are not normally distributed (). Thus, parametric statistics (e.g., SD) may not fully capture changes in ENSO behavior and, moreover, may not sample the full range of ENSO signatures preserved in fossil coral records. Analysis of quantiles has thus recently expanded as a promising tool to track the full spectrum of seasonal and longer forms of variability using foraminifera (–). For example, previous work demonstrates the utility of Q-Q profiles for visualizing changes in ENSO variability using individual foraminiferal analysis at sites in the western and eastern equatorial Pacific (). Here, we introduce Q-Q analysis for coral paleoclimatology using the distribution of CESM1 monthly Niño 3.4 SSTA for 0 and 9 ka (Fig. 4A and fig. S8) as an illustrative example (see also fig. S7). Unlike the foraminiferal studies that compare geochemical populations between two time periods, our Q-Q analysis contrasts the coral data with standard normal quantiles to highlight departures from normality. This technique advantageously does not require continuous time series to investigate changes and eliminates the need to identify discrete ENSO events in the coral data (an ongoing challenge due to dating uncertainties). When the quantiles of the 0- and 9-ka SSTA distributions are compared with the equivalent quantiles for a standard normal distribution (Fig. 4B), the resulting Q-Q distribution is not linear and reveals that the differences occur in the tails of the distribution. At ENSO-sensitive locations, the degree of nonlinearity and curvature of the tails in the Q-Q distribution is shaped by the magnitude and frequency of extreme El Niño events. Higher SD for the 0-ka time interval (Fig. 2) yields a steeper slope compared to 9 ka in the Q-Q distribution (Fig. 4B), demonstrating that changes in ENSO variability affect both the tails and slope of the Q-Q distribution.

Fig. 4.

Simulated and coral-inferred changes in extreme ENSO events.

(A) PDFs of CESM1 simulated Niño 3.4 monthly SSTA for the 0-ka (dark red) and 9-ka (pink) intervals. (B to D) Q-Q plots for monthly (B) Niño 3.4 SSTA (0 and 9 ka). (C) Q-Q plots for coral δ18O anomalies younger than 1 ka (blue) and older than 1 ka (gray) at Kiritimati versus standard normal quantiles. (D) Same as in (C) for Fanning, with 1- to 0-ka corals in purple. (E) PDF of the δ18O anomalies for a 4.3-ka fossil coral (blue) from Kiritimati () and nonoverlapping 30-year windows (gray). PDF of modern coral δ18O from Kiritimati (red). (F) Q-Q plots for the 4.3-ka (blue) and modern (red) coral δ18O anomalies versus standard normal quantiles. (B to D and F) Solid gray line connects the first and third quartiles of the data. Dashed gray reference line extends the solid interquartile line to the ends of the data. Each data set in (B) to (D) and (F) is divided into 100 quantiles. If the distribution of SSTA or δ18O anomalies is normal, then the Q-Q distribution will be linear. A change in the distribution (e.g., the SD and skewness) will change the slope of the Q-Q distribution. Please see fig. S7 for an additional Q-Q plot interpretation guide.

Grouped coral δ18O anomalies at Kiritimati (Fig. 4C) and Fanning (Fig. 4D) are consistent with Q-Q profiles from the simulated Niño 3.4 SSTA (Fig. 4B). The older time intervals contain reduced slopes and less pronounced departures from linearity in the tails, in accordance with reduced ENSO variability in the past. The PDFs at both sites are negatively skewed with heavy tails (Fig. 3, B and D) due to the occurrence of warmer and/or less saline sea-surface conditions during extreme El Niño events. Whereas the difference in slopes is consistent with reduced SDs and lower ENSO amplitude, Q-Q patterns in the negative tails of both sites (Fig. 4, C and D) provide evidence for extreme El Niño events occurring in the past. A closer examination of the distribution of monthly δ18O anomalies from the Kiritimati modern and long 177-year fossil coral record from ~4.3 ka () (Fig. 4, E and F) corroborates this finding using individual records. The distribution of monthly coral δ18O anomalies at 4.3 ka is notably narrower than the modern record at Kiritimati (Fig. 4E) and is more linear in the Q-Q distribution (Fig. 4F). However, similarly to the modern record, long negative tails are also found for two of the five nonoverlapping 30-year windows (Fig. 4E), which, along with the character of the PDFs, suggests that extreme El Niño events still occurred in the past.

DISCUSSION

This work presents a systematic model-data comparison that incorporates uncertainties related to the coral archive and provides constraints on the range of ENSO variability over the Holocene epoch. Together, the Kiritimati and Fanning coral δ18O records provide an independent test of the simulated changes in Holocene ENSO variability by CESM1. The simulated monotonic increase in ENSO variability from the early to late Holocene (Fig. 2A) is driven by a change in the magnitude and frequency of extreme El Niño events (Fig. 2C). The CESM1 time slices show a linear relationship between SD and the number of extreme events (Fig. 2C) that is further supported by diminished nonlinearity in the Q-Q plot at 9 ka versus 0 ka (Fig. 4B).

We hypothesize that the simulated reduction in extreme El Niño events is caused by seasonal changes in the mean climate of the equatorial Pacific in response to precessional forcing (, ). CESM1 simulates a seasonal reduction in the strength of the trade winds in the western Pacific during early to late spring under 9-, 6-, and 3-ka boundary conditions (February–May; fig. S2). In response, the thermocline deepens, a signal that propagates across the basin and results in a deeper thermocline in the east and reduced stratification during boreal summer (fig. S2). The change in the thermocline, in turn, contributes to warmer eastern Pacific SSTs during boreal fall and a weaker cold tongue. On the basis of previously documented mechanisms (), these changes weaken the upwelling feedback in the eastern Pacific and contribute to a reduction in extreme El Niño events.

The comparison between measured coral data and PSM-processed simulations (Fig. 3) shows that both lines of evidence point to intervals with reduced ENSO variability during the Holocene, a result consistent for both the SD (a traditional metric for ENSO amplitude) and the Q-Q analysis. Similar to CESM1, the grouped 7- to 1-ka fossil coral data at Kiritimati and Fanning both have a lower SD compared to the 1- to 0-ka interval (Fig. 3) and reduced nonlinearity in their Q-Q profiles (Fig. 4, C and D), corresponding to fewer and/or weaker extreme events. This result is also consistent for the modern and long 4.3-ka fossil coral records from Kiritimati (Fig. 4F).

However, despite the similarities between the CESM1 time slices and the coral proxy data, large differences exist between simulated and reconstructed magnitudes of the Holocene changes in ENSO. These persist even when we (i) consider multicentury simulations spanning the Holocene and (ii) use an intermediate-complexity PSM to place corals and the model in the same reference frame. The median percent change in simulated coral SD for the 3- and 6-ka time slices compared to 0 ka is far less (−5.1 and −10.2%, respectively) than the difference between the median of the Kiritimati fossil coral data within the 5- to 3-ka interval (−50.9%).

It is unclear whether the discrepancy in magnitude between simulated corals and measured corals over the past 7000 years stems from shortcomings in the model, or from a lack of coral data, as the coral data are only a subset of true ENSO variability with a relatively low signal-to-noise ratio. Potential sources of uncertainty capable of producing the discrepancy between the simulated coral and measured coral data over the past 7000 years are further discussed below.

Climate model and PSM uncertainties

The time slice simulations used in this work are from a single model such that model biases or transient climate phenomenon may hinder our comparison with the coral data. Model biases can lead to divergent projections of future ENSO, and these same biases could contribute to divergent simulations of ENSO variability in the past (). CESM1 simulates an ENSO that is too active compared to observations (), and this, in combination with other biases (), could drive divergent trends in SD between the simulated and observed Holocene coral data. In addition, within the time slice simulations, centennial- to millennial-scale variability associated with volcanic aerosols, changes in solar output, or other feedbacks within the climate system that operate on longer time scales are not included. For example, documented millennial-scale ENSO variability reconstructed in sediment records (, , ) would also affect coral archives and thus lead to discrepancies between the simulations and coral data.

In addition, there may be unknown structural uncertainties, for example, nonclimatic processes, or assumptions regarding data extraction that are not accounted for in the coral PSM. Such assumptions could systematically affect coral-inferred estimates of interannual variability and contribute to the simulated coral-measured coral data disagreement. Future process-based, in situ studies will help elucidate coral response to nontemperature influences, such as hydroclimate. For example, recent advances in isotope-enabled modeling will reduce uncertainty in the [δ18O of seawater—salinity] relationship used in our PSM and allow us to better understand the temperature versus hydrological response to ENSO at coral proxy locations.

Coral paleoclimate data uncertainties

The largest limitation of the data compilation of central Pacific corals is that, despite more than three decades of critical data collection, there still are insufficient data spanning the Holocene to quantify the forced ENSO signal amidst a wide range of natural ENSO variability. Although CESM1 simulates a statistically significant trend in ENSO variability from 9 to 0 ka (2%/ka, P < 0.01), there is large internal variability within the individual time slices (Fig. 2A) that may hinder the ability to detect a small forced change in ENSO variability using short (several decades or less) coral records. The SD of coral δ18O anomalies (Fig. 3) yields a heterogeneous picture, a result consistent with previous coral compilation studies (, ). We interpret millennia with a limited number of years of coral data (e.g., 3 to 1 ka) with caution, as there are likely insufficient data to capture the full range of natural variability; therefore, the data may not converge on the median of the forced signal. There are, however, millennia with sufficient coral data from a single site (e.g., 5 to 3 ka at Kiritimati) that consistently show reduced interannual variability compared to the most recent 1000 years. Reproducing and extending this interval of reduced variability at Fanning Atoll and at other ENSO-sensitive locations will provide important validation and help further constrain the magnitude of the percent reduction. While there is better model-data agreement at Fanning (Fig. 3C) compared to Kiritimati (Fig. 3A), additional coral records will allow us to further investigate these site dependencies.

Although this study focuses on the central Pacific where the amount of fossil coral data is the largest for model-data comparison purposes, coral data from other ENSO-sensitive sites will also advantageously enable future studies that seek to investigate changes in the spatial pattern of ENSO during the Holocene. CESM1 simulates the largest reductions in interannual variability in the eastern Pacific (Fig. 1, B to D) and an associated decrease in extreme El Niño events (Fig. 2, C and D), suggesting that there may be changes in the spatial pattern of ENSO, consistent with previous modeling work (). A recent synthesis of coral and mollusk data from the tropical Pacific tentatively suggests that there may be changes in the spatial pattern of ENSO (). However, at present, it is difficult to validate climate-model inferred changes in the spatial pattern of ENSO due to an insufficient amount of proxy data.

Acquisition and measurement of additional fossil coral data spanning the past 2000 years, as well as corals older than 7 ka, are also necessary to overcome uncertainties in our trend analysis. A lack of coral data for the early Holocene prevents us from conclusively establishing whether there is an early Holocene increase in ENSO variability relative to a local minimum during the mid-Holocene, or a monotonic increase in ENSO variability throughout the Holocene as simulated by CESM1 and inferred from other proxy records [see, e.g., ()]. Furthermore, the coral-inferred trend in ENSO variability over the past 7000 years at Kiritimati is sensitive to the choice of baseline and the inclusion or exclusion of the modern or short (<20-year-long) records (table S5). Coral data predating 7 ka from these and other ENSO-sensitive sites would also provide insights surrounding early Holocene ENSO variability with which to compare to the 9-ka time interval when the CESM1 simulated corals suggest larger forced changes to ENSO.

Last, our analysis underscores the need for fossil corals that are at least 30 years long, and ideally more than 100 years long. Uncertainties surrounding the removal of the climatology from a short record may lead to uncertainties in the resulting δ18O anomalies. Furthermore, a short coral record may fall within an active or inactive interval of natural variability, in which the SD of a single coral record is not a good estimate for the overall SD of the past time interval. Short records also prevent us from establishing a threshold for extreme El Niño events using the internal properties of each coral δ18O anomaly series as we are able to do with CESM1 (Fig. 2, C and D) and with long coral records. Long fossil coral records will reduce uncertainty in the thresholds for extreme events, although we acknowledge the great complexities and challenges that the collection of such samples entails.

Implications for paleo-ENSO

The ability to characterize past and future changes in ENSO variability relies on studies that incorporate information from both climate models and paleoclimate proxy records. Previous work has attributed changes in Holocene ENSO variance to changes in seasonality that arise from orbital forcing, with an emphasis on the importance of boreal summer and fall insolation (, , –). Many previous modeling studies have exclusively focused on the mid-Holocene and preindustrial time intervals. However, many of these simulations are less than 250 years long, making it difficult to separate out internal ENSO variability from the forced response.

Thus, the multicentury CESM time slice simulations presented in this work address this limitation, and the comprehensive model-data comparison shows that reduced ENSO variability in both CESM1 and the central Pacific coral proxy data during the Holocene is driven by a change in extreme El Niño events. Furthermore, the time slice simulations provide an important perspective on long-term changes in ENSO variability that may arise due to changes in external forcing. Our study also enhances previous model–coral data comparison efforts for Holocene ENSO variability (, , ) by using a process-based coral PSM in conjunction with the multicentury CESM1 time slice simulations. The PSM allows us to directly compare the CESM time slices and the coral geochemical data on a level plane (Fig. 3). Together, this represents a comprehensive model-data comparison of Holocene ENSO. Both CESM1 and the central Pacific coral proxy data show intervals with reduced ENSO variability during the Holocene driven by extreme El Niño events. However, too few long coral records and uncertainties from internal variability do not allow us to quantify the relationship between SD and extreme El Niño events in the coral data over the past 7000 years.

In closing, the model simulations and coral PSM framework presented here provide context for characterizing the range of natural ENSO variability, as well as how ENSO responds to changes in external forcing during different background climate states. Future studies involving time slice or transient simulations from multiple models will provide an opportunity to investigate the role of differences in model physics on the simulation of ENSO and explore the response of ENSO to changes in external forcing through time. Furthermore, ongoing coral data collection efforts from ENSO-sensitive locations will provide expanded constraints on ENSO variability during times with different climate forcings. As anthropogenic greenhouse gas emissions continue to warm the Earth at an accelerating pace, rigorous data synthesis and model-data comparison studies are needed to (i) validate and improve climate model simulations and (ii) help reduce uncertainties in future climate projections.

MATERIALS AND METHODS

Holocene climate model simulations

The climate model simulations used in this study are from version 1.2.1 of the Community Earth System Model, which uses Community Atmospheric Model Version 5, CESM1 (CAM5) (). CESM1 has ~1° horizontal resolution in the ocean and the atmosphere (). Our analyses use surface temperature and salinity output for four time intervals spanning the Holocene (9, 6, 3, and 0 ka). Each simulation was run with the boundary conditions corresponding to the interval, including changes in the distribution of seasonal and latitudinal insolation, greenhouse gas concentrations, and sea level.

The 1500-year-long simulation for 0 ka is analogous to a preindustrial long control simulation (), as outlined in the fifth and sixth phases of the Coupled Model Intercomparison Project [CMIP5/CMIP6; (, )]. The 0-ka simulation was run with external forcings (solar irradiance, orbital, greenhouse gas, dust and other aerosol emissions, ozone, and land use) set at constant preindustrial (1850 CE) values. The 0-ka simulation was initialized from an existing simulation () and run until the climate reached a new equilibrium. All analyses using the 0-ka simulation were performed using output from the past 1500 years of the equilibrated climate.

Each of the remaining Holocene simulations was initialized from the 0-ka simulation and run until the surface climate reached equilibrium. All analyses using the 9-, 6-, and 3-ka simulations were performed using output from the past 500 years of the equilibrated climate. All simulations are of sufficient length to capture a range of ENSO behavior ().

Changes in precession-related insolation variability (, ) are the dominant forcing over the Holocene simulations (fig. S1), although greenhouse gas concentrations change by small amounts (table S1). The 9-ka simulation includes a 20-m decrease in global mean sea level, but the impact of sea-level changes at 9 ka is minimal in the context of large changes in orbital forcing. The concentrations of greenhouse gases (CO2, CH4, and N2O; table S1) were altered for each simulation based on ice core measurements (). The concentrations of other greenhouse gases, such as chlorofluorocarbons, were set to zero. Ozone concentrations were prescribed at 1850 values for all simulations. CESM1 was run with its prognostic aerosol module, which requires prescribed aerosol and dust emissions. The 9-, 6-, and 3-ka simulations were run with the same dust emissions as in the 0-ka simulation. The solar flux was kept constant at 1365 W/m2 in all simulations. In each simulation, vegetation was also prescribed to be the same as in the 0-ka control simulation. The remaining surface properties, such as albedo or surface roughness, are computed by the land component of CESM1 based on the soil and plant properties and passed to the atmosphere component via the coupler. The vegetation phenology, including the total leaf and stem area indices and canopy heights, was prescribed and does not respond to changes in climate.

SSTAs averaged across the Niño 3.4 region in the central equatorial Pacific (5°N to 5°S, 120°W to 170°W) are canonically used to define the occurrence of ENSO events () and quantify the range of ENSO variability. The simulated changes in the SD of monthly Niño 3.4 SSTA (figs. S3 and S4) are used to quantify changes in ENSO variability in this study. To investigate changes in interannual variability at sites with coral proxy data, we use CESM1 SST and SSS output for the regions encompassing Kiritimati Atoll (2°N, 157°W) and Fanning Atoll (4°N, 160°W) in the central equatorial Pacific.

Central equatorial Pacific coral geochemical data

We compile published monthly resolved modern and fossil coral oxygen isotope data from Kiritimati and Fanning Atolls in the central equatorial Pacific (tables S2 to S4) (, , , , , –). The monthly 18O/16O variations are reported as δ18O in per mil (‰) relative to Vienna Pee Dee Belemnite (VPDB). Data from Palmyra Atoll are excluded as there are insufficient coral data earlier than the last millennium with which to compare with the CESM1 Holocene time slices.

Kiritimati and Fanning are both located within the Niño 3.4 region and collectively have around 1200 years of published coral proxy data spanning the most recent 7000 years of the Holocene. Kiritimati and Fanning respectively have 901 and 333 years of fossil coral data. Previous work by (, ) demonstrates the high correlation between modern coral δ18O at these sites and the Niño 3.4 SST Index on interannual time scales [−0.92 for Kiritimati and −0.85 for Fanning ()]. The selected central Pacific coral data and their corresponding 230Th ages are publicly archived by the National Centers for Environmental Information (www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/coral-sclerosponge) and organized as MATLAB files on GitHub (https://github.com/CommonClimate/EmileGeay_NatGeo2015) and National Oceanic and Atmospheric Administration (www.ncdc.noaa.gov/paleo-search/study/22415) as part of previous paleo-ENSO studies (, ). The majority of the fossil coral records are short, ranging between 7 and 82 years in length (tables S2 and S3), although there is one 177-year-long coral δ18O record from ~4.3 ka from Kiritimati (). This study also includes published modern coral δ18O records generated from living corals for both Kiritimati (1938–2016) (, , , ) and Fanning (1949–2005) (). Although the responses at Kiritimati and Fanning are similar, with both sites tracking the changes observed for the broader Niño 3.4 region (Fig. 2A versus Fig. 3), the two locations are analyzed separately to avoid any site-specific responses. This approach also facilitates a more direct comparison between the local simulated corals and the observed proxy data at a specific site.

Modeling simulated coral δ18O

To investigate how large-scale changes in simulated ENSO variability manifest at specific sites, we use a process-based coral PSM framework developed by (, , ) to transform the CESM1 Holocene simulations into estimates of past ENSO variability that are more directly comparable to coral geochemical data. The surface temperature and SSS output for the model grid points closest to Kiritimati and Fanning Atolls are transformed into mean-removed (i.e., zero-centered, Δ) simulated δ18O anomalies for each Holocene time slice using the multivariate coral sensor model of ()

We use a slope −0.22‰/°C for a1, based on the inverse SST dependence that arises from thermodynamic fractionation (). SSS and the δ18O of seawater (δ18Osw) are often assumed to be linearly proportional and are both affected by precipitation, evaporation, and advection (). We approximate a2 using the observed δ18Osw-SSS slope of 0.27‰/practical salinity unit determined from basin-scale regression analysis for the tropical Pacific (). The resulting Δδ18O values are then subsequently converted from the Vienna Standard Mean Ocean Water (VSMOW) reference standard to the VPDB reference for carbonates using the following equation: δ18OVPDB = 0.97001 × δ18OVSMOW ().

The resulting simulated coral Δδ18O time series for 9, 6, 3, and 0 ka for Kiritimati and Fanning are then subsequently perturbed (n = 100 realizations) by a process-based coral PSM () that includes uncertainties inherent to coral-based paleoclimatology that may alter the ENSO signal of interest. The coral PSM incorporates the impact of (i) variable growth rates experienced when sampling a coral along the maximum growth axis, (ii) analytical errors typical of laboratory analytical precision (0.1‰, ±1σ), and (iii) the age modeling assumptions associated with transforming coral geochemical data from the depth to the time domain. The combination of the three algorithms listed above is herein referred to as the full coral PSM. We direct the reader to () for an in-depth discussion of the coral PSM algorithms and how various nonclimatic processes and assumptions may affect a simulated coral’s ability to capture changes in ENSO-related variance.

In summary, the growth rate algorithm incorporates the effect of how variable annual coral growth rates affect a coral geochemical time series when the coral archive is sampled at a fixed resolution (e.g., 1-mm increments) along the maximum growth axis. The age model algorithm transforms coral geochemical data from the depth to the time domain. For coral δ18O, temperature and salinity variations often constructively interfere such that more negative extrema indicate warmer and/or less saline conditions, while more positive extrema indicate cooler and/or more saline conditions (, ), although exceptions may occur. When constructing an age model, the peaks and/or troughs in the coral geochemical data are assigned specific calendar months based on knowledge about the climatology at the site. Once the annual geochemical extrema are identified, the coral data are interpolated to achieve evenly spaced monthly resolution. Previous work by () developed a MATLAB algorithm to standardize coral age modeling. This version identified both the annual peaks and troughs in the geochemical data before linearly interpolating to monthly resolution. As part of this study, we added additional functionality to the age model to optionally use a single chronological tie point, as is commonly done for sites with a small annual cycle. Here, we use a one chronological tie point per year, where the most negative δ18O value every year is assigned a calendar month based on the simulated climatology of the site. The PSM algorithms used in this study, including the updated age model algorithm, are publicly available to the broader community via a GitHub repository: https://github.com/lawmana/coralPSM.

The simulated change in Holocene ENSO variability is consistently evident at each site in (i) the SST output from the model, (ii) a simple bivariate coral sensor model using SST and SSS to model δ18O (, ), and (iii) the full coral PSM using SST and SSS to model simulated coral δ18O (fig. S6) (). There is an expected reduction in interannual variability as quantified by the shift in SD values between the original and perturbed model output (fig. S6). This loss of variance is predominantly attributed to the age modeling assumptions. Previous work shows that the impact of age modeling assumptions is particularly pronounced in the central equatorial Pacific, where the amplitude of the annual cycle is small and the impacts from ENSO are large (). Furthermore, the age modeling subcomponent of the coral PSM is particularly relevant for studies of Holocene ENSO variability using corals given that orbital precession changes seasonality (fig. S2), and this may affect the variability of simulated coral δ18O at Kiritimati and Fanning Atolls.

Statistical analysis: Quantifying changes in past ENSO variability

We perform a series of mathematical operations to isolate interannual (>1- to 9-year) variability in the coral δ18O records and the CESM time slices, the result of which preserves more variance compared to a traditional 2- to 7-year or 2- to 8-year band-pass filter (). The SST output for the Niño 3.4 region (fig. S3) and the simulated coral δ18O time series for Kiritimati and Fanning are first filtered with a 9-year high-pass filter to remove variability at decadal and longer time scales. Second, the climatology for the full-length time series is removed to generate monthly SST or δ18O anomalies. Third, a 5-month running mean is applied to the monthly anomaly series to smooth out subseasonal variability (i.e., “weather”) (). To facilitate an objective model-data comparison, we apply the same mathematical operations to all coral δ18O records from Kiritimati and Fanning Atolls.

To quantify the range of variability, we calculate the 30-year running SD of Niño 3.4 SST, simulated coral δ18O anomalies, and measured coral δ18O anomalies. A 30-year window is used to reflect the average duration of a coral record. We report the SD of the full-length time series if a coral record is <30 years long. The running SD of the climatology-removed anomalies is presented as a more robust metric for capturing temporal changes in interannual variability and has been previously used for fossil coral records spanning the mid- to late-Holocene (, , ). In summary, the overarching objective is to process both the model output and coral data the same way such that we (i) isolate interannual variance, (ii) quantify the range of natural variability, and (iii) develop a technique that is equally well equipped to handle short and/or discontinuous fossil coral records.

For each CESM time slice simulation, the occurrence of an extreme El Niño event (Fig. 2, C and D) is defined when the average November-December-January (NDJ) SSTA in the Niño 3.4 region exceeds the 95th percentile (p95) of monthly SSTA for that time interval (other thresholds are provided in fig. S5 for comparison). The 95th percentile thresholds for extreme El Niño events, which contribute to the overall change in simulated ENSO variability, are 1.9°, 1.8°, 1.8°, and 1.4°C for 0, 3, 6, and 9 ka, respectively (fig. S3). The use of a nonparametric threshold for extreme El Niño avoids any assumption of normality, which is appropriate given the skewed SSTA distributions (Fig. 2B). Given that the 0-ka simulation is three times longer (1500 years) than the other time intervals (500 years), we report the average number of extreme events per century (Fig. 2C, solid symbols).

All K-S tests () are reported at the 5% significance level. Trend analysis is performed using ordinary least squares linear regression for CESM1 and bivariate weighted linear regression () for the coral data. Errors are taken to be 0.05‰ for each SD value with an error of 1/N, where N is the length of the window. For longer records in which we compute the 30-year running SD, the weighting is 1/30 for each SD value. For <30-year records, the weighting is proportional to the length of the full record such that longer records are given greater weight in the trend analysis.

Q-Q analysis is used to visualize differences in the Niño 3.4 and coral δ18O anomaly distributions compared to a standard normal distribution (Fig. 4). We standardize the number of quantiles (N = 100 quantiles) in all comparisons. An additional Q-Q plot interpretation guide is provided in fig. S7.