Inferred changes in El Niño-Southern Oscillation variance over the past six centuries

It is vital to understand how the El Niño–Southern Oscillation (ENSO) has responded to past changes in natural and anthropogenic forcings, in order to better understand and predict its response to future greenhouse warming. To date, however, the instrumental record is too brief to fully characterize natural ENSO variability, while large 5 discrepancies exist amongst paleo-proxy reconstructions of ENSO. These paleo-proxy reconstructions have typically attempted to reconstruct the full temporal variability of ENSO, rather than focusing simply on its variance. Here a new approach is developed that synthesizes the information on common low frequency variance changes from various proxy datasets to obtain estimates of ENSO variance. The method is tested using 10 surrogate data from two coupled general circulation model (CGCM) simulations. It is shown that in the presence of dating uncertainties, synthesizing variance information provides a more robust estimate of ENSO variance than synthesizing the raw data than identifying its running variance. We also examine whether good temporal correspondence between proxy data and instrumental ENSO records implies a good represen- tation of ENSO variance. A signiﬁcant improvement in reconstructing ENSO variance changes is found when combining several proxies from diverse ENSO-teleconnected source regions, rather than by relying on a single well-correlated location, suggesting that ENSO variance estimates provided derived from a single site should be viewed with caution. Finally, identifying the common variance signal in a series of existing 20 proxy based reconstructions of ENSO variability over the


Introduction
The El Niño-Southern Oscillation (ENSO) is characterized by variations in sea surface temperature (SST) in the eastern tropical Pacific, causing changes in ocean currents and atmospheric circulation patterns globally. Related shifts in wind and rainfall patterns can lead to changes in extreme events including flooding, droughts, and trop-5 ical cyclone activity (Chan, 1985;Larkin and Harrison, 2002;Nicholls, 1985;Power et al., 1999). ENSO has been shown to exhibit significant multi-decadal variability in its strength and frequency throughout the instrumental period (Power, et al., 1999;Timmermann et al., 2003;Zhang et al., 1998), with additional longer-term variability reported in proxy data (Li et al., 2011;Wolff et al., 2011;Emile-Geay et al., 2013b). 10 Characterizing ENSO's long-term changes in frequency, magnitude and duration has been hampered by the fact that reliable instrumental records cover a period of less than 150 yr. This period is too brief to capture the range of long-term changes in ENSO frequency, magnitude and duration (Wittenberg 2009). Multi-century paleo climate reconstructions derived from monthly to annually resolved tree rings, ice cores, 15 lake sediments and coral records can be used to extend the observational record and to further quantify ENSO's sensitivity to external radiative perturbations in the "backdrop" of its internally generated variability (Federov and Philander 2000). Numerous attempts have been made to reconstruct ENSO variability back in time using various paleo-proxy archives (e.g., see Table 1 (2006). CM2.1 has a robust ENSO which includes realistic multidecadal fluctuations in amplitude, an irregular period between 2-5 yr, and SST anomalies that are skewed toward warm events as observed. The evolution of ENSO subsurface temperatures is quite realistic, as are ENSOs teleconnections as represented by correlations with precipitation anomalies outside the tropical Pacific. CM2.1s tropical climate and 10 variability is extensively described in Wittenberg et al. (2006).

Testing the running variance implications of a good temporal correlation
The fact that there is such a wide spread between the running variances of the individual ENSO reconstructions ( Fig. 1) raises the question whether a good temporal 15 correspondence between a given regional climate variable and ENSO can be used to imply that the variable will also provide a good representation of ENSO variance. We use output from the CCSM4 and CM2.1 models to address this question. Calculating the correlation between the simulated CM2.1 running variance of T S with the running variance of the modeled ENSO (Fig. 2a, shading), represented by the Nino 3.4 area-20 averaged (5 • S-5 • N, 120 • W-170 • W) sea surface temperature anomalies (SSTA), and comparing it with the absolute value of the correlation between anomalies of T S and modeled ENSO (Fig. 2a, contours) reveals a strong similarity between the two maps of correlation coefficients. Note, here and elsewhere in this manuscript that the running variance is calculated in a 30 yr sliding window. This similarity is reflected by the spa-25 tial correlation coefficient of 0.79 (r 2 = 0.63) calculated between the two spatial maps. Similar results are also obtained when the same calculations are performed on the CCSM4 simulation output (Fig. 2c) and again this is reflected by the spatial correlation of 0.70 (r 2 = 0.49) calculated between the two spatial maps. This implies that a large absolute value of the correlation between T S and ENSO provides a clear indication that the running variance of T S will track that of ENSO reasonably well.

5
Carrying out the same analysis for CM2.1 and CCSM4 precipitation data reveals some interesting differences between the running variance (precipitation running variance -ENSO running variance) correlation patterns and the raw (precipitation -ENSO) correlation patterns in the tropical Pacific ( Fig. 2b and d). These differences between the correlation maps are reflected by the spatial correlations (r) of 0.67 (r 2 = 0.45) and 0.55 (r 2 = 0.31) respectively. Figure 3 displays the 5th, 50th and 95th percentiles of the correlation coefficients calculated between the running variance of precipitation for all grid points and the running variance of ENSO binned according to the correlation between precipitation and ENSO. As expected from the significant differences between the two maps displayed in Fig. 2, 15 the 5th percentile curve reveals that having a strong correlation between precipitation and ENSO at a particular location does not guarantee a strong correlation between the running variance of precipitation and the running variance of ENSO. For instance, in CM2.1 if a precipitation signal is selected that has an r 2 value of between 0.6 and 0.7 when compared with the time series of ENSO SSTAs, there is a 10 % chance that 20 the running variance of that precipitation time series will have no significant correlation (r < 0.31 or r 2 < 0.1) with the running variance of ENSO (Fig. 3a). This indicates that ENSO may influence the sign and timing of the rainfall change at this location, however unrelated processes influence the magnitude of that change. We find that a strong correlation between the common precipitation signal, averaged 25 over two or more different locations and ENSO is a much better indicator for a high correlation between the running variance of the median (common) precipitation data and the running variance of ENSO (Fig. 3a)  precipitation signal has an r 2 of between 0.6 and 0.7 with ENSO, there is only a 1 % chance that that the running variance of that common precipitation time series will have no significant correlation (r < 0.31 or r 2 < 0.1) with the running variance of ENSO. This is 10 times less likely than the case with precipitation data only sourced from one location. This result is consistent with CCSM4 data which suggests a common precipitation 5 signal, from two geographic locations that have r 2 of > 0.7 when compared to ENSO SSTA, will make it 3.5 times less likely that the running variance of that signal will have no significant correlation (r < 0.3 or r 2 < 0.1) with the running variance of ENSO (Fig. 3b). We note that there are several issues with this type of analysis. The first is that 10 a large number of the proxies are sensitive to both temperature and rainfall (e.g., Cobb et al., 2003); the second issue is that all models have biases in the mean state, and ENSO's modeled teleconnections can be spatially shifted compared to observationspossibly implying biases in the spatial locations of the regions where the anomaly and running variance correlations with ENSO differ. We note, however, that both simulations 15 suggest that these regions (i.e., where rainfall is moderately correlated with ENSO, while the running variance displays little correlation with the running variance of ENSO) are in, or very close to, ENSO's centre of action in the western tropical Pacific. Regardless of these model issues, these results suggest: (i) a strong correlation between T S and ENSO provides a clear indication that the running variance of T S is 20 likely to track that of ENSO; (ii) a strong correlation between precipitation and ENSO at a particular location does not guarantee a strong correlation between the running variance of precipitation and the running variance of ENSO. As such, ENSO variance estimates provided by reconstructions derived from precipitation sensitive proxies at a single site should be viewed with caution; and (iii) a multi-site common precipitation 25 signal which is strongly correlated with ENSO provides much better indication that the common precipitation signal's running variance will be related to the running variance of ENSO.

Identifying the common variance signal
A previous study has shown that the linear combination of many of the above listed individual eastern equatorial Pacific reconstructions can further increase the signal-tonoise ratio compared with the individual reconstructions (McGregor et al., 2010). However, it was noted that the variance changes of this combined proxy product appeared 5 to show the effects of dating errors towards the beginning of the 350 yr record. Similar dating errors have also been noted in studies that attempt to reconstruct observed ENSO over the last 1000 yr (Emile-Geay et al., 2013a,b). This is a concern because proxies displaying small dating errors can act to cancel out the common signal when combined; hence artificially reducing the combined signal variance.

Methods for identifying the common variance signal
Here we will examine whether identifying the median of the running variance signals (MRV) of numerous source proxies (e.g., median (σ 2 (P x (t)))) is equivalent to identifying the running variance of the median signal (RVM) (e.g., σ 2 (median(P x (t)))) in the absence of dating errors. The source proxy timeseries are abbreviated as P x (t). This will 15 be accomplished by analyzing the output from the CCSM4 and CM2.1 models. To this end, we select x (x = 1, 2, 3, . . ., 9) time series of either precipitation or T S from random global locations, with the only constraint on the location being that its time series must have a absolute correlation coefficient > 0.3 when compared to modeled ENSO variability (again ENSO variability is calculated as average SSTA in the Niño 3.4 region).

20
We then calculate the MRV and RVM of the x time series. Note again that the running variance is calculated in a 30 yr sliding window. The above procedure is repeated until we obtain 10 000 MRVs and 10 000 RVMs for each x representing either precipitation or surface temperature. Plotting the 5, 50 and 95 percentiles of the correlation coefficients calculated be- the percentiles of the two derivation methods are very similar for both fields in both models (Fig. 4). This is also confirmed by plotting the correlation coefficients calculated between ENSO running variance and the RVM against the corresponding ENSO running variance and MRV correlation coefficients ( Fig. 5a and c). This indicates that while both methods will provide different estimates of the ENSO running variance time 5 series, MRV is roughly equivalent to RVM in the absence of dating errors. Furthermore, all percentiles increase as x increases, indicating that as number of proxies (x) increases the common signal is more likely to provide a skillful estimate of ENSO running variance. This confirms one of the initial assumptions used in many paleo climate studies that a network of proxies helps to reduce the effects of biases, 10 noise, and weaknesses in the individual indicators (Mann et al., 2000).

The effect of dating errors on the common running variance signal
Here we examine whether calculating the median individual proxy running variance time series (MRV), as opposed to calculating the running variance of the median of the proxy time series (RVM), reduces the effects of dating uncertainties in the reconstruc- 15 tions of ENSO variance changes as proposed by McGregor et al. (2010). To this end, five sets of precipitation time series are extracted from the CM2.1 simulation, which include all precipitation time series that have the absolute value of the correlation coefficients > 0.3 with modeled average SSTA in the Niño 3.4 region. In four of these five sets, the entire time series is shifted by between 1-5 yr to mimic a dating error. What 20 varies between these four sets, however, is the ratio of the time series that is subject to the introduced temporal shift, changing from 1/5, 1/4, 1/3, and 1/2.
Using these data and the approach outlined above in Sect. 3.2.1, we select x (x = 1, 2, 3, . . ., 9) precipitation time series from one of these data sets and then calculate the MRV and RVM. The 5, 50 and 95 percentiles of the correlations calculated between 25 ENSO running variance and the MRV for each of the five data sets as a function of x, reveal that the introduced temporal shifts have virtually no effect on the estimates of ENSO's running variance time series identified by the MRV method (Fig. 6a) the other hand, plotting the same percentiles for the correlations calculated between ENSO running variance and the RVM signal (e.g., σ 2 (median(P x (t)))) reveals that the introduced temporal shifts have a large effect on the estimates of running variance identified by the RVM method (Fig. 6b). This clearly indicates that the MRV method provides a much more robust estimate of 5 ENSO's running variance than the RVM method when dating errors are present in the input time series. We note that this result remains consistent regardless of the model (CM2.1 or CCSM4) or the output field used (i.e., precipitation or T S ).

Application to existing ENSO proxy reconstructions
Here we examine the MRV of the 14 existing ENSO reconstructions listed in Table 1 10 and presented in Fig. 1 1 . Each of these reconstructions is based on current and paleoproxy data sources from a range of geographic regions, with a variety of statistical methods used in the derivation of the relevant indices (see Table 1 and references therein). Considering that we have shown that the multi-proxy approach minimizes the impact of weaknesses in any individual proxy (Sect. 3.2.1), it is expected that the use of 15 this wide network of ENSO proxies will allow us to reconstruct the long-term changes in past ENSO variance with greater statistical skill than available from individual reconstructions alone.
Here we focus on the interannual ENSO frequency by high-pass filtering (10 yr cutoff period) each of the 14 ENSO reconstructions. Each of these filtered time series are 20 then normalized using the base period 1900-1977 for the calculation of standard deviation and a time series of thirty-year-sliding-window running variance is calculated for each reconstruction. The running variance time series are then adjusted by adding a constant, such that their mean running variance over the period 1900-1977 matches CPD 9,2013 Inferred changes in El Niño-Southern Oscillation variance that of the high-pass filtered observed N34 index over the same period. As shown in Fig. 1, there is a large spread between the running variances that highlight the discrepancies amongst the ENSO reconstruction variances (Fig. 7). Comparing the MRV of the ENSO proxy reconstructions with the 30 yr running variance of the high pass filtered (10 yr cutoff) ENSO signal derived from direct observa-5 tions during the instrumental period reveals a remarkable correspondence (Fig. 7). This further increases the confidence in the ensemble averaging technique applied here to the running variance. Figure 7 also reveals that the observed amplitude in the most recent 30 yr of the 20th Century , indicated by the red star) appears to have been larger than at any time during the preceding multi-century MRV record.
We note that the number of ENSO time series taken into account when calculating the MRV signal varies with time, depending on the length and span of the individual records ( Fig. 7, Table 1). Error bars are calculated using two separate pseudo-proxy methods, such that the length of the error bar varies in relation to the number of available ENSO proxies (see Appendix A). For instance, with fewer proxies available (e.g. 15 during the period 1400-1500) the estimators of the running variance become more uncertain than for a period (e.g. the 1800s) where many more proxy time-series overlap. At any point in time the most conservative (widest) error bar of the two pseudo proxy methods is displayed in Fig. 7.
The increased variance of ENSO for the most recent 30 yr is significantly larger (ex-20 ceeding the 95 % confidence level) than at any time during the past ∼ 400 yr (i.e., the observed variance, indicated by the red star, is larger than the median proxy signal error bars in the period from the year 1590 through to the start of the instrumental data). Prior to this period, the uncertainty range of the MRV signal, estimated by the error bars, is too large to be conclusive.

Independence of existing ENSO proxy reconstructions
We note that some of the input proxy networks overlap among these 14 ENSO reconstructions (see Appendix B for a more in depth discussion), such that they are not 2939 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | always completely independent in a geographical sense. In spite of geographical overlap for some of the reconstructions, usually different statistical methods are applied to tease out the ENSO signal. These methods range from simply calculating an Empirical Orthogonal Function (EOF) of the input proxy network (e.g., Braganza et al., 2009;McGregor et al., 2010;Li et al., 2011) to the more complex method of regressing proxy 5 data onto the leading spatio-temporal eigenmodes of observed data (e.g., Mann et al., 2000;Evans et al., 2002;Cook 2008). Despite the methodological differences, this apparent lack of independence raises the concern that the noise, i.e. the non-ENSO component of these overlapping input proxies, may distort or mask the common ENSO signal and this could be incorrectly 10 attributed to changes in ENSO. If this were true, we would expect the noise influencing tropical coral-dependent reconstructions to be distinct from that of North American tree rings. As such, if the running variance of either or both of these groups' reconstructions were dominated by noise we would expect to see significant differences in their common variance signal. However, this is not the case. Identifying the com-15 mon variance signal in only those reconstructions with common tropical corals reveals a post-1700 signal that is extremely similar to that identified in only those reconstructions with common North American tree rings (Fig. 8), thus lending further support of the RMV method applied here. In addition, both of the above signals resemble the common variance signal identified in the full set of ENSO reconstructions. This gives 20 us confidence that the common variance signal identified in this study is not strongly distorted or dominated by independent noise. We note that reconstructions 3 and 9 include common tropical corals and North American tree rings; however, neglecting these two reconstructions from both of the above subsets leaves the common variance signals virtually unchanged (figure not shown). Introduction

Application to observed single station proxies
It is also noted here that by using existing ENSO reconstructions that exhibit high correlations with the N34 observations during the instrumental period, we have attempted to consolidate the common information contained within these reconstructions rather than defining another reconstruction to add to the literature. However, as the recon-5 structions in most cases have already incorporated information from numerous time series, artificial variance reduction due to dating errors in the source proxies may already be incorporated into our analysis. As such, we may not be getting the best out of the methodological advantage of combining running variances when applying it to these synthesized time series.

10
Thus, here we also examine the common running variance signal in a range of single location proxies and compare the results with those derived above (Table 3). Again we focus on the interannual ENSO frequency by high-pass filtering (10 yr cutoff period) each of the proxies. We utilize all tree ring proxies listed by Emile-Geay et al. (2013a, b) and all tropical Pacific coral proxies available from the NOAA National Climate Data 15 Centre (NCDC) that continuously span the period 1800-1980. However, prior to their inclusion in Table 3, these proxies were pre-screened to ensure they have a correlation of at least 0.3 when compared with the annual mean (July-June) observations of ENSO from any of the four SST products used in Table 4. Each of these filtered time series are then normalized over the 1900-1977 period and a time series of running variances, 20 in a thirty year sliding window, is calculated for each reconstruction. The variance time series are then adjusted by adding a constant, such that their mean running variance over the period 1900-1977 matches that of the high-pass filtered observed N34 index over the same period.
Calculating the common variance signal of all proxies listed in Table 3 (Table 1). In fact, comparing this common running variance time series with that of the 14 existing ENSO reconstructions reveals a strong similarity through most of the 600 yr period (Fig. 9a). Using only the 12 tree ring proxies (sourced from Emile-Geay et al., 2013; Table 3) we calculate the MRV signal. Comparing this running variance signal with that of the 5 observed high pass filtered ENSO signal reveals a good correspondence (Fig. 8b inset). Again, this common variance time series falls well within the error estimates of the 14 existing ENSO reconstructions, and directly comparing this common running variance time series with that of the 14 existing ENSO reconstructions, reveals a strong similarity through most of the 600 yr period (Fig. 8b).

10
Calculating the MRV signal of the 9 coral proxies from Table 3, and comparing this running variance signal with that of the observed HPF ENSO signal reveals some correspondence ( Fig. 8c inset), although the correlation appears much weaker than for the individual tree-ring records (Fig. 8b inset). The correspondence is most notable prior to 1960; after this time the observed RV begins to steadily increase while the coral MRV 15 signal remains relatively constant. The facts that (i) the coral MRV signal does not pick up the increase in ENSO variance since the mid-1900s; and (ii) there is very little correspondence between this signal and the ENSO proxy reconstruction MRV time series outside of the period 1850-1950 (Fig. 8c); reduces our confidence in this signal, which is namely sourced from the southwestern tropical Pacific, as it relates to ENSO. As the 20 southwestern tropical Pacific is home to the South Pacific Convergence Zone (SPCZ), a possible explanation for these discrepancies could be that there may have been slow shifts in the tropical Pacific temperature/rainfall covariance that could have affected some of the δ 18 O-based coral proxies in this region. However, in spite of this, the coral common variance time series falls mostly within the error estimates of the ENSO proxy 25 reconstruction MRV signal. Thus, consistent with the analysis of existing ENSO reconstructions presented in Fig. 7, the common running variance signal of the observed amplitude in the most recent 30 yr of the 20th Century (indicated by the red star) appears to have been larger than at any time during the preceding multi-century record.

Summary and conclusions
The main goal of this study was to synthesize existing ENSO reconstructions to arrive 5 at a better estimate of past ENSO variance changes. This is an important task since there are considerable discrepancies in ENSO variance estimates derived from the 14 available ENSO reconstructions that exhibit good correlations with the instrumental data ( Figs. 1 and 7). We first tested an implicit assumption of many ENSO reconstruction studies, namely, 10 that if a proxy (whether it represents precipitation or surface temperature) represents ENSO variability well, it will also represent ENSO running variance well. Our results from the analysis of two independent multi-century CGCM simulations suggest that this assumption is valid for proxies of temperature variability, but not for single proxies of precipitation variability. That is, while local precipitation may be modulated by 15 ENSO, its running variance might bear a much weaker correlation with ENSO's running variance. Thus, we recommend caution when inferring ENSO running variance from precipitation-sensitive proxy data from a geographically confined area, even if they are located in the central/eastern equatorial Pacific. Furthermore, we demonstrated that this effect is significantly reduced when a common precipitation signal is sourced from 20 a range of geographic regions. That is, a common precipitation signal, which is strongly correlated with ENSO, provides a much better indication that the correlation between the rainfall running variance and ENSO running variance will also be strong. We then showed that in the absence of dating errors, the median running variance signal (MRV) is roughly equivalent to calculating the running variance of the median Introduction much more sensitive to dating errors than the former. Importantly, these dating errors can lead to an artificial variance reduction due to the cancellation of the ENSO signal.
The results of our analysis also show that as the diversity of input proxies increases, the resulting common variance signal is more likely to provide a skillful estimate of ENSO running variance. This confirms one of the initial assumptions used in many 5 paleo climate studies that a broad network of proxies helps to reduce the effects of biases and weaknesses in the individual indicators (Mann et al., 2000). Thus, we expect the MRV method to increase the signal to noise ratio of the ENSO variance changes, while minimizing the effects of dating uncertainties in the source proxies.
There are two important limitations of this study. Firstly, although our combined vari-10 ance estimate is supported by the running variance signal calculated from uncalibrated single station corals (see Sect. 5), we do not assess the sensitivity of individual source ENSO reconstructions to their derivation and calibration techniques. A recent study has suggested that reconstructed centennial scale variability is particularly sensitive to the historical SST data set used during calibration (Emile-Gaey et al., 2013a,b). Secondly, 15 we assume that the spatial structure of ENSO SSTA and ENSO's teleconnections do not change over time. This is a concern, as the single station proxies (Table 3) are largely sourced from regions considered teleconnected to ENSO and do not include any proxies from ENSO's centre of action in the eastern/central tropical Pacific. Thus, while the common running variance signals of these single station proxies support the 20 synthesis of predefined ENSO reconstructions, the differences between the common running variance time series could reflect changes of teleconnection patterns over time ( Fig. 9b and c). With these potential limitations in mind, the synthesis of the fourteen predefined ENSO reconstructions (Table 1), along with the analysis of 21 single station proxy 25 records ( Within the current methodology, sources of uncertainty could be reduced by (i) increasing the number of ENSO proxy reconstructions representing ENSO variability prior to 1600 CE; and/or (ii) increasing the number of annually resolved single station 5 proxies that represent ENSO variability. The increase in single station proxies should preferably occur outside of the two regions most commonly used in ENSO reconstructions (Table 3), namely central/north America (tree rings) and the southwest tropical Pacific (corals).

Error bar calculation
The main goal of this pseudo-proxy technique is to assess how accurately the median of multiple pseudo-proxies can estimate the variance of the original proxies from which the pseudo-proxies are generated. Here we define the original source proxy as the "true" signal. To this end, each of the 14 normalized high-pass filtered proxy timeseries 15 p i (t)(i = 1, . . ., 14) (Table 1) is used as a "true" signal. Pseudo-proxies pp i j (t) = p i (t) + sz j (t)(j = 1, . . .J; J = 600) are generated by adding Gausian white noise z j (t) to the annual mean "true signal". The variance of z j (t) is randomly modulated by s, such that the correlation coefficient between the pseudo-proxy pp i j (t) and the original proxy p i (t) (the "true" signal) is contained within the range of correlation coefficients between the 20 original proxies and the instrumental SSTA signal (see Table 2). By definition, selecting numerous pseudo-proxies and averaging (mean(pp i (t)) = J −1 S j pp i j (t)) should allow for the "true" signal p i (t) to be identified, as the additive Gaussian white noise would cancel out if enough pseudo-proxies were utilized. If none of the proxy records had dating errors, the 30 yr running variance of mean(pp i (t)) would match that of the "true Introduction averaging cannot completely attenuate the noise. Furthermore, the available proxy data will likely have dating errors leading to an artificial damping of the combined signal variance.
To bypass this signal cancellation issue, the 30 yr running variance of pp i j (t) is calculated, giving pp rv i j (t). The variance of these pp rv i j (t) time series are then adjusted by 5 adding a constant, such that their mean running variance over the period 1900-1977 matches that of the "true proxy" over the same period. From this point two methods can be used to estimate the error bar for the running variance time series: (1) the first method generates n = 5000 groups for each proxy (i ), that contains m = 3 randomly selected pseudo proxy running variance time series (pp rv i j (t)); and 5000 groups that 10 contain m = 4 randomly selected pseudo proxy running variance time series (pp rv i j (t)); etc with m going up to 14). Hence the total number of realizations of pseudo-proxy running variances in this ensemble is 14 × 102 × 5000, where 102 = 14 m=3 m. The median of this set is then calculated with respect to m; e.g. for m = 3, 4,. . . , 14, and for any given i we obtain the median over these m realizations in each of the 5000 groups 15 => mpp rv i mn (t) giving a 14 × 12 × 5000 matrix. Then for each i we calculate the error of these 5000 groups of median running variances with respect to the corresponding true running variance giving ε rv i mn (t) = p rv i t − mpp rv i mn (t). The matrix ε rv i mn (t) is also of dimension 14 × 12 x 5000. To obtain an estimate of the uncertainty in running variances that can be applied to the median ENSO proxy variance presented here, and to main-20 tain positive definiteness of the variance range, the following procedure is applied: the true proxy running variance as a function of time p rv i (t) is split into p = 8 bins of equal width (going from 0.25 to 2) depending on the value of p rv i (t); also the n corresponding error estimates (ε rv i mn (t)) are binned into the same 8 equal bins => binned-ε rv i mp (q) is a 14×12×8×p matrix where the number of error estimates in each bin (q) is a multiple 25 of n. For each bin and each proxy i the 5 and 95-percentiles of the realizations is computed within this bin. For each i and m there are now 8 sets of error bars, which depend on the magnitude of the corresponding true running variance signal. A conservative estimate of the error bars is obtained by taking the minimum or maximum of these 5  = 3, . . ., 14), the corresponding proxy-variance magnitudedependent conservative error bars are selected for the corresponding m to get the final proxy variance uncertainty range.
(2) The second method simply calculates the median of the pp rv i j (t) over the i di-5 mension, giving mpp rv j (t). As the number of original proxy time series at any point in time is consistent with the magenta line in Fig. 3, this calculation directly provides 600 pseudo-estimates of the proxy median time series displayed in Fig. 2a (black line). The error bars calculated via this method implicitly incorporate the number of proxies used to estimate the variance at each time step and the magnitude of the signal it is trying to 10 represent, so error bars are then calculated at each point in time by simply identifying the 5 % and 95 % levels.
As the addition of white noise to the pseudo-proxies acts to increase the mean variance of each pseudo-proxy time series, the shifting of each of the pseudo-proxy running variance time series by adding a constant (such that the mean of the last 15 100 yr matches that of the "true" signal p rv i (t)) may act to shift the positive definite pseudo-proxy error bars towards more negative values; hence making it possible that the pseudo-proxy error bars include negative values. Such values will not be considered here, because of the positive definiteness of variances.

Proxy data independence
In Sect. 4, we noted that some of the input proxy networks used to develop these 14 ENSO reconstructions overlap. Neglecting reconstruction 10, which utilizes reconstructions 1-9 in its derivation along with the Urvina Bay corals (Dunbar et al., 1994), this input proxy data overlap can be broken down into two main groups, those utilizing: Introduction  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Landrum, L., Otto-Bliesner, B. L., Wahl, E. R., Conley, A., Lawrence, P. J., Rosenbloom, N., and Teng, H.: Last Millennium Climate and Its Variability in CCSM4, J. Climate, 26, 1085-1111, doi:10.1175/JCLI-D-11-00326.1, 2013 and Cold (La Niña) event life cycles: ocean surface anomaly patterns, their symmetries, asymmetries, and implications, J. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Quinn, T. M., Crowley, T. J., Taylor, F. W., Henin, C., Joannot, P., and Join, Y.: A multicentury stable isotope record from a New Caledonia coral: interannual and decadal sea surface temperature variability in the southwest Pacific since 1657 AD, Paleoceanography, 13, 412-426, 1998. Rayner, N., Parker, D., Horton, E., Folland, C., Alexander, L., Rowell, D., Kent, E., 5 and Kaplan, A.: Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century, J. Geophys. Res., 108, 4407, doi:10.1029/2002JD002670, 2003  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Table 3. The single station (SS) temperature and/or rainfall proxies employed in this study. Note that these proxies have at least annual resolution and are pre-screened such that they have a correlation of at least 0.3 with the annual mean (July-June) observations of ENSO, represented by Niño 3.4 area-averaged (hereafter N34, 5 • S-5 • N, 120 • W-170 • W) SST anomalies (SSTA; Table 4). This correlation cutoff ensures that the proxies provide a reasonably skillful estimate of ENSO variability over the instrumental epoch.    Table 1  ENSO. The color of the lines indicates the number of source locations utilized in the 5 precipitation signal, with black, blue and red lines corresponding to one, two and five source 6 locations. The inset histogram displays the distribution of squared correlation coefficients 7 calculated between precipitation running variance and the running variance of ENSO from the 8 identified x-axis (grey shading) bin. b) as in a) but for CCSM4 model output. Note the different 9 Y-axis for panels a) and b). running variance (N34RV) are displayed in red, while those of the RVM surface temperature are 5 displayed in black. b) as in a) but for the CM2.1 precipitation data. c) and d) as in a) and b) but 6 for CCSM4 surface temperature and precipitation data respectively. The x-axis in each panel 7 displays the number of grid point time series used to calculate the common signal.  the correlation coefficients calculated between the CM2.1 rainfall MRV and ENSO running 4 variance, while b) displays the percentiles of correlation coefficients calculated between CM2.1 5 rainfall RVM and ENSO running variance. The black lines indicate those percentiles using data 6 with no introduced temporal shifts, while the red, yellow, green and blue line respectively 7 represent those percentiles using data with 1/5, 1/4, 1/3, 1/2 of the time series including an 8 introduced temporal shift. plus x's) of the correlation coefficients calculated between the CM2.1 rainfall MRV and ENSO running variance, while (b) displays the percentiles of correlation coefficients calculated between CM2.1 rainfall RVM and ENSO running variance. The black lines indicate those percentiles using data with no introduced temporal shifts, while the red, yellow, green and blue line respectively represent those percentiles using data with 1/5, 1/4, 1/3, 1/2 of the time series including an introduced temporal shift. 9,2013 Inferred changes in El Niño-Southern Oscillation variance reconstructions with common tree ring source proxies (Proxy numbers correspond to those in 4 Table 1), overlaid with the corresponding ensemble median running variance (thick dashed red 5 line). b) the 30-yr running variance of each of the 6 HPF (10-year cutoff) ENSO reconstructions 6 with common tropical coral reconstructions (Proxy numbers in legend correspond to those in 7 Table 1), overlaid with the corresponding ensemble median running variance (thick dashed red 8 line). The underlying thick black line in both panels represents the median running variance of 9 all 14 HPF (10-year cutoff) ENSO reconstructions presented in Fig. 6, while the thin black lines  Table 1), overlaid with the corresponding ensemble median running variance (thick dashed red line). (b) The 30 yr running variance of each of the 6 HPF (10 yr cutoff) ENSO reconstructions with common tropical coral reconstructions (proxy numbers in legend correspond to those in Table 1), overlaid with the corresponding ensemble median running variance (thick dashed red line). The underlying thick black line in both panels represents the median running variance of all 14 HPF (10 yr cutoff) ENSO reconstructions presented in Fig. 7  The running variance of the 12 single station tree ring proxies (light green dots) with the median (thick green dashed line) overlaid, while (c) displays the running variance of the 9 single station coral proxies (light purple dots) with the median (thick purple dashed line) overlaid. The thick black line in both panels represents the median running variance of all 14 HPF (10 yr cutoff) ENSO reconstructions presented in Fig. 7, while the thin black lines represent the accompanying running variance signal error bars. The blue line represents the observed ensemble median running variance and the red star indicates its most recent value, while the thin red line just extends this most recent value back through time for comparison with the ensemble median proxy running variance.