A Century of Reduced ENSO Variability During the Medieval Climate Anomaly

Climate model simulations of El Niño–Southern Oscillation (ENSO) behavior for the last millennium demonstrate interdecadal to centennial changes in ENSO variability that can arise purely from stochastic processes internal to the climate system. That said, the instrumental record of ENSO does not have the temporal coverage needed to capture the full range of natural ENSO variability observed in long, unforced climate model simulations. Here we demonstrate a probabilistic framework to quantify changes in ENSO variability via histograms and probability density functions using monthly instrumental and coral‐based sea surface temperature (SST) anomalies from 1900–2005 and 1051–1150 CE. We find that reconstructed SST anomalies from modern corals from the southwest Pacific capture changes in ENSO variability that are consistent with instrumental SST data from the central equatorial Pacific. Fossil coral records indicate 100 years of relatively lower ENSO variability during part of the Medieval Climate Anomaly. Our results demonstrate that periods of reduced ENSO variability can last a century, far longer in duration than modern observations in the instrumental record of ENSO, but consistent with results from unforced climate model simulations.


Introduction
The El Niño-Southern Oscillation (ENSO) is a coupled ocean-atmosphere climate phenomenon with global impacts on temperature and precipitation patterns (Bjerknes, 1969;Ropelewski & Halpert, 1987). ENSO is the leading mode of interannual (>1-9 year) climate variability, but instrumental observations are of insufficient length (Deser et al., 2010) to characterize the full range of natural variability (Wittenberg, 2009). Given the wide range of ENSO behavior simulated in the absence of forcings external to the climate system (Deser et al., 2012;Wittenberg, 2009), it is critical to ascribe the degree to which anthropogenic warming and internal climate variability are each contributing to future projections of ENSO in climate models (Collins et al., 2010;DiNezio et al., 2013). This motivates the use of paleo-ENSO reconstructions as out-of-sample tests of climate model simulations (Cobb et al., 2013;Gagan et al., 2000;Schmidt et al., 2014).
Isolating the unforced and forced components of ENSO variability remains an ongoing challenge in paleoclimate science, particularly for different mean climate states when forcings were different from today (e.g., the mid-Holocene or the Last Glacial Maximum) (Liu et al., 2014;Masson-Delmotte et al., 2013). Focusing on ENSO variability over the last 2,000 years (the Common Era, CE) provides context for the understanding of natural ENSO variability with current and near-future climate change. The Medieval Climate Anomaly (MCA: 950-1250 CE (Masson-Delmotte et al., 2013)) is identified as an interval with Northern Hemisphere surface temperatures similar to the modern (Masson-Delmotte et al., 2013), but our understanding of paleo-ENSO variability is inadequate due to a limited number of subannually resolved proxy records over the last millennium (Emile-Geay et al., 2013a, 2013b. Furthermore, given that the magnitude of orbital (Bertrand et al., 2002), solar (Bard et al., 2000), and volcanic (Crowley, 2000;Gao et al., 2008) forcing during the MCA is both small and similar to the modern (preindustrial), sustained changes in ENSO variability are likely dominated by processes internal to the climate system (McGregor et al., 2013;Rustic et al., 2015).
Coral records of surface ocean conditions extend our knowledge of interannual tropical climate variability to places and times when there is no (or limited) instrumental data (Corrège, 2006;Fairbanks et al., 1997;Gagan et al., 2000). Traditionally, coral-based ENSO reconstructions use the standard deviation of a band-pass-filtered time series (2-to 7-year window) of coral geochemical proxies as a metric of past ENSO variability (Cobb et al., 2003(Cobb et al., , 2013Emile-Geay et al., 2016;, but this approach (1) filters out, by mathematical construction, important ENSO variance that has a period of less than 2 years and (2) necessitates many decades and longer continuous datasets. Many fossil coral records, particularly older Holocene or Last Glacial Maximum corals, are short (several decades or less) or discontinuous and thus ill-suited for traditional filtering and data analysis methods (Cobb et al., 2013;Tudhope, 2001). To address these challenges, we extend the procedure suggested in Trenberth (1997) and use descriptive statistics in tandem with probability theory by assessing histograms (Trenberth, 1997) and probability density functions (PDFs) (Parzen, 1962) of monthly resolved coral data to quantify changes in ENSO variability.
The Niño 3.4 sea surface temperature (SST) index is a well-recognized record of ENSO variability (Trenberth et al., 2002); however, conditions in other regions of the tropical Pacific, notably the southwest Pacific, also accurately record changes in ENSO variability . Departures from the long-term monthly mean SST (SST anomalies; SSTA) from the Niño 3.4 region (5°N to 5°S, 120-170°W, Figure 1, box) in the central equatorial Pacific are canonically used to define the occurrence of ENSO events (Trenberth, 1997). During El Niño (La Niña) events, the central and eastern tropical Pacific experience positive (negative) SST anomalies that peak during boreal winter, while the western tropical Pacific experiences negative (positive) excursions (Trenberth, 1997) (Figure 1). Many paleo-ENSO studies target the Niño 3.4 region (Cobb et al., 2003;Nurhati et al., 2009), but other regions of the Pacific, like the tropical southwest Pacific in the South Pacific Convergence Zone, are also sensitive to ENSO variability, with ENSO detection skill broadly similar to the Niño 3.4 region (60-70% skill) (Hereid, Quinn, & Okumura, 2013). The tropical southwest Pacific is also advantageously home to abundant, high-quality modern and fossil corals, making this region a suitable location for paleo-ENSO studies (Gorman et al., 2012;Hereid, Quinn, Taylor, et al., 2013;Linsley et al., 2006;Quinn et al., 1996;Quinn et al., 2006). We also concentrate our efforts on reconstructing decadal to interdecadal changes in paleo-ENSO variability, rather than reconstructing the month-to-month changes of SST in the Niño 3.4 region, as this is difficult to reconstruct back in time due to age uncertainties (Emile-Geay et al., 2013a, 2013b. Here we use modern corals from Vanuatu, an archipelago in the southwest Pacific (Figure 1, star), to document ENSO variability during the twentieth century, and fossil corals to determine ENSO variability during the MCA. The recent decades of instrumental data  indicate that 71% of the variance of southwest Pacific SST anomalies is explained by ENSO. Despite the varied amplitude of the SST response to different types of ENSO events (Capotondi et al., 2015;Vincent et al., 2011), the southwest Pacific experiences a consistent SST response during ENSO events, albeit of the opposite sign to the Niño 3.4 region ( Figure S1): cooler SSTs during El Niño events ( Figure 1a) and warmer SSTs during La Niña events ( Figure 1b). The consistent SST response at Vanuatu during the most recent ENSO events increases our confidence in using coral records from the tropical, southwest Pacific for paleo-ENSO studies. Fieldwork at Vanuatu identified and recovered abundant, high-quality modern and fossil Porites lutea corals well suited for ENSO variability studies (section 2.1). Due to the tectonic activity of south Pacific islands like Vanuatu (Taylor et al., 1987), the rate of uplift outpaces sea-level rise, which exposes fossil corals above present day sea level. Another unique feature of our study site is that the fossil coral heads are all in situ (Thirumalai et al., 2015), allowing us to better understand the morphology of the reef flat and use estimates of the uplift rate to constrain the water depth in which the corals lived. We first demonstrate our data analysis technique by quantifying instrumental SST variability in the Niño 3.4 region and then apply our methods to replicated coral Sr/Ca-SST records from the southwest Pacific.

Coral Selection and Sampling
We located pristine, well-preserved, in situ fossil P. lutea coral heads spanning the last two millennia from an uplifted reef offshore of Tasmaloum, Vanuatu (TMV: 15.6°S, 166.9°E). The cores were drilled in 2011 using a Stihl chainsaw equipped with a Pomeroy Gear-Reduced Core Drill and diamond coring bits. All coral cores were uranium-thorium (U-Th) dated at the High-precision Mass Spectrometry and  1972-1973, 1982-1983-1998 Average NDJ SSTA for the 1988-1989, 1995-1996, 1998-1999 La Niña events. SST data are from the Met Office Hadley Centre HadISST product (Rayner et al., 2003). SSTA in this study are computed by applying a 9-year high-pass filter to monthly SST data, removing the climatology, and calculating the 5-month running mean (Trenberth, 1997) SSTA (section 3.3). The Niño 3.4 region (5°N to 5°S, 120-170°W) in the central equatorial Pacific is outlined by a white box (a, b). The modern coral site at Sabine Bank, Vanuatu, in the southwest Pacific (15.9°S, 166.0°E) is indicated with a star (a, b). Environment Change Laboratory, National Taiwan University, using multicollector inductively coupled plasma mass spectrometry Shen et al., 2012). Table S1 in the supporting information  provides a summary of the properties for the selected modern and fossil corals. Tables S2 and S3 provides the U-Th information for the 230 Th age calculation. This study uses fossil corals 11-TM-S5 ( 230 Th ± 2σ age: 1127.1 ± 2.7 CE, 36.4-cm depth) and 11-TM-I1 ( 230 Th ± 2σ ages: 1125.7 ± 6.2 CE, 30.5-cm depth; 1142.6 ± 4.9 CE, 18.0-cm depth; 1149.0 ± 4.1 CE, 10.7-cm depth). Based on estimates of the uplift rate (~5.5 mm/year) (Taylor et al., 1990), the selected fossil coral heads grew approximately 1-3 m below the sea surface during the MCA.
To provide modern climatological context, the analysis incorporated core 06-SB-A1 collected from a live P. lutea coral head at 8 m water depth at Sabine Bank, Vanuatu (SBV: 15.9°S, 166.0°E),~90 km to the southwest of TMV. Core 06-SB-A1 was collected in 2006 using the French Research Institute for Development vessel R/V Alis. Core 06-SB-A1 has been previously analyzed for oxygen isotopes (Gorman et al., 2012), but this is the first study to present Sr/Ca data from the same core. To quantify the replication uncertainty in our modern coral reconstructions, we also incorporated 50 years (1941( -1990 of published modern coral Sr/Ca data from Malo Channel, Vanuatu (MCV: 15.7°S, 167.2°E) (Kilbourne et al., 2004).
X-ray images of 5-mm slabs extracted from the coral cores ( Figure S2) highlighted the annual density banding and the optimal sampling paths along the maximum growth axis. All slabs were sonicated in distilled water and air dried prior to sampling. The coral slabs were micromilled at approximately monthly resolution (12 points per year) following established protocols (Alibert & McCulloch, 1997;DeLong et al., 2013;Fairbanks & Dodge, 1979;Marshall & McCulloch, 2002). The sampling resolution varied from 0.5-1.0 mm depending on the average growth rate of each respective coral (Table S1). The coral slabs were X-ray imaged a second time after micromilling ( Figure S2) to confirm that the sampling paths were parallel to the growth direction of individual corallites and along the central axis of a radially extending corallite fan (DeLong et al., 2013).
To develop a reliable Sr/Ca-SST record, it is critical to ensure that the coral is sampled along the maximum growth axis. We therefore considered how the coral growth architecture in three dimensions is projected in the 2-D plane of the cross-sectional slabs. In the case of fossil coral 11-TM-S5, we extracted additional 5-mm slabs from the core ( Figure S2) to generate a continuous record and ensure that the resultant Sr/Ca composite passed all of the quality control metrics outlined in DeLong et al. (2013). Sections with visible stress banding in the X-ray images and a lack of clearly defined theca walls (e.g., the bottom of the 11-TM-I1 replication fossil coral; Figures S2 and S4) are also suboptimal as this can impact the annual cycle in the geochemistry (Marshall & McCulloch, 2002) and/or make it difficult to identify the maximum growth axis for sampling. Suboptimal sampling can lead to unreliable climate reconstructions, so we excluded sampling paths that did not pass the quality control metrics of Delong et al. (2013) and conservatively limit our climate interpretations to the final Sr/Ca composites presented herein.

Coral Sr/Ca Analyses
Elemental ratio analyses were conducted using a Perkin Elmer Optima 8300 inductively coupled plasma-optical emission spectrometer located at UT Austin. All Sr/Ca measurements were corrected for plasma drift using standard-sample bracketing techniques (Schrag, 1999) with an internal reference solution gravimetrically prepared to have Ca, Sr, and Mg proportions similar to that of a coral. For each analysis, 113-262 μg of carbonate powder was dissolved in 2 wt. % nitric acid such that the Ca 2+ concentration in each sample was approximately 20 ppm, and within our 8-to 32-ppm calibration range for Ca 2+ .

Coral Sr/Ca Composites and Age-Modeled Time Series
X-ray images of the micromilled coral slabs ( Figure S2) provided clear constraints on the amount of overlap between two sampling paths, and the strong seasonal cycle observed in coral Sr/Ca was used to align peaks and troughs over the common period of overlap and generate the final Sr/Ca composite records. For a given year in the coral time series, the highest Sr/Ca value indicates the climatological coldest month, whereas the lowest Sr/Ca value indicates the climatological warmest month. To convert Sr/Ca versus depth to Sr/Ca versus time, we used a MATLAB® algorithm to identify the Sr/Ca peaks and troughs. The identified Sr/Ca maxima were assigned the climatological coldest month at Vanuatu (August), and the Sr/Ca minima were assigned the climatological warmest month (February). Once all of the annual peaks and troughs were identified, the Sr/Ca data were linearly interpolated to achieve monthly resolution (12 points per year). The relative age model for the modern SBV coral was converted to calendar years by counting back from the date of collection, whereas the relative chronologies for the fossil corals were converted to calendar years using four 230 Th ages as tie points (Tables S1-S3). The fossil coral Sr/Ca time series were shifted within the analytical error of the four 230 Th ages ( Figure S4) such that the resulting overlap between 11-TM-S5 and 11-TM-I1 achieved the highest Pearson correlation coefficient (Pearson, 1920) (r = 0.81, p < 0.01). We interpolated the published MCV modern coral Sr/Ca (Kilbourne et al., 2004) versus time data to 12 points per year using a piecewise cubic hermite interpolating polynomial (Fritsch & Carlson, 1980). We used a cubic interpolation scheme as it better preserved the sinusoidal shape of the original MCV modern coral Sr/Ca data, as compared to linear interpolation.

Instrumental SST Data
All instrumental SST data are from the Met Office Hadley Centre 1°latitude × 1°longitude gridded product (HadISST) (Rayner et al., 2003). SST data for the Niño 3.4 region was averaged over the (5°S to 5°N, 120-170°W ) domain1. The twentieth century historical ENSO events are based on the Oceanic Niño Index (NOAA: http://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php) and the multivariate ENSO index (Wolter & Timlin, 2011). For Vanuatu, we computed an unweighted average of the SST data from the two nearest grid points to the coral sites (15.5°S, 166.5°E) and (16.5°S, 166.5°E). The average of the two SST series yielded a resulting SST series that better represented the SST annual cycle observed in 7.75 years (November 1999 to July 2007) of in situ SST logger data from Sabine Bank, Vanuatu (Ballu et al., 2013), as compared to using only a single HadISST grid point for Vanuatu.

Proxy Calibration (Sr/Ca-SST)
We used the modern SBV coral Sr/Ca composite and gridded SST for Vanuatu to perform a calibration-verification exercise (Quinn & Sampson, 2002) ( Figure S3). We applied the following linear calibration to all age-modeled modern and fossil coral Sr/Ca measurements: SST (°C) = −20.73 × Coral Sr/Ca (mmol/mol) + 210.53. We also provide the inverse of the calibration equation to facilitate comparison with other coral geochemical studies: Coral Sr/Ca (mmol/mol) = −0.05 × SST (°C) + 10.16. To determine this calibration equation, we performed a weighted bivariate regression (Thirumalai et al., 2011) of the SST annual cycle (defined as the maximum SST − minimum SST) versus the SBV coral Sr/Ca annual cycle over the 1985-2005 CE calibration window. We performed a regression using the annual cycle to minimize age uncertainty on the monthly timescale. Although monthly resolution was targeted when sampling (section 2.1), one sample of coral powder may average 2-3 weeks (−2σ) of time when the coral is growing faster, or 5-6 weeks (+2σ) when the coral is growing slower. This yields some uncertainty in every "monthly" Sr/Ca value. A month-to-month calibration with instrumental SST would incorporate some of these errors. For this reason, we instead perform a regression using the annual cycle because we are most confident in identifying the clear peaks and troughs in the Sr/Ca data. This practice thus minimizes age uncertainty on the monthly timescale. When performing the weighted, bivariate regression, we conservatively used 0.1°C as the uncertainty in the instrumental SST, and the analytical uncertainty (0.012 mmol/ mol, ±2σ) as the uncertainty in coral Sr/Ca. The regression was force fit through the origin to yield a slope value of −20.73°C/mmol/mol. The y intercept for the calibration equation was empirically determined such that the median coral Sr/Ca-SST equaled the median instrumental SST over the calibration interval.
We used 1985-2005 CE as the calibration interval because this window includes the most recent 20 year of modern coral data that overlaps with the most reliable subset of the instrumental record that incorporates satellite observations (Rayner et al., 2003). The 1985-2005 CE calibration window maximized the correlation with instrumental SST (r = 0.87, p < 0.01) and also minimized the residual sum of squares over a 30-year verification window (1955( -1984 as compared to other possible calibration intervals. We note that four other potential calibration intervals (1950( -2005( , 1971( -2000( , 1980( -2005( , 1985 yielded similar slope values (range: −19.55 to −21.67°C/mmol/mol) that were all within the range of published coral Sr/Ca-SST calibrations (Corrège, 2006;DeLong et al., 2010). The coral trend of 0.83°C per 50 years at Sabine Bank ( Figure 4a) resembles the trend of 0.5-0.75°C per 50 years seen in observations for the southwest Pacific (Cravatte et al., 2009).

SST Data Processing
This subsection describes a series of mathematical operations that are used to isolate variability at interannual (>1-9 year) timescales, the result of which preserves more variance than a 2-to 7-year band-pass filter. We first applied a 9-year high-pass filter to the monthly instrumental and monthly coral-based SST to remove decadal and longer variability. We then removed the climatology to generate monthly SST anomalies calculated as deviations from the 1961-1990 CE climatology for the tropical Pacific composite map (Figure 1), the Niño 3.4 and Vanuatu instrumental SST, and the modern coral Sr/Ca-SST. Monthly SSTA were calculated as deviations from the 1126-1145 CE climatology for the MCA fossil corals. The climatology reference intervals were selected to maximize the temporal overlap between contemporaneous records. Lastly, we computed a 5-month running mean of monthly SSTA to smooth out intraseasonal variations (Trenberth, 1997). All PDF estimates of the monthly SST and monthly SSTA data were computed using a kernel density estimation method (Parzen, 1962). All instrumental and coral-based monthly SST anomaly results are reported as 9-year high-pass filtered, climatology removed, and 5-month running mean SSTA.

Uncertainty Analysis
We quantify changes in variability using the extreme percentiles (p 2.5 and p 97.5 ) of monthly SST and SSTA distributions (sections 5.1 and 5.2). We performed a Monte Carlo simulation (n = 10,000) to quantify the analytical and calibration uncertainties (Thirumalai et al., 2014) in the 2.5 (p 2.5 ) and 97.5 (p 97.5 ) percentiles for the coral-based monthly SST and monthly SSTA distributions. We also used replicated coral Sr/Ca records, to incorporate the effects of "geological" uncertainty, which is due to all other sources including the oceanographic setting and sampling. To calculate the uncertainties in the extreme percentiles, we first subset all contemporaneous coral Sr/Ca records to their common period of overlap prior to performing the Monte Carlo simulation (modern coral overlap: 1941-1990 CE; fossil coral overlap: 1126-1145 CE). As part of the Monte Carlo simulation, we perturbed each Sr/Ca data point in the original 50-(modern) or 20-year (fossil)-long time series n times with values randomly sampled from a normal distribution with mean 0 and ±2σ equal to the analytical uncertainty (±0.012 mmol/mol, ±2σ). We then transformed each Sr/Ca realization into SST taking uncertainty in the proxy calibration into account. For each realization of the Sr/Ca time series, a slope value was randomly sampled from a normal distribution centered on the empirically determined slope for Vanuatu (−20.73°C/mmol/mol; section 3.2) and a ±2σ range that approximately spanned the range of published coral Sr/Ca-SST calibration slopes (Corrège, 2006;DeLong et al., 2010). A corrective factor was applied to the y intercept such that the linear transformation for a given slope produced a y intercept that yielded the mean SST for the unperturbed time series. The resulting n SST realizations generated for each modern and fossil coral record include the impact of both analytical and calibration errors.
We next applied two filters to the coral records in order to remove long and short term variability and isolate interannual variability. First, a 9-year high-pass filter was applied to the n SST realizations prior to removing the climatology to remove variability on decadal and longer timescales. The 5-month running mean of the climatology-removed SST anomaly series were then computed to smooth out intraseasonal variations as defined in section 3.3 above. For each realization of the SST and SSTA time series we computed the p 2.5 and p 97.5 values. The overall uncertainty is the ±2σ range based on n realizations. In the event that the p 2.5 and p 97.5 ± 2σ values slightly differed for a given coral, we averaged the two values as the combined effect of analytical and calibration uncertainty. The analytical and calibration uncertainty quantification for the modern SBV coral is provided as an example ( Figure S5). We quantify the "geological" uncertainty in the coral records using replication. The total uncertainties that include the effect of replication in the 2.5 and 97.5 percentiles of the SSTA distributions are reported as the root-mean-square error of the percentile uncertainties for each respective modern and fossil coral determined by the Monte Carlo simulation with analytical and calibration uncertainty discussed above.
To test whether coral Sr/Ca-SSTA come from the same distribution, we performed Kolmogorov-Smirnov (K-S) tests (Massey, 1951) at the 1% significance level. Given that the SSTA time series are serially correlated, the effective degrees of freedom were considered when assessing the significance (Hu et al., 2017) of the K-S tests. Adjusting the effective degrees of freedom (ν eff ) makes it more difficult to reject the null hypothesis that the two data sets are from the same distribution at a specified significance level. The statistical significance of the K-S tests is further described and discussed in the results section. We also further explored the fidelity of the differences in the p 2.5 and p 97.5 values taking the total uncertainty into account. We used the results from our uncertainty analysis (e.g., Figure S5) to examine the overlap between the p 2.5 and p 97.5 distributions (additional details provided in sections 5.1 and 5.2 and Figures S6-S8).

SST Variability in the Niño 3.4 Region
We use instrumental SST data from the Niño 3.4 region (Rayner et al., 2003) (Figure 2a) as a test case to demonstrate how a probabilistic framework quantifies previously identified changes in twentieth century ENSO variability (Torrence & Compo, 1998;Trenberth, 1976). We choose a time domain subset of 1900-2005 CE to temporally match the modern coral climate record from Vanuatu. Prior research has often quantified ENSO variability using power spectra (Quinn et al., 1996;Wittenberg, 2009) and the standard deviation of band-pass-filtered SSTA (Cobb et al., 2003(Cobb et al., , 2013Emile-Geay et al., 2016). In this study, we use an alternative statistical approach, involving histograms and PDFs (Trenberth, 1997), which does not require continuous time series to characterize changes in ENSO variability. Another advantage of using these techniques is that they eliminate the need to identify discrete ENSO events, an ongoing challenge for paleo-ENSO records due to dating uncertainties (Comboul et al., 2014;Emile-Geay et al., 2013b;Hereid, Quinn, Taylor, et al., 2013). Instead, we characterize variability based on the distribution of observations over a time interval, an approach analogous to the analysis of individual foraminifera preserved in marine sediment (Thirumalai et al., 2013), albeit with more accurate annual chronology. The technique of looking at changes in ENSO over windows in the past has also been previously employed using corals of Holocene ages from the central Pacific (Cobb et al., 2013).  (Rayner et al., 2003) include the total variability in the SST record, including annual, interannual, and decadal and longer timescales. We note that for shorter records, linear or nonlinear detrending could be used in lieu of a high-pass filter to isolate interannual (subdecadal) variability. SSTs in the Niño 3.4 region document interdecadal variability in both the frequency and magnitude of ENSO events during the instrumental record, with an increase of extreme events over the last 40 years (Cai et al., 2014Trenberth & Hoar, 1996). For example, the 1982-1983 and 1997-1998 El Niño events (Figure 2c) are two out of the three most extreme ENSO events on record (Santoso et al., 2017). We quantify changes in ENSO variability during the twentieth century using two statistical metrics computed in moving windows (Okumura et al., 2017) (Figure 2e). Both the ±2 standard deviation range (±2σ) (Figure 2e, dashed line) as well as the difference between the 2.5 (p 2.5 ) and 97.5 (p 97.5 ) percentiles (Figure 2e, solid line) show lower values during the early twentieth century, and higher values during the late twentieth century that correspond with changes in the magnitude and/or frequency of ENSO events, that is, a change in ENSO variability. We note that if the data are normally distributed (Gaussian), the p 2.5 to p 97.5 interpercentile range is approximately equal to the ±2σ range.
These twentieth century changes in ENSO variability observed in the time domain (Figure 2e) are also captured using a probabilistic framework. We target intervals of enhanced and suppressed ENSO variability (Wittenberg, 2009) as end-members to demonstrate that histogram width quantifies changes in variability (Figure 2f). Extreme ENSO events yield SST anomalies that fall into the tails of the SSTA distribution as defined by the 2.5 and 97.5 percentiles (Figure 2d). A change in ENSO variability will correspondingly increase or decrease the p 2.5 to p 97.5 interpercentile range (herein referred to as the width of the distribution). The increase in ENSO variability during the late twentieth century extends the overall width of the SSTA PDF (Figure 2f, red PDF) as compared to the interval with less ENSO variability (Figure 2f, blue PDF).
Although both negative (La Niña) and positive (El Niño) SST anomalies become more extreme with increased ENSO variability, the relative increase in El Niño-related positive SST anomalies is larger, corroborating documented (Cai et al., 2014Trenberth, 1997) increases in the magnitude and frequency of strong El Niño events during the late twentieth century, that is, an increase in skewness.
The distribution of SST anomalies from the Niño 3.4 region demonstrates that histograms and PDFs quantify interdecadal changes in ENSO variability during the twentieth century. Although we use a continuous time series as demonstration, the histogram and PDF technique advantageously does not require long and/or continuous climate records. We subsequently apply the same statistical techniques to replicated SST records developed from modern and fossil corals from the tropical southwest Pacific.

Modern SST Variability
While it is common practice to use instrumental observations to characterize modern SST variability, observational coverage (Deser et al., 2010) in the southwest Pacific is limited in the first half of the twentieth century (Figure 3a), which leads to uncertainty in the magnitude of changes in ENSO-related SST variability. The SSTA signal in the observational product at Vanuatu (Figure 3d) is smaller compared to the Niño 3.4 region, but Vanuatu expectedly cools and warms during known historical ENSO events ( Figures S1, 1,  and 3d). However, instrumental SSTA at Vanuatu does not document a clear difference in interdecadal  : 1920-1939CE) and more (red: 1980-1999 shading in (c, e) highlight the intervals in (f) with less (more) ENSO variability. Numerical values above the PDFs in (d, f) denote the 2.5 and 97.5 percentiles for the designated subset interval. The horizontal bars above the PDFs in (d, f) indicate the p 97.5 -p 2.5 interpercentile range. Monthly SSTA calculated using the same methodology as Figure 1 (section 3.3). PDFs in this and all subsequent figures are based on a kernel density estimation method (Parzen, 1962). ENSO variability (Figure 3f) as observed in the Niño 3.4 region (Figure 2e). The PDF of instrumental monthly SSTA (Figure 3g) for Vanuatu shows an increase in width during the late twentieth century interval; however, the difference between the more and less variable intervals is small and not statistically significant (section 3.4). The SSTA distributions (Figures 3e and 3g) have a large concentration of weak anomalies, and the lack of interdecadal changes in ENSO variability is most likely due to the statistical infilling of the climatological mean (i.e., zero anomaly by construction in the data product) SST when observations are lacking (Rayner et al., 2003;Reynolds & Smith, 1994) (Figure 3a). Care must be taken when interpreting the variability gridded SST products for data-sparse regions, as interpolation, infilling, and other signal-processing techniques lead to loss of variance in these areas. Given a well-documented lack of variability in gridded SST products for observation-limited regions (Rayner et al., 2003), our modern coral-based  (Kennedy et al., 2011a(Kennedy et al., , 2011b grid box (17.5°S, 167.5°E) that includes Vanuatu. Note the logarithmic scale for the number of observations per month and the horizontal gray lines that indicate the number of observations required to achieve specified temporal coverages. The amplitude of SST variability outside of the monthly mean-removed climatology is loosely tied to the number of observations, such that more observations lead to more interannual, and even decadal, SST variability. (b) Instrumental monthly SST (Rayner et al., 2003) and (c) : 1920-1939CE) and more (red: 1980-1999 ENSO variability observed in the Niño 3.4 region (Figure 2). Monthly SSTA computed the same as Figure 1 (section 3.3). The x-axis in (e) and (g) is reversed to reflect the inverse relationship between SSTA at Vanuatu and SSTA in the Niño 3.4 region during ENSO events. The numerical values and horizontal bar above the PDF in (e) and (g) indicate the~±2σ range as defined by the percentiles.
SST reconstruction augments limited SST products and demonstrates the need for additional modern coral climate records from data-sparse regions.
The coral skeleton used to reconstruct and characterize modern monthly SST variability in the southwest Pacific was collected from a live P. lutea coral head at 8 m water depth at Sabine Bank, Vanuatu, during a trip in 2006 (Gorman et al., 2012) (15.9°S, 166.0°E; section 2.1). We apply the same analytical techniques used for instrumental SST to SST derived from coral Sr/Ca (Figure 4a) to test how corals from the southwest Pacific record changes in ENSO variability.
The Sabine Bank modern coral Sr/Ca-SSTA reconstruction (Figures 4 and S3) faithfully captures known interdecadal changes in ENSO variability over the last 100 years, as observed in the Niño 3.4 region. The SSTA estimates from the corals (Figure 4c), which are a point source, are larger than that for the spatially interpolated instrumental SST product (Figure 3d), which also contains a known loss of variance (Rayner et al., 2003). More importantly, the coral reconstruction agrees with the instrumental data that individual historical ENSO events alter SST in the region (Figure 4c). Another notable distinction is that unlike instrumental SST, the SSTA reconstruction from the coral archive faithfully captures known interdecadal changes in interannual variability over the twentieth century (Figure 4e). Reconstructed interannual SST variability at Vanuatu (Figure 4e) tracks the pattern of lower variability during the early twentieth century and a late twentieth century increase in variability as observed in the Niño 3.4 region (Figure 2e). This change is also captured by the PDFs for the select intervals with more and less ENSO variability (Figure 4f). La Niñarelated positive excursions become more extreme during the interval with enhanced ENSO variability (Figure 4f; p 97.5 less variable/more variable: 0.63 vs. 1.00°C). El Niño-related negative SST anomalies also become more extreme (p 2.5 less variable/more variable: −0.83 vs. −1.26°C). The larger values for the variability metrics in moving windows (Figure 4e) in conjunction with the increased width of the SSTA distribution for the interval with more ENSO variability (Figure 4f, red PDF) demonstrate that southwest Pacific modern coral Sr/Ca-SST captures the changes in ENSO variability observed in instrumental SST for the Niño 3.4 region (Figures 2e and 2f).
Our uncertainty quantification explores whether the Sr/Ca estimates of SST result in larger changes compared to instrumental data, and whether the changes in ENSO variability are statistically significant. The reported uncertainty in the 2.5 and 97.5 percentiles for modern coral SSTA is ±0.21°C for the monthly SSTA (Figures 4d and 4f) based on a Monte Carlo error propagation algorithm with analytical, calibration, and replication uncertainties (Thirumalai et al., 2014) (Figure S5 and section 3.4). Replication of the coral records quantifies the term that we refer to as "geological" uncertainty, which is due to all other sources, including the oceanographic and geologic setting and sampling. We quantify this geological uncertainty by comparing two modern corals over a common interval (1941( -1990 from the southwest Pacific (section 3.4). Our second, replication coral comes from Malo Channel, Vanuatu (Kilbourne et al., 2004) (15.7°S, 167.2°E),~120 km to the northeast of Sabine Bank. The ±2σ range for the Sabine Bank modern coral SSTA distribution  replication interval) is ±0.16°C including analytical and calibration uncertainty ( Figure S5). The ±2σ range for the Malo Channel modern coral SSTA distribution  replication interval) is ±0.14°C including analytical and calibration uncertainty. The resultant root-mean-square error is thus ±0.21°C as reported above. The modern coral SSTA distributions (Figure 6a) are reproducible over their common period of overlap and come from the same continuous distribution (passes the K-S test (Massey, 1951) at the 1% significance level regardless of the effective degrees of freedom; section 3.4). The intervals of more and less ENSO activity reconstructed from the Sabine Bank modern coral SSTA (Figure 4f) are significantly different at an effective degree of freedom 65.8% less (ν eff = 82) than the total number of months in each interval (n = 240; section 3.4).
We further explored the fidelity of the differences in the p 2.5 and p 97.5 values taking the total uncertainty into account. We used the results from our uncertainty analysis ( Figure S5) to examine the overlap between the p 2.5 and p 97.5 distributions for the intervals with more and less ENSO variability ( Figure S6). The amount of overlap between two percentile distributions highlights the similarity or difference between the reported percentile values. Our analysis confirms that the known changes in twentieth century ENSO variability are recorded by corals from the southwest Pacific, that the estimates are larger than instrumental estimates of SST changes, and that the changes in ENSO variability are large compared to the calculated uncertainty of the coral-based SSTA reconstruction.

MCA SST Variability
After demonstrating that modern corals from the southwest Pacific capture observed changes in ENSO variability, we next apply our statistical approach to corals from the MCA. The tectonic activity at Tasmaloum, Vanuatu (15.6°S, 166.9°E), yields pristine, well-preserved in situ fossil coral heads above present-day sea level. Our cores are collected from an uplifted reef, so we have the unique opportunity to sample multiple, contemporaneous coral heads and quantify the uncertainty in our fossil coral climate reconstruction via replication. The monthly SST anomalies for both fossil corals are calculated with respect to the 1126-1145 CE climatology since this interval is common to both corals. High precision U-Th dating Shen et al., 2012) confirms that the century-long fossil coral 11-TM-S5 ( 230 Th age: 1127.1 ± 2.7 CE, ±2σ) and the shorter, 24-year-long replication coral 11-TM-I1 ( 230 Th ages: 1125.7 ± 6.2, 1142.6 ± 4.9, 1149.0 ± 4.1 CE, ±2σ) selected for this study were alive~900 years ago during the MCA (section 2.1 and Table S1 ).
We apply our statistical techniques to the SST record derived from fossil coral Sr/Ca ( Figure 5a) as a test of ENSO variability during a century of the MCA (1051-1150 CE). The MCA fossil coral SST (Figure 5a) encompasses the total variability in the record (annual, interannual, as well as decadal and longer timescales), and shows similar overall variability (Figure 5b) compared to the modern period (Figure 4b). The SSTA time series (Figure 5c) shows both positive and negative excursions that correspond to ENSO events. However, the number of large SSTA excursions is smaller and leads to a narrower SSTA distribution (Figure 5d). While the modern SSTA shows a unidirectional increase in ENSO variability from the early to late twentieth century (Figure 4e), interannual SST variability during the MCA fluctuates between intervals with more and less variability (Figure 5e). For example, 1070-1090 CE has more interannual variability, while 1120-1140 CE has less variability as quantified by the ±2σ and interpercentile ranges. However, neither variability metric for the MCA (Figure 5e) exceeds the values for the last two decades of the twentieth century (Figure 4e).
We can compare the SSTA distribution for the MCA (Figure 5d) to either a century of modern coral data (Figure 4d) or discrete windows of time in the modern era (Figure 4f). The values for ENSO variability, as quantified by the percentiles of the SSTA distribution (Figure 5d), are 0.76°C for La Niña-related SSTA (p 97.5 ) and −0.82°C for El Niño-related SSTA (p 2.5 ), similar to what we observe during the earlier part of the twentieth century (Figure 4f). We choose to show an entire century of data for the MCA, but note that our results are consistent if we choose a subset of the fossil coral time series and generate the PDFs in 20-year moving windows ( Figure S9). The populations of SSTA for the two MCA fossil corals (Figure 6b) are drawn from the same distribution over their common interval of overlap, as they pass the K-S test (Massey, 1951) at the 1% significance level (regardless of the effective degrees of freedom; section 3.4). The uncertainty in extreme monthly SSTA values (Figure 5d, 2.5 and 97.5 percentiles) is ±0.24°C based on our algorithm with analytical, calibration and geological uncertainty (section 3.4). Incorporating the total uncertainty, we find that ENSO variability during the MCA, as recorded by the fossil corals, is within the range of ENSO variability observed in the modern and is significantly different than the interval with more ENSO variability during the late twentieth century ( Figure S7). Moreover, we cannot distinguish the SSTA population that records ENSO variability during the MCA as significantly different than the modern interval with less ENSO variability (based on K-S tests (Massey, 1951) at the 1% significance level and Figure S8). That said, we can distinguish the SSTA population during the MCA as significantly different than the SSTA population during modern interval with more ENSO variability (via the K-S test), and the new coral records from Vanuatu show less ENSO variability during the MCA as compared to the late twentieth century.

Discussion and Conclusions
Previously published proxy records of ENSO (Cobb et al., 2003;Moy et al., 2002;Newton et al., 2006;Rein et al., 2004;Rustic et al., 2015) for the MCA and the Little Ice Age (LIA: 1450-1850 CE (Masson-Delmotte et al., 2013)) often show conflicting results, indicating large uncertainties in the proxy records, in our estimations of the range of natural variability, or both. Numerous records provide evidence for a strengthened SST gradient across the equatorial Pacific and/or an inferred reduction in ENSO variability during the MCA relative to the LIA (Cobb et al., 2003;Newton et al., 2006;Rustic et al., 2015). In contrast, other proxy records (Moy et al., 2002;Rein et al., 2004) indicate a peak in ENSO variability during the MCA. Two recent compilations of ENSO-sensitive records actually find no statistically significant change in ENSO variability between the MCA and the LIA and highlight the need for additional high-resolution proxy records to fully characterize the range of ENSO variability over the Common Era (Emile-Geay et al., 2013b;Henke et al., 2017).
Internal climate variability contributes a large source of uncertainty in detecting forced changes in ENSO variability over the Common Era. Our results show a prolonged period of low variability during a time with external forcings similar to preindustrial values (Bradley et al., 2016). We interpret our results in tandem with the compilation studies to indicate that devoid of strong external climate forcing, internal variability within the climate system can produce a wide range of responses in the variability of ENSO. The PDF for coral data from the MCA indicates that ENSO variability over a full century is statistically indistinguishable from two decades with less ENSO variability observed during the twentieth century (Figures 7 and S8). We show a century of data for the MCA and 20 years of data for the modern, but our results are consistent if we choose a subset of the fossil coral time series and generate the PDFs in 20-year moving windows ( Figure S9). Furthermore, even when including the total uncertainty in the coral reconstructions (±2σ analytical, calibration, and geological uncertainty), ENSO variability during the MCA and the early twentieth century is statistically different and lower than the recent decades with larger ENSO variability ( Figures S6 and S7). Thus, we conclude that while the MCA contained lower ENSO variability, such ranges have been observed in the historical record.
Although our study focuses on ENSO variability during the Common Era, the histogram and PDF technique we present here is broadly applicable to other paleoclimate studies that seek to reconstruct variability across a variety of timescales. Quantifying the range of natural variability is critical as it may complicate our ability to detect a forced ENSO response from short records during times with different background states such as the mid-Holocene, the Last Glacial Maximum, future climate scenarios, and the most recent decades of instrumental data. Only by collecting more paleoclimate proxy data can we establish a baseline to determine if changes in ENSO variability during these other times are outside the bounds of natural variability. Our findings provide new insight to this challenge by replicating the bounds of low ENSO variability from two different time periods and showing that intervals of low ENSO variability can last for a full century, consistent with multidecadal to centennial intervals of reduced ENSO variability simulated in unforced climate models (Deser et al., 2012;Wittenberg, 2009).

Data Availability Statement
All coral Sr/Ca data from Sabine Bank and Tasmaloum, Vanuatu, produced from this study is archived in the paleoclimatology dataset repository in the National Centers for Environmental Information, NOAA database (https://www.ncdc.noaa.gov/paleo-search/study/28590).

Code Availability
The MATLAB® codes that have contributed to the analysis and results in this study are available upon request from the lead author (A. E. L: alawman@utexas.edu). The age model algorithm used to transform the coral geochemical data from the depth to the time domain is publicly available on the GitHub repository for the lead author (https://github.com/lawmana/coralPSM).