Was the Little Ice Age more or less El Niño-like than the Medieval Climate Anomaly ? Evidence from hydrological and temperature proxy data

The El Niño–Southern Oscillation (ENSO) is the most important source of global climate variability on interannual timescales and has substantial environmental and socio-economic consequences. However, it is unclear how it interacts with large-scale climate states over longer (decadal to centennial) timescales. The instrumental ENSO record is too short for analysing long-term trends and variability and climate models are unable to accurately simulate past ENSO states. Proxy data are used to extend the record, but different proxy sources have produced dissimilar reconstructions of long-term ENSO-like climate change, with some evidence for a temperature–precipitation divergence in ENSO-like climate over the past millennium, in particular during the Medieval Climate Anomaly (MCA; AD ∼ 800–1300) and the Little Ice Age (LIA; AD ∼ 1400–1850). This throws into question the stability of the modern ENSO system and its links to the global climate, which has implications for future projections. Here we use a new statistical approach using weighting based on empirical orthogonal function (EOF) to create two new large-scale reconstructions of ENSOlike climate change derived independently from precipitation proxies and temperature proxies. The method is developed and validated using model-derived pseudo-proxy experiments that address the effects of proxy dating error, resolution, and noise to improve uncertainty estimations. We find no evidence that temperature and precipitation disagree over the ENSO-like state over the past millennium, but neither do they agree strongly. There is no statistically significant difference between the MCA and the LIA in either reconstruction. However, the temperature reconstruction suffers from a lack of high-quality proxy records located in ENSO-sensitive regions, which limits its ability to capture the large-scale ENSO signal. Further expansion of the palaeo-database and improvements to instrumental, satellite, and model representations of ENSO are needed to fully resolve the discrepancies found among proxy records and establish the long-term stability of this important mode of climatic variability.


Introduction
The El Niño-Southern Oscillation (ENSO) is the most influential source of interannual variability in the modern climate.The warm El Niño state is characterised by a weaker sea surface temperature (SST) gradient across the equatorial Pacific and a shift in precipitation from the western Pacific toward the central Pacific, while the cool La Niña state is roughly the opposite.Although ENSO originates in the tropical Pacific, it has far-reaching effects through teleconnections on some regions in higher latitudes, and El Niño years are generally anomalously warm on a global scale.However, it is unclear whether there is a link between anomalously warm or cool periods and the two ENSO states on decadal to centennial timescales.Given the severe socio-economic consequences of ENSO events (Hjelle and Glass, 2000;Page et al., 2002;Kovats et al., 2003;Badjeck et al., 2010), and a warmer future under continued anthropogenic warming, it is important Published by Copernicus Publications on behalf of the European Geosciences Union.L. M. K. Henke et al.: ENSO reconstructions to understand the natural long-term ENSO and its interaction with the climate.It allows for an evaluation of the effects of anthropogenic impacts on recent and future ENSO behaviour (Collins, 2005;Guilyardi et al., 2009;Vecchi and Wittenberg, 2010;Bellenger et al., 2014).
Recent multi-model studies of projected changes in ENSO under anthropogenic warming suggest robust changes to ENSO-driven temperature and precipitation, including an increase in extreme El Niño (Cai et al., 2015a) and La Niña (Cai et al., 2015b) events and changes in the ENSO SST pattern and ENSO-driven precipitation variability (Power et al., 2013).However, most current general circulation models (GCMs) cannot simulate many aspects of the modern-day ENSO accurately, often overestimating the western extent of the Pacific cold tongue and failing to correctly simulate central Pacific precipitation anomalies, ENSO feedbacks, and ENSO amplitude (Bellenger et al., 2014;Collins, 2005;van Oldenborgh et al., 2005).This translates into uncertainty over simulations of past ENSO-like climate change, calling for alternative sources of climatic information to supplement, complement, and corroborate the model and instrumental data.This is done using proxy climate records such as tree ring widths, tropical ice cores, sediment cores, and corals (Jones and Mann, 2004).
There are very few annually resolved proxy records available longer than ∼ 500 years (Mann et al., 2008).The issue with high-resolution proxies is that they tend to be short in length; trees and corals, for example, rarely live beyond a few centuries (Jones and Mann, 2004).Some such highly resolved records are available for the more distant past, but these generally offer snapshots rather than continuous records (McCulloch et al., 1996;Corrège et al., 2000;Abram et al., 2009).However, there are several long, lowerresolution proxy records of ENSO variability on decadal to millennial scales, often derived from lake sediments (see Conroy et al., 2010), marine sediments (see Barron and Anderson, 2011), or speleothems (see Maupin et al., 2014).While these are unable to capture the interannual frequency and amplitude of individual ENSO episodes, they provide an insight into longer-term ENSO-like climate states and average ENSO behaviour.Although there are some endeavours to combine some low-resolution proxies, often to capture spatial gradients (Conroy et al., 2010;Yan et al., 2011b;Anderson, 2012), there has not, to our knowledge, been a comprehensive effort to systematically merge a large set of such lowresolution records to create a long-term reconstruction of ENSO-like climate variability.Doing this could shed light on the long-term stability of ENSO and its links with the wider climate, for example by examining ENSO behaviour under different dominant cool or warm climate states, which in turn can inform our understanding of potential future ENSO-like changes in a warmer world.
A number of proxy, instrumental, and modelling studies investigate links between ENSO and global climate variability on interannual (Klein et al., 1999;Wang et al., 1999), decadal (Nelson et al., 2011), centennial (Mann et al., 2005;Trouet et al., 2012), and millennial (Cane, 2005;Ivanochko et al., 2005;Moy et al., 2002;Shin et al., 2006) timescales.A wide range of proxy records and modelling studies point to substantial shifts in the ENSO-like state of the climate linked to changes in solar variability on orbital timescales (Clement et al., 2001;Barron and Anderson, 2011), movement of the Intertropical Convergence Zone (ITCZ; Partin et al., 2007;Gomez et al., 2004;Carré et al., 2005;Nelson et al., 2011), and changes in ocean circulation linked to sea level rise (Wanner et al., 2008;McGregor et al., 2008).In the more recent past, changes in solar irradiance and stratospheric aerosol loadings due to volcanic activity have played significant roles in modulating the hemispheric to global scale climate (Mann et al., 2009).The so-called Medieval Climate Anomaly (MCA) and the Little Ice Age (LIA), generally defined to fall within ca.AD 900-1300 and AD 1300-1850, respectively (Jones and Mann, 2004), are two periods of climate upheavals based on Northern Hemisphere (NH) climate variability.The MCA was a period of anomalously warm conditions, while the LIA was anomalously cool, at least in the Northern Hemisphere (Jones and Mann, 2004;Mann et al., 2009).These two periods are often used for exploring past behaviour of climatic phenomena as they represent relatively large and sustained excursions from the long-term mean.However, a comparison of hemispheric temperature reconstruction (Neukom et al., 2014) and continental-scale temperature reconstructions (PAGES2k Consortium, 2013) finds no evidence of a globally coherent MCA and only partial evidence for a global LIA.Neither of these two studies focus specifically on the equatorial tropics, leaving open the question of the strength of a potential MCA and/or LIA in these latitudes.This is an important knowledge gap because proxy evidence for the expression of ENSO-like climate change over these periods appears to be ambiguous.A range of proxies point to a more northerly ITCZ (Haug et al., 2001;Sachs et al., 2009;Tierney et al., 2010) during the MCA, which is characteristic of La Niña-like conditions and is in agreement with warming patterns found in multi-proxy reconstructions of hemispheric and global-scale temperature (see Mann et al., 2009).Langton et al. (2008) similarly infer a reduction in El Niño-like activity during the MCA based on ocean basin ventilation changes in Indonesia.In contrast, a Southern Oscillation Index reconstruction based on two proxy records (Yan et al., 2011b) shows an El Niño-like state during the MCA (defined as AD 800-1200) and a La Niña-like LIA (AD 1300(AD -1850)).This seems to be supported by a number of other precipitation proxies from the West Pacific (Yan et al., 2011a) and East Pacific (Moy et al., 2002).Other precipitation proxies indicate a highly variable ENSO during the LIA, including two multi-decadal droughts in Java (Crausbay et al., 2006), high-amplitude rainfall fluctuations in Madagascar (Brook et al., 1999), and three southerly ITCZ excursions (Haug et al., 2001).
The ENSO-like state of the climate may be linked to the Interdecadal Pacific Oscillation (IPO).The IPO is characterised by a Pacific SST anomaly pattern resembling ENSO (Power and Colman, 2006), which oscillates on decadal timescales.The North Pacific section (20-45 • N) is often referred to as the Pacific Decadal Oscillation (PDO).The interactions between the PDO/IPO and ENSO are still not well understood beyond a statistical relationship (Wang et al., 2014b), but there is evidence for interactions between ENSO behaviour and PDO state (Yuan Zhang et al., 1997) on multidecadal timescales over the instrumental period (Rodgers et al., 2004;Schopf and Burgman, 2006;Newman et al., 2003) and further back in time (Verdon and Franks, 2006;McGregor et al., 2010).Some studies suggest that these decadal oscillations are essentially integrated long-term expressions of ENSO (Newman et al., 2003) and that they can be explained by stochastic ENSO fluctuations on decadal timescales (Power and Colman, 2006).A phase change analysis of palaeoclimate data over the past 400 years (Verdon and Franks, 2006) finds a tendency for more frequent El Niño events during positive PDO phases and more frequent La Niña events during negative PDO phases.Over Australia, the IPO phase modulates the strength of ENSO influence on rainfall extremes, with a more pronounced effect of La Niña during negative IPO (Power et al., 1998).The spatiotemporal robustness of the IPO is questionable as South Pacific coral records (Linsley et al., 2004) find low spatially coherent IPO related prior to the mid-1880s, suggesting that the spatial pattern may have changed.Additionally, documentary rainfall data from China (Shen et al., 2006) indicate that while the PDO has existed for at least ∼ 530 years, its periodicity changed from 75-115 years pre-1850 to 50-70 years post-1850.
The discrepancies in long-term ENSO-like variations between proxy records raise two important questions.The current dynamical understanding of ENSO is underpinned by the strong relationship between temperature and rainfall observed today and the relationships between the ENSO source region in the tropical Pacific and teleconnected regions, which largely fall between 40 • S and 40 • N. As Yan et al. (2011b) highlight, however, temperature and precipitation proxies appear to disagree on the ENSO-like states of the MCA and the LIA.Thus, to what extent does the modernday precipitation-temperature relationship in the source and teleconnected regions continue to exist in the past?The second question concerns the relation of ENSO to the wider climate; is there a link between global temperatures and longterm ENSO state on multi-decadal to centennial timescales?A comparison between the MCA and the LIA can give some insight into this and may hold some clues to what we can expect under anthropogenic climate change.
The use of proxy archives can contribute valuable insights on past climate variability by extending the instrumental records back in time, but substantial uncertainties remain.This is because all reconstructions have inherent limitations and ambiguities that must be identified and dealt with appropriately.These include resolution, dating errors, noise, limited and/or skewed spatial coverage, and nonlinear responses to the climatic variable of interest (Jones et al., 2009).Various statistical techniques have been employed to create multi-proxy reconstructions of climatological phenomena, broadly falling into the categories "composite plus scaling" (CPS) or "climate field reconstruction" (CFR) (Jones et al., 2009).CPS encompasses any method that involves combining standardised proxy records into a single reconstruction, which is subsequently calibrated to a known time series of the target variable (e.g.instrumental temperature record) to provide a quantitative estimate of the variable.CFRs, on the other hand, aim to reconstruct large-scale spatial patterns of climatic change using covariance between proxies and instrumental data.Within both methods there is a wide variety of approaches; see Jones et al. (2009) for detailed descriptions and examples of both CPS and a range of CFR methods.The focus of this study -comparing the climate signals in temperature and precipitation proxy records separately -calls for a slightly different approach.
Here we create two new ENSO reconstructions, one derived from temperature proxies and one from precipitation proxies, using a new method for assessing the stability of the modern-day ENSO patterns in the source region and the wider teleconnected regions.In a fashion similar to Braganza et al. (2009), for example, proxy records are not tuned to instrumental data other than a simple location-dependent weighting.While this precludes direct quantitative comparisons, it removes the bias towards high-frequency trends that stems from calibrating to the relatively short (∼ 150 year) instrumental record (or indeed any short record; Cook et al., 1995;Jones and Mann, 2004).The method amplifies the ENSO component of proxy records and simultaneously attempts to quantify uncertainty related to noise and incomplete spatio-temporal data coverage, whilst maximising the use of a wide range of tropical proxies.With this, we aim to answer two questions: 1. Do temperature and precipitation proxies show consistent long-term ENSO behaviour over the last millennium?
2. Do the LIA and the MCA differ significantly in their mean ENSO state?
2 Data description

Proxy records
For this study, a comprehensive effort was made to collect all published proxy precipitation and temperature records between 40 • S and 40 • N that cover the last 2000 years.The large majority of records were accessed from the NOAA Paleoclimatology and Pangaea Databases (https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets).In addition, over 200 tree ring records were taken from the dataset used by Mann et al. (2008) and hence were subject to their criteria, including length, intra-site signal coherence, and sample density (see Appendix A for details).A set of coral records was taken from the dataset made available by Tierney et al. (2015); these were largely temperature proxies, although some were assigned as precipitation proxies or excluded altogether based on information in the original publications (see Appendix A for details).
After collection, all records were screened for a maximum dating error of 60 years.Although somewhat arbitrary, this cut-off was decided by taking double the averaging bin width of 30 years that was applied to the data prior to analysis (Appendix B3).This is a step towards addressing the issue of dating uncertainty while allowing a wider range of proxies to be utilised.Proxies with larger dating errors generally have lower (multi-annual to multi-decadal) resolution but are also usually much longer and are arguably more useful for capturing long-term trends that may be less evident or reliable in annual resolution proxies (cf.Cook et al., 1995;Esper et al., 2002;Mann et al., 2008, on low-frequency trends in tree rings).Other quality judgements regarding temporal resolution, record length, and proxy location are accounted for by the method set out in Sect. 3 and Appendix B.

Modern climate datasets
Instrumental climate data are the best available in terms of dating accuracy, calibration, and physical basis (Jones and Mann, 2004).However, their spatial coverage is not complete and sharply decreases back in time.The nature of the method used in this study calls for full spatial coverage over a long period; therefore, reanalysis products are more suitable.These are combinations of instrumental and satellite data interpolated using models.The 20th Century Reanalysis Version 2c (20CRv2c) is the longest global dataset of atmospheric circulation available, spanning AD 1851-2014.It is based on surface pressure, temperature, and sea ice distribution data, filled in with a deterministic ensemble Kalman filter (EKF).It has a spatial resolution of 2 • latitude × 2 • longitude × 24 vertical pressure levels and a temporal resolution of up to 6 h.It has been demonstrated that the 20CRv2c is competent at representing the global tropospheric circulation as well as the mean state and variability of the hydroclimate; for a detailed description and evaluation of the product see Compo et al. (2011).Comparison with ENSO indices indicates that 20CRv2c accurately represents the temporal evolution of ENSO (Table B1), and a recent comparison with other reanalyses (ERA-20C, JRA-55, and ERA-Interim) indicates that 20CRv2c shows similar variations in various climate indices, including NINO3.4 and SOI (Poli et al., 2016).The monthly mean surface air temperature and precipitation rate datasets were downloaded from the NOAA/ESRL/OAR-PSD website (http://www.esrl.noaa.gov/psd/).The 20CRv2c data were regridded to 2 • × 3 • to be comparable to the model data described below (Sect.2.3).The climatology (for the period 1851-2014) was removed to produce monthly anomalies, which were then averaged to annual resolution.
As noted in Sect. 1, the ENSO-like temperature pattern is similar to the IPO.Regardless of the directionality of influence between the two modes, the long-term ENSO-like state may thus be similar to the IPO/PDO phasing.The ENSO-like state (and the associated spatial pattern, here derived from 20CRv2c) is here referred to as ENSO-like.

General circulation model simulations
There are several comprehensive modelling projects with the aim of improving comparability between GCMs produced by different teams.GCMs taking part in these projects perform a set of simulations with standardised forcings and boundary conditions.For this study, the pre-industrial control (piControl; pre-1850 parameters, no external forcings) and historical (AD 1850(AD -2000) ) runs from the Coupled Model Intercomparison Project 5 (CMIP5; Taylor et al., 2012) were used, in addition to the last millennium (past1000; AD 850-1850) runs from the Paleoclimate Model Intercomparison Project 3 (PMIP3; Braconnot et al., 2011Braconnot et al., , 2012)).Of the six GCMs that have all three runs available, two were chosen for their similarity to 20CRv2c in terms of spatial ENSO representation (see Appendix B2) for precipitation (climate modelling groups in brackets): -CCSM4 (National Center for Atmospheric Research), -GISS-E2-R (NASA Goddard Institute for Space Studies).
All model datasets were regridded to 2 • × 3 • , converted to anomalies, and degraded to annual resolution to enable comparison between models and 20CRv2c.

Methodology
The method used in this study consists of two stages: first, a temperature and a precipitation ensemble of proxy networks is created based on GCM-derived pseudo-proxy experiments; then, these network ensembles are used to produce two independent reconstructions of ENSO-like climate change using the real proxy data weighted by spatial temperature and precipitation ENSO patterns from 20CRv2c reanalysis data.This approach attempts to take into account the effects of proxy selection, the temporal limitations of individual proxies, and non-climatic noise.A brief overview of the method is given here, with a more extensive description in Appendix B.

Pseudo-proxy creation
Pseudo-proxies are simulated proxies that attempt to mimic various sources of uncertainty inherent in the real proxy records.This ranges from adding white Gaussian noise with a prescribed signal-to-noise (SNR) ratio approximating nonclimatic random noise to more sophisticated process-based additions that take into account effects such as dating error, nonlinear and multivariate responses of the proxy sensor to the climate variable, and sampling biases (Mann, 2002;Smerdon, 2012).The utility of any pseudo-proxy exercise lies in the fact that the answer to the question is known as it can be derived directly from the original model dataset.By putting the "signal plus noise" pseudo-proxies through a method to make a reconstruction of the signal, it allows for inferences about the stability and limitations of the method, estimates of uncertainties due to noise and other proxy record characteristics, and, as highlighted by the method used here, it provides a way of objectively and systematically selecting the most appropriate data (see Smerdon, 2012, for an introduction to pseudoproxy experiments).
Ideally, all pseudo-proxies in this study would be created from proxy system models (PSMs; e.g.Dee et al., 2015).Unfortunately, many proxies (e.g.δ 18 O from corals or speleothems) would require isotope-enabled GCMs, which are not available for the model runs used in this study.Tree ring widths (TRWs), however, can easily be simulated using the VS-Lite Model, which is a freely available PSM designed to estimate TRW based on minimal climatic input (Tolwinski-Ward et al., 2011;Breitenmoser et al., 2014;Evans et al., 2014).It is a simplified version of the full Vaganov-Shashkin model (Vaganov et al., 2006;Anchukaitis et al., 2006;Vaganov et al., 2011) and only requires temperature, precipitation, and latitude information.All TRW records in the proxy dataset used in this study were thus represented using VS-Lite-derived pseudo-TRW series.For the rest of the proxies, the GCM raw temperature or precipitation series were used.
Coral δ 18 O is known to be particularly vulnerable to showing multivariate responses, depending on both the δ 18 O of the seawater (which can be altered through changes in sea surface salinity (SSS) linked to ocean circulation or precipitation amount) and SST, which is usually the desired variable (Evans et al., 1998;Russon et al., 2013).Since there is no simple method of determining the extent of this contamination of the SST signal without an isotope-enabled GCM and coral PSM, the decision has been made here to exclude all coral δ 18 O records.Other coral proxies such as Mg / Ca or Sr / Ca have been retained as they are more directly associated with SST variations.Exclusion of these coral records equates to the loss of 12 potential proxies that meet the minimum length criterion.Testing showed that the exclusion of these proxies did not alter the conclusions of this paper (not shown).
White Gaussian noise was added to all model temperature, precipitation, and TRW series.Where the original publications (or in the case of the coral records, Tierney et al., 2015) provided an indication of the strength of the proxy-climate relationship (e.g.R value), the SNR was calculated.For all other records, a prescribed SNR of 0.4 (which corresponds to an R value of ∼ ± 0.38) was used as this has been shown to be a realistic average (Smerdon, 2012).The median calculated SNR was 0.63 overall, 0.62 for precipitation, and 0.58 for temperature proxies, suggesting that the prescribed SNR of 0.4 is a conservative estimate.The pseudo-proxies were also degraded to reflect the length, time span, and resolution of the real proxies.

Network ensemble creation
The creation of a network ensemble was done using a calibration-validation scheme to assess the ability of various pseudo-proxy combinations to reconstruct multi-decadal ENSO-like precipitation or temperature variations over the past millennium.As 20CRv2c only covers ∼ 160 years, it is not suitable for testing over these time spans and the past1000 runs from the GCMs were used instead.To account for the fact that the ENSO patterns in the GCMs are different from the pattern in 20CRv2c (and the real ENSO pattern is likely to be different again), two different GCMs are used for the calibration and validation stages.This means that the sensitivity of the networks to the shape of the ENSO-like pattern is tested.GCM selection is described in Appendix B2.The precipitation and temperature calibration GCMs (GCM cal ) were CCSM4 and ISPL-CM5A-LR respectively; the corresponding validation GCMs (GCM val ) were GISS-E2-R and BCC-CSM1-1.
Using the GCM cal past1000 dataset, 1000 proxy networks were created via an "add-one-in" pseudo-proxy algorithm that automatically builds a network based on how much each proxy improves the reconstructive power of the network (Fig. 1).Similar to a forward stepwise regression procedure, GCM cal -derived pseudo-proxies are gradually added to a base network of zero proxies, testing the quality of the network with each addition, until all proxies have been inwww.clim-past.net/13/267/2017/Clim.Past, 13, 267-301, 2017 corporated (Fig. 1a, b).The final optimal network is the one that performed best over all steps (Fig. 1c).By repeating the process 1000 times, adding different random noise series to the pseudo-proxies each iteration, it addresses the influence of stochastic processes on the ability of a proxy network to optimally reconstruct the large-scale ENSO pattern.
The reconstruction process itself is a weighted average approach, where the proxy weights are based on the ENSOlike empirical orthogonal function (EOF) pattern of the GCM cal past1000 run (for the add-one-in network building) or 20CRv2c (for the final reconstruction).First, all (pseudo)proxies in the network were normalised to a common period of at least 100 years (see Sect. 3.3) and transformed using an inverse transform sampling (ITS; Emile-Geay and Tingley, 2016) to approach a normal distribution.The EOF values w at the proxy locations i were scaled such that their absolute sum at each time step t is 1 (Eq.B3).This EOF scaling deals with the fact that the number of proxies available changes over time, preventing more proxy-dense periods from being amplified.Each proxy series p i was then weighted by its corresponding (scaled, time variable) EOF value w i , and the nweighted proxy series were summed to create a single reconstruction series ENSO r (Eq.B4).This is essentially a sparse reconstruction of the PC, which is the dot product of the full raw dataset and EOF.The quality of a network was assessed by comparing ENSO r with the principal component (PC) using the Pearson rank correlation R (see Appendix B4).
Validation was performed in two steps.First, the 1000 networks produced with the GCM cal were used to make reconstructions using GCM val past1000 data and the GCM cal EOF.Using data from a different GCM ensures complete separation between the calibration and validation periods and tests the sensitivity of the networks to the spatial stability of the EOF pattern as the ENSO-like EOFs from each model and from 20CRv2c are different.The switch from GCM cal to GCM val thus mirrors the switch between the GCMs and 20CRv2c.Each network was reconstructed 1000 times using GCM val pseudo-proxies, again adding different noise realisations for each iteration.Validation test scores were calculated to check the quality of these validation reconstructions compared to the validation PC (calculated using EOF cal and GCM val past1000 data).Second, critical values for the validation test scores were calculated by repeating the first validation step but using the GCM val piControl run.If the validation R value of a network failed to exceed R crit , its reconstruction was deemed to be no better than random noise and was consequently discarded.

Final reconstruction
The remaining networks were used to create an ensemble of ENSO r s using the real proxy data.The proxy records were first normalised to account for the different units (Eq.B2) and were subjected to an ITS transform before undergoing the same reconstruction process used in the network creation.All ) is selected and the associated proxy is moved from the test proxies to base.This is repeated until all test proxies are incorporated into base.(c) The optimal network is selected by cutting base where max n was highest.The entire process is repeated 1000 times, with new noise realisations being added to the pseudo-proxies at the start of each run.
ensemble members (ENSO r ) were then re-normalised to the reference period 0-650 yr BP to ensure comparability.The final ENSO reconstruction was taken as the ensemble mean, while the ensemble range represents part of its uncertainty envelope.
The last step in creating the reconstruction was calculating the final error range.RMSEs were calculated for each network during validation, providing 1000 error estimates for each ensemble member.The 95th percentile for each member was calculated from this and added and subtracted from the ENSO r series to find the maximum and minimum error limits.The uncertainty envelope around the final ENSO reconstruction (i.e. the ensemble mean) is thus a combination of the reconstruction ensemble range and the error ranges for individual ensemble members.

LIA-MCA difference analysis
The absence of a known reference period to which the reconstructions can be calibrated precludes any absolute comparison of the result with recent trends.However, it is possible to ascertain whether the MCA and the LIA differ significantly in how El Niño-like they are.Evaluating the LIA-MCA difference also directly removes the bias introduced by taking any reference period (Mann et al., 2009).To do this, the means over the two periods were taken and the MCA mean was subtracted from the LIA mean.If the difference is significantly greater than zero, the LIA is more El Niño-like than the MCA; if the difference is significantly less than zero, the MCA is more El Niño-like than the LIA.

Results: ENSO reconstructions
The final ENSO reconstructions for precipitation and temperature and their proxy network density are shown in Figs. 2  and 3.The final number of networks included in the precipitation ensemble is 1000, of which 999 are unique (Fig. 2a).The total number of proxies used is 48, with a maximum of 40 for a single network.Proxy availability increases steadily throughout time, save a slight drop off in the most recent period (Fig. 2b).Although there is spread in the ensemble, there are clear peaks and troughs visible.The within-ensemble coherence was tested by correlating 1000 randomly chosen pairs with each other.This confirmed that there is generally good agreement over the full period (100-1500 yr BP) as well as during the MCA and LIA individually (Fig. 4).
The temperature reconstruction ensemble consists of 182 optimal networks, all of which are unique (Fig. 3a).The total number of proxies is 211, with a maximum of 104 for a single network.Despite the higher number of proxies available for temperature, the median proxy coverage (Fig. 3b black line) is lower compared to the precipitation reconstruction; while roughly the most recent 1000 years of the precipitation reconstruction are based on a median of eight or more proxies, this is only true for the last ∼ 330 years of the temperature reconstruction.Prior to ∼ 480 yr BP, half the temperature ensemble members rely on a single proxy.This also accounts for the high within-ensemble correlations for the full and MCA periods but more variable correlations during the LIA (Fig. 4).With few long temperature proxies, most ensembles likely rely on the same data for the pre-LIA reconstruction.The steep increase in Fig. 3b.reflects the high number of tree ring and coral series, all but a handful of which are less than 600 years long (see Mann et al., 2008;Tierney et al., 2015) and most of which are clustered in North America.The addone-in method has mitigated some of the risk of co-varying non-white noise in a subset of the proxies skewing the resulting reconstruction; testing showed that when all North American tree ring records were added to the reconstruction, a regional climate trend obscured the ENSO-like signal (not shown).However, the relatively poor spatial coverage elsewhere and the lack of long proxies leaves the reconstruction prone to spurious noise-driven trends in the earlier period.The behaviour of the temperature ensemble members is less coherent and shows more spurious temporal variation compared to precipitation, particularly in the early period (Fig. 4).Nevertheless, for both temperature and precipitation the error from proxy noise is overshadowed by the uncertainty associ- ated with the choice of network -the ensemble spread makes up the bulk of the uncertainty envelope.
Figures 5 and 6 show the proxy locations plotted onto the precipitation and temperature EOF patterns respectively.The proxies included in the precipitation ensemble members are distributed well over the western and eastern side of the Pacific, though missing good coverage of the central Pacific.The relatively uniform size of the bubbles suggests that there is no immediate preference of any one proxy over the others.The spatial distribution of the temperature proxy locations, in contrast, is highly skewed towards North America, where most of the tree ring records are located, and the central and equatorial East Pacific lack coverage.The combination of this poor spatial coverage, low temporal coverage (Fig. 2b), and wide ensemble range leads to the expectation that the temperature reconstruction is of lower quality than www.clim-past.net/13/267/2017/Clim.Past, 13, 267-301, 2017    the precipitation reconstruction.This is further supported by the low fraction of networks that passed the validation (182 compared to 1000 for precipitation).There is no clear preference of any combination of proxies, with most proxies being selected equally often (i.e.equal bubble sizes).The fact that the similarity of the EOF patterns of GCM cal and GCM val to the 20CRv2c EOF pattern was lower for temperature than for precipitation (Appendix B2) further reduces confidence in the temperature reconstruction.
Figures 5 and 6 illustrate the benefit of using the pseudoproxy approach in creating the optimal networks.There is no direct correlation between proxy weighting (indicated by the bubble colour) and frequency of use, suggesting that other aspects such as resolution, length, and the relationship to other proxy locations played a significant role in determining the usefulness of a proxy that would be difficult to judge from the outset.The fact that the choice of proxy network is the dominant source of error is further evidence of the utility of the pseudo-proxy optimal network method.The high clustering of temperature tree ring records in North America is an example of where the add-one-in method has worked to reduce the risk that some co-varying non-white noise in a subset of the proxies skews the resulting reconstruction; testing showed that when all North American tree ring records were added to the reconstruction, a regional non-ENSO trend obscured the ENSO signal (not shown).

Comparing precipitation and temperature
Figure 7 shows the correlation between 1000 randomly chosen combinations of temperature and precipitation ensemble members as an indication of the agreement between the two climate variables.There is no correlation -positive or negative -apparent between the temperature and the precipitation reconstructions, neither over the entire 1500 years nor over the MCA or LIA individually.Whether this is a true physical phenomenon or simply a reflection of the high uncertainty in the reconstructions is difficult to separate.Therefore, it is not possible to categorically determine a systematic difference between the ENSO signals in temperature and precipitation proxies.
The definitions of the MCA and the LIA used here are based on those given by Yan et al. (2011b); there are many alternative definitions, however (Jones and Mann, 2004).To test the sensitivity of the results to the definition of these periods, we recalculated the LIA-MCA difference using two widely used alternative definitions: from Mann et al. ( 2009 , 2013). Figure 8 shows the range of LIA-MCA differences for the individual members within the precipitation and temperature ensembles, using the three definitions.Using the Yan et al. (2011b) definition, the precipitation interquartile range indicates that the LIA is more El Niño-like than the MCA, though the difference is statistically insignificant (p = 0.22).For temperature, there is no evidence of any difference between the ENSO-like state of the MCA and the LIA, with a median value very close to 0 (p = 0.48).There are also many more outliers (i.e.values outside the 95 % confidence interval) compared to precipitation, again reflecting the high uncertainty in the temperature reconstruction.Applying the definitions used by Mann et al. (2009) or Masson-Delmotte et al. ( 2013) only makes minor, and largely statistically insignificant, differences.For precipitation, the difference between the two periods is more pronounced for the alternative definitions, with a weakly significantly more El Niño-like LIA (p < 0.1).For temperature there is very little change; although the median is negative for the alternative definitions, the interquartile range still encompasses zero.The precipitation reconstruction thus qualitatively suggests that the LIA was more El Niño-like than the MCA, but our conclusion that there is no evidence for any precipitation-temperature correlation stands.

Discussion
An important question addressed in this study is whether the modern-day links between ENSO-like temperature and precipitation persist back in time.There is no concrete evidence in this study of any correlation between the precipitation and temperature reconstructions, whether positive (as seen in the modern day) or negative (as suggested by Yan et al., 2011b).This is contrary to expectations based on instrumental and modelling data, which both show a strong relationship between ENSO-like precipitation and temperature.To test the robustness of the ENSO-like temperature-precipitation coupling in GCMs, we calculated the correlation between the temperature and precipitation ENSO-like EOF time series (PCs) in the six CMIP5 GCMs listed in Appendix B2.
Four of the six models show a significant (p < 0.05) positive precipitation-temperature correlation over the study region (40 • S-40 • N) at annual and 30-year resolution for the historical run (0.74 ≤ R 2 ≤ 0.98), and five out of six for the past1000 run (0.30 ≤ R 2 ≤ 0.92; not shown).This is similar to coupling in the 20CRv2c data (R 2 ≥ 0.76).The fact that the palaeodata apparently do not display this relationship over the past 2 millennia (see Yan et al., 2011b, and this study) is thus interesting from a physical dynamical point of view as it contradicts our conventional understanding of long-term ENSO-like climate change.
There is also no evidence in the two reconstructions presented here that there was any significant difference in the mean ENSO-like climate state during the MCA and the LIA.This is also contrary to the findings of Yan et al. (2011b).They create a SOI reconstruction (SOI pr ) from two precipitation proxies from the Galápagos (Conroy et al., 2008) and the Indo-Pacific warm pool (Oppo et al., 2009), weighting them according to the relationship of local rainfall to the instrumental SOI.Interestingly, the SOI pr shows broad trends opposite to the precipitation reconstruction presented here, with a more La Niña-like LIA compared to the MCA.While the two proxies used were considered for this study, they were both rejected due to high dating errors (average around 100 years).Several other precipitation (Tierney et al., 2010;Yan et al., 2011b) as well as temperature (Conroy et al., 2010) proxies supporting conclusions of the Yan et al. (2011b) study were similarly rejected due to high dating errors, as were proxies supporting the opposite conclusion (Partin et al., 2007;Conroy et al., 2009).Testing showed that applying the method described here using only the two proxies used by Yan et al. (2011b) produced highly similar results to their SOI pr reconstruction, suggesting that it is not a methodological difference but rather it is related to proxy selection.A reconstruction based on only two (poorly dated) proxies is likely to be more vulnerable to spurious noise or other climatic influences distorting the signal; this is evidenced by the degradation back in time of the reconstructions presented in this study as the number of proxies declines.This highlights the need for more accurately dated proxy records, which remains an issue for low-resolution but long proxy archives such as marine sediments (Jones and Mann, 2004).
A number of other (non-temperature or precipitation) ENSO-sensitive proxies that were not included in our reconstructions provide evidence for a more La Niña-like climate state in the MCA compared to the LIA, although the mean state of the LIA appears inconsistent.Sedimentary sterol concentrations in marine sediment off the Peru coast (Makou et al., 2010) suggest the MCA coincides with a reduction in El Niño activity, with both El Niño and La Niña activity increasing from the late MCA onwards.Based on a range of North American proxies, Graham et al. (2007) conclude that the MCA was characterised by arid conditions in western North America consistent with a La Niña-like state, fol-lowed by a wetter LIA.A basin ventilation record from the Western Pacific warm pool (WPWP) (Langton et al., 2008) agrees particularly well with the earlier part of our precipitation reconstruction.It shows a peak in El Niño activity at ∼ 1150 yr BP and a distinctive minimum during the MCA, followed by a more El Niño-like LIA characterised by a steady decline in activity.This decline is not apparent in our reconstruction, but it is reflected in some other multimillennial proxy records (Moy et al., 2002;Stott et al., 2004;Conroy et al., 2008).
Most multi-proxy reconstructions of ENSO variability are temperature-based and focus on NINO regions (El Niño regions in the equatorial Pacific).A NINO3 region (90-150 • W and 5 • S-5 • N) temperature reconstruction by Mann et al. (2009) shows a slow millennial-scale warming trend (to a more El Niño-like state) from AD 1100 onwards, with relative cooling during the MCA compared to the LIA, consistent with a La Niña-like state during the MCA.In contrast, Emile-Geay et al. ( 2013b) are unable to detect a systematic difference between the MCA and LIA in their Boreal winter NINO3.4 (120-170 • W and 5 • S-5 • N) SST reconstruction, which is consistent with the findings of this study.The discrepancy between the two ENSO reconstructions may be due to the difference in proxy networks, particularly the use of lower-resolution proxies here and by Mann et al. (2009), which contribute a substantial part of the signal due to the slightly different definition of the NINO regions or the target season (Boreal winter versus annual).Other reasons may be related to the methodology or target instrumental dataset, particularly for low-frequency variability and amplitude.Work by Emile-Geay et al. (2013b) indicates that the results of many temperature reconstruction methods are sensitive to the target SST dataset used for calibration, and Wang et al. (2015) find that the La Niña-like pattern in the MCA evident in Mann et al. (2008) is not a robust feature across CFR methods.The fact that multi-proxy reconstructions are less likely to show strong differences between the ENSO-like state of the MCA and the LIA again highlights the potential sensitivity of individual records to non-physical trends, and suggests that conclusions drawn from single proxy records must be considered with caution.
The lack of a coherent and strong difference between the ENSO-like state of the MCA and the LIA in the reconstructions presented here and among other multi-proxy reconstructions may indicate that there is no direct relationship between long-term regional to global temperature anomalies and ENSO-like state of the climate.Conversely, this result may be an indication that these periods indeed were not characterised by a significant climate anomaly in ENSOsensitive regions.An examination of separate NH and Southern Hemisphere (SH) temperature reconstructions (Neukom et al., 2014) only indicates a global cool period between AD 1594 and 1677, while the PAGES2k Consortium (2013) report no globally synchronised multi-decadal climate shifts evident in their compilation of continental-scale temperature reconstructions except for generally cool conditions between AD 1580 and 1880.Both these cool periods coincide with the definition of LIA used in this study (AD 1400(AD -1850)).Neither study finds strong evidence for a global MCA, instead suggesting an initial NH warm period during AD 800-1100 (PAGES2k Consortium, 2013) followed by warming in South America and Australasia between AD 1160 and 1370 (PAGES2k Consortium, 2013) and the SH overall for AD 1200-1350 (Neukom et al., 2014).These continentalhemispheric warm phases are also encapsulated in the MCA definition used here (AD 800-1300).The lack of MCA or LIA indications in ENSO-like behaviour is partly consistent with the multi-proxy NINO3.4 reconstruction by Emile-Geay et al. (2013b), which finds no systematic difference between the two periods, and with the La Niña-like pattern during the MCA found in the multi-proxy reconstruction by Mann et al. (2009).Conversely, numerous reconstructions based on one or two proxy records (see Oppo et al., 2009;Yan et al., 2011a;Rustic et al., 2015) do report a significant difference, although as Yan et al. (2011b) highlight, the sign of the difference varies.This is another indication of the possible vulnerability of single-proxy reconstructions to over-interpretation, the danger of over-extrapolation, and "cherry-picking" of supporting evidence.The fact that multiproxy reconstructions are less likely to show strong differences between the ENSO-like state of the MCA and the LIA again highlights the potential sensitivity of individual records to non-physical trends and suggests that conclusions drawn from single-proxy records must be considered with caution.Further research into the physical and mechanistic interactions between global temperature and ENSO dynamics over long timescales is needed to elucidate this problem.Meanwhile, comparative studies of the MCA and the LIA must thus be approached with caution for low-latitude and SH study areas as their definitions are not necessarily rooted in real climatic events outside the NH.

Interactions with other climatic phenomena
There may be alternative interpretations of the reconstructions presented here due to the fact that ENSO interacts with various other climatic phenomena.The ITCZ is a zone of low-level atmospheric convergence that lies over the warmest water and creates substantial positive precipitation anomalies.It follows a seasonal migratory pattern, moving southwards towards the Equator during the boreal winter and retreating northwards during the summer months.As mentioned in the Introduction (Sect.1), there is a significant link between the movement of the ITCZ and ENSO state of the climate.This is thought to be true on interannual (Clement, 1999) to millennial timescales (Cane, 2005).On decadal to centennial timescales, proxy (Haug et al., 2001;Hodell et al., 2005;Sachs et al., 2009;Rustic et al., 2015) and modelling (Broccoli et al., 2006) evidence points to a potential NH temperature forcing on the ITCZ, causing it to shift southwards in response to cooler NH temperatures.This hypothesis has been invoked to explain an apparent southward shift of the ITCZ during the LIA (Hodell et al., 2005;Sachs et al., 2009;Rustic et al., 2015).Similarly, there is evidence for a northward shift of the ITCZ during the MCA (Wang et al., 2005;Zhang et al., 2008;Sachs et al., 2009;Oppo et al., 2009;Tierney et al., 2010).A southward shift is analogous to an El Niño-like state (Deser and Wallace, 1990;Haug et al., 2001;Ivanochko et al., 2005), suggesting that the MCA was characterised by a La Niña-like climate state and the LIA by an El Niño-like state.However, recent studies hint at an expansion and contraction of the ITCZ in the West Pacific during the MCA and LIA respectively (Yan et al., 2015;Griffiths et al., 2016), accompanied by a respective weakening and strengthening of the Pacific Walker Circulation (PWC).El Niño is characterised by a weakening in the PWC, which would imply a La Niña-like state during the LIA.
The most frequently used precipitation proxies (Fig. 5) are located in the West Pacific, which is an important region for ENSO, partly due to the relatively large latitudinal migration of the ITCZ.The precipitation reconstruction presented here suggests a slightly El Niño-like LIA (Fig. 8), which appears congruous with a southward shift of the ITCZ (and a weakening of the PWC), but not with a contraction of the ITCZ.However, inspection of the West Pacific proxy records indicates that all but two records are assigned negative EOF weights (i.e.their locations experience drought during El Niño), including the most frequently used and most heavily weighted proxies (Fig. 5).If the ITCZ did contract during the LIA, this could show up in the ensemble as an El Niñolike shift as these heavily weighted proxies would have experienced drought.Nevertheless, the ensemble is not solely based on these Western Pacific proxies, instead incorporating information from across the Pacific Basin.Thus, while it is possible that a contraction of the ITCZ influenced the results, this is (at least partially) mitigated by proxies in other locations.Further examination would be needed to conclusively determine the extent to which the ITCZ plays a role in driving the precipitation reconstruction here.
As mentioned in Sect. 1, the interpretation of the reconstructions is potentially confounded by interactions between ENSO and the PDO/IPO.The two (pseudo)-oscillations are known to share many similarities in their spatial pattern (Deser et al., 2016), and the temperature EOF used in this study (Fig. 6) is indeed very similar to the IPO pattern as reported in Henley et al. (2015).Many proxy records used here have decadal or lower resolution (Appendix A), and the 30-year averaging applied here brings the reconstructions into the realm of IPO/PDO variability.Moreover, many records, for example the North American tree rings, lie in PDO-sensitive North America (MacDonald and Case, 2005).This again raises the issue of how individual proxies are interpreted, potentially changing proxy-climate relationships at different timescales.Many proxy records of low resolution have been interpreted as ENSO reconstructions (see Moy et al., 2002;Nelson et al., 2011;Yan et al., 2011b) based on modern-day climate-proxy interactions, just as the reconstructions presented in this paper are based on modern-day (interannual) spatial ENSO patterns.
Comparison of the precipitation reconstruction in this study with a reconstruction of the PDO based on North American tree rings over AD 993-1996 (MacDonald andCase, 2005) shows a slight tendency for the precipitation and PDO series to have the same sign over the MCA and LIA separately, but the relationships are not statistically significant.There is no indication of any relationship between the PDO and the temperature reconstructions, despite the fact that many of the temperature proxies are located in potentially PDO-sensitive areas (most notably North America).MacDonald and Case (2005) find a strongly negative PDO during the MCA (roughly equivalent to a La Niña-like spatial pattern), which corresponds to the qualitatively La Niña-like tendency of the precipitation reconstruction presented here.D'Arrigo and Wilson (2006) find no significant correlation between Boreal winter NINO3 SST and a 9-year smoothed reconstruction of the Asian expression of the PDO (based on East Asian tree rings as opposed to the North American tree rings used by, for example, MacDonald and Case, 2005), suggesting that the ENSO-PDO link may be spatially variable.There is also no statistically significant correlation between the temperature or precipitation reconstructions and this Asian PDO reconstruction.The lack of co-variability between the reconstructions presented here and these PDO reconstructions may be due to the fact that this study relies on a wider proxy network that takes into account non-Pacific locations.Another possibility is that the IPO/PDO oscillation and spatial pattern is not robust further back in time, as suggested by Linsley et al. (2004) and Shen et al. (2006).The resolution of proxy data employed in this study prevents a robust separation of variability on ENSO and IPO/PDO timescales, leaving open the question of how these oscillations influence each other.
Not addressed in this study is the role of different "flavours" of ENSO patterns.A different type of ENSO pattern, first defined by Ashok et al. (2007) and dubbed ENSO Modoki, differs from the traditional (canonical) ENSO pattern in the shift of positive SST anomalies from the West Pacific (mainly in NINO3 and NINO3.4) to the central Pacific (NINO4; 170-120 • W and 5 • S-5 • N) and its mid-latitude teleconnections.It is sometimes defined as the EOF2 of detrended SST data (note that the data in this study were not detrended; hence, here it would be EOF3; Ashok et al., 2007;Cai et al., 2015a, b) or a combination of EOF1 and EOF2 (Takahashi et al., 2011;Karamperidou et al., 2015).Some modelling studies suggest that ENSO Modoki will increase in frequency compared to canonical ENSO as a result of anthropogenic climate change (Yeh et al., 2009;Kim and Yu, 2012;Cai et al., 2015a, b), and there is some model evidence that ENSO Modoki was also more common in the mid-Holocene (Karamperidou et al., 2015).The difference in equatorial spatial pattern and teleconnections has implications for the interpretation of the proxy reconstructions in this study as ENSO-Modoki-like climate change may appear here as a reduction in ENSO-like activity.
The poor quality of the temperature reconstruction, which limits the statistical robustness of the precipitationtemperature comparison, is likely due to the low number and unequal distribution of available data locations.Most temperature proxies are located in teleconnected regions outside the ENSO source region, which have been shown to be subject to more temporal variability in precipitation-temperature relationships (Wilson et al., 2010;Coats et al., 2013;Gallant et al., 2013;Lewis and LeGrande, 2015).A multi-region tree ring reconstruction of ENSO variability displays substantial variability in the strength of ENSO teleconnections over time and space (Li et al., 2013).The authors find that the Pacific Northwest and Texas-Mexico regions show highly unstable teleconnections (although there is no discussion on whether this is related to the different ENSO flavours).This may explain the lack of signal in the temperature reconstruction presented here as many of the temperature proxies are located in these teleconnected regions (Fig. 6).If the strength of the teleconnection has indeed changed over time, the weightings based on modern-day ENSO patterns would not reflect this; thus, this reconstruction should be regarded as an indication of change of the modern-day ENSO-like climate pattern only.Without proxies located in the centre of action or more robustly teleconnected areas, the loss of signal due to unstable teleconnections can be expected to be substantial, as suggested by the results presented here.
An interesting observation of the EOF maps presented here (Figs. 5 and 6) is the low to no correlation between EOF weighting (bubble shading) and how often proxies are used (bubble size).Correlations of temperature proxy frequency of occurrence versus GCM-derived and 20CRv2c-derived EOF weighting are 0.23 and 0.20 respectively (p < 0.05); there is no significant correlation for temperature proxies (p > 0.54).Two possible explanations for this are (i) that the climatic noise in the high-occurrence but low-weighted areas is less spatially correlated with the noise elsewhere than in the low-occurrence but higher-weighted areas or (ii) that the length and resolution of the proxy records have a more important effect on a proxy's utility than its weighting.For instance, the proxies off the Australian coast and in the western Pacific Islands are mostly short (< 500 years) coral records; while several of them have high weighting, their frequency of occurrence is very low.The proxy at the southern tip of Australia, in contrast, is a ∼ 3600-year-long tree ring record and is the most frequently used temperature proxy.Overall, however, temperature proxy length is only very weakly correlated to occurrence frequency (R = 0.17, p < 0.05).The SNR assigned to the temperature pseudo-proxies is similarly only very weakly correlated to frequency (R = 0.13, p < 0.05).In the precipitation case, neither length nor SNR are significantly correlated to how often the proxies occur in the ensemble.The selection process is likely driven by a combination of these factors rather than any single factor and is modulated by the number of proxies available.More detailed analysis will be needed to elucidate this.

Reflections on the method
The method set out in this study is one of few that attempt to take into account the effect of real spatial and temporal patterns of proxy records, thus increasing our confidence in their ability to accurately evaluate the effectiveness of the networks.To our knowledge, this is the only study in which realistic temporal proxy resolution in addition to their length has been taken into account by the pseudo-proxies.This is an improvement of the pseudo-proxy design used by, for example, Wang et al. (2014a), who take into consideration the declining proxy availability back in time but not their resolution.The authors find that this already significantly impacts the quality of multi-proxy reconstructions; thus, the inclusion of proxy resolution as done here is likely to have further impacts.More extensive research is needed to quantify this however.
The optimal network creation still has scope for improvement.Although we have screened for maximum dating errors, its effect on the included proxies is not explicitly assessed.This issue is often neglected in (multi)-proxy reconstructions (but see Comboul et al., 2014, for a recent effort to address it systematically).Moreover, the noise simulation used here is relatively simplistic; the use of a wider noise spectrum (including red, and possibly even blue, noise) may alter the composition of the networks (Smerdon, 2012, and references therein).However, the issue remains that there is no easy way to determine the real noise spectrum of the proxies.With the advent of more isotope-enabled GCM simulations, further improvement could come from the use of more proxy system models to more accurately estimate the proxy-climate relationships for all types of proxies (see Conroy et al., 2008;Evans et al., 2013;Russon et al., 2013;Stansell et al., 2013;Steinman et al., 2012;Sturm et al., 2010;Thompson et al., 2011;Tierney et al., 2011;Dee et al., 2015).
The choice of dataset from which to derive an EOF is also a source of uncertainty (see Emile-Geay et al., 2013b) as differences in the EOF pattern will affect the weighting of the proxies.This is particularly pertinent for the precipitation reconstruction as the modern-day ENSO precipitation signature is much less well-established than for temperature due to less and lower-quality instrumental data.This is partially tested by using different GCMs to calibrate and validate the proxy networks.However, the true ENSO-like pattern has been non-stationary over time, as has been shown to be true in 20CRv2 for the North Atlantic Oscillation (NAO) and Pacific North American pattern (PNA), for example (Raible et al., 2014).We tested the stability of the 20CRv2c temperature EOF used in this study by recalculating it for a running 30-year window and found substantial variability in the spa-tial pattern and amount of variance captured by the EOF.Further investigation is necessary to explore whether this result is an artefact of internal variability, is due to uncertainties in the reanalysis dataset, or reflects real changes in the nature of ENSO.Nevertheless, it highlights the vulnerability of the majority of ENSO reconstructions (including ours) to the assumption that the modern-day ENSO is a good analogue for the past.

Conclusions
Two reconstructions of ENSO-like climate change are presented based on temperature-and precipitation-sensitive proxies.The quality of the reconstructions degrades further back in time as there is less proxy data available, which is particularly detrimental to the temperature reconstruction.The main implications of these reconstructions are that 1. we find no evidence that temperature and precipitation proxies disagree over the ENSO-like state of the climate during the past 2 millennia.The two reconstructions in fact show little to no correlation, which is surprising as there is a strong relationship between temperature and precipitation ENSO behaviour at interannual timescales in instrumental-reanalysis data and GCMs.
2. the precipitation reconstruction shows a tendency for a more El Niño-like LIA compared to the MCA, but the difference is not statistically significant and is not apparent in the temperature reconstruction.This result is insensitive to the choice of definition for the MCA and LIA.
3. a major limitation on our ability to accurately reconstruct ENSO-like climate change back in time is the lack of high-quality long proxy records in the tropical and subtropical latitude bands, and we reiterate the need for continued efforts to collect such data.The discrepancies between the two series presented here and many other interannual and (multi)-decadal ENSO reconstructions are more likely to be reconciled with denser proxy networks in the ENSO source region, along with resampling of existing locations to increase the signal-tonoise ratio (Wang et al., 2014a).The pseudo-proxy experiments described in this paper can quite easily be adapted to search for optimal locations from which additional proxy information would be the most beneficial, as previously done specifically for corals by Evans et al.

Appendix A: Proxy data details
Tables A1 and A2 provide an overview of the proxy records collected for this study.Where a proxy was rejected, the reason is given.AOI refers to a proxy not being selected for any networks by the add-one-in algorithm.This could be due to poor ability to capture the EOF pattern related to location or time resolution.Records from the NOAA Paleoclimate Database are identified by original publications.Where there are multiple time series from one publication, an identifier suffix is added.The naming of this identifier is based on the naming in the original database files or the proxy type (e.g.δ 18 O, Sr / Ca).
Most tree ring records were taken from the dataset used by Mann et al. (2008), which is a reduced set derived from the International Tree Ring Data Bank (ITRDB, version 5.03; www.ncdc.noaa.gov/paleo/treering.html).The naming for these series has not changed from the original (an abbreviated location followed by a core number).The tree ring series were subject to the following selection criteria (Mann et al., 2008): "(i) series must cover at least the interval 1750 to 1970, (ii) correlation between individual cores for a given site must be 0.50 for this period, (iii) there must be at least eight samples during the screened period 1800-1960 and for every year used.Series that were redundant with other compilations [used in the Mann et al. (2008) study] were not included.Four other series were not included because of sitespecific problems [. . . ].Of the remaining series, [some] had to be replaced because of format errors in the chronology file on the ITRDB [. . .], or because sample depth data were missing from the chronology file.[. . .] When sample depth data were absent, the raw ring-width data from ITRDB were used to recalculate the chronology using program ARSTAN (Version 6.05P), with the following settings: (a) a single detrending fitting a cubic spline of 50 % variance reduction function at 66.6 % the length of each sample, no stabilisation of variance or autoregressive modelling, indices computed as ratio, that is measurement divided by curve, and chronology calculated as biweight robust mean estimate." The Tierney et al. (2015) coral dataset is a comprehensive compilation of coral data covering the last ∼ 400 years.As explained in Sect.3.1, only coral records other than δ 18 O are included here.The remaining records in this database were used in this study as temperature proxies, with the following exceptions: -DeLong et al. (2012,2013,2014) were replaced by the coral-derived SST series as presented in original publications.
-Goodkin et al. ( 2008) δ 18 O and Sr / Ca were replaced by the coral-derived SST series as presented in the original publication.
- Additional tree ring and coral records were retrieved from the NOAA Paleoclimatology Database (http://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets) and were used as presented there.Where available, temperature or precipitation reconstructions were used (i.e.where the raw proxy series has already been converted).This was done to minimise biases due to nonlinearities in the raw proxy data, which are accounted for in the conversion process.In some cases, two precipitation series were available for different seasons; these were summed (or averaged, depending on the type of data) to get a better approximation of an annual signal.While this may not be entirely accurate, annual signals are more desirable for the purpose of this study.Moreover, summing the records as opposed to treating them as individual records makes very little difference due to the nature of the method (weighting and summing the series).For corals, spliced records were used where available to maximise their length.
Proxy records with dating errors > 60 years were excluded from the analysis.The dating error was taken from the original publications where it was reported; otherwise it was derived from the age model results in the raw data files.In the latter case, the maximum error reported for the past 2000 years was used; larger errors further back in time were thus not taken into account as they are irrelevant for this study.
Table A1.Precipitation proxy details. "SNR" gives the signal-to-noise ratio used to make the pseudo-proxies."EOF" gives the (unscaled) 20CRV2c EOF value used to weight the proxy."Start" and "End" are starting and ending years in years BP."Res" is proxy resolution, rounded to the nearest integer; sub-annual proxies are listed as having a 1-year resolution."Dating" refers to dating error."Included?"indicates whether the proxy contributed to the final ENSO reconstruction; if not, the reason for exclusion is listed ("Dating error" and "Length" are a priori conditions, which the proxies failed to meet."AOI" indicates that it passed preprocessing screening, but was not selected during the add-one-in process).Table A2.Precipitation proxy details. "SNR" gives the signal-to-noise ratio used to make the pseudo-proxies."EOF" gives the (unscaled) 20CRV2c EOF value used to weight the proxy."Start" and "End" are starting and ending years in years BP."Res" is proxy resolution, rounded to the nearest integer; sub-annual proxies are listed as having a 1-year resolution."Dating" refers to dating error."Included?"indicates whether the proxy contributed to the final ENSO reconstruction; if not, the reason for exclusion is listed ("Dating error" and "Length" are a priori conditions, which the proxies failed to meet."AOI" indicates it passed preprocessing screening, but was not selected during the add-one-in process).Results for the highest-scoring PCs are shown in Table B1.Correlations with all three indices were extremely high for PC2 (R 2 > 0.62, p < 0.001); thus, EOF2 was selected as a basis for weighting the proxies.The precipitation and temperature EOF2 describe 11.57 and 9.75 % of variance respectively.These EOFs were used for the final reconstructions of ENSO-like climate change.

B2 GCM data
For the creation of the network ensemble, GCM data were employed.The objective of the add-one-in method is to create networks that accurately reconstruct the long-term (30-year averaged) ENSO signal.20CRv2c covers less than 200 years, which is too short for meaningful evaluation of low-frequency change, particularly if the dataset is to be partitioned into calibration and validation sets.GCMs, meanwhile, offer much longer datasets.Although they cannot simulate the real temporal climate change and variability of the past 1000 years, they are still useful for the building of proxy network ensembles, which asks only that they simulate realistic modes of spatio-temporal variability (i.e.EOF patterns).For the pseudo-proxy experiment, the past1000 runs from two different GCMs were used for calibration (GCM cal ) and validation (GCM val ).
The following six GCMs (followed by their climate modelling groups in brackets) were considered as they have past1000 runs available: All model datasets were regridded to 2 • × 3 • .Each GCM produces a different ENSO-like EOF pattern, with varying biases (Bellenger et al., 2014;Collins, 2005); there is further variation among the different model runs.For this study, the most accurate GCM past1000 runs were selected in the sense that their EOF values at proxy locations were most similar to the corresponding 20CRv2c values, which is the most realistic EOF pattern of modernday ENSO available here.The real modern ENSO pattern will be different again, as will the real ENSO pattern of the past 1-2 millennia.By calibrating and validating the networks on datasets with slightly varying realisations of this ENSO pattern, the sensitivity of the networks to these variations is tested.EOF analysis was performed separately on the piControl and past1000 runs of the six GCMs for precipitation and temperature.For each analysis, the first three EOFs were retained for comparison to the 20CRv2c ENSOlike EOF.The GCM EOF values of grid boxes with proxies were compared to the corresponding 20CRv2c grid boxes via Pearson rank correlation.The two GCMs with the highest correlations with 20CRv2c across the two runs were chosen separately as GCM cal and GCM val for precipitation and temperature.As a result, CCSM4 (R ≥ 0.79 in piControl and past1000) and GISS-E2-R (R ≥ 0.75) were chosen as GCM cal and GCM val for precipitation; ISPL-CM5A-LR (R ≥ 0.46) and BCC-CSM1-1 (R ≥ 0.44) were chosen as GCM cal and GCM val for temperature.

B3 Pseudo-proxies
Since there is no straightforward way of assessing the quality or relevance of a proxy beyond the selection criteria already discussed, a pseudo-proxy approach can aid in making a more objective and refined decision on how to optimise the use of available proxy data (Smerdon, 2012).Since for pseudo-proxies the real world is known (in this case the GCM-derived EOFs and PCs), it is possible to quantify the skill of the reconstruction.While a blanket approach (in which every available record is used) may sound attractive, it increases the risk that some co-varying non-white noise in a subset of the proxies skews the resulting reconstruction.This is, for example, pertinent in North America where there is high clustering of tree ring records.Testing showed that when all records were added to the reconstruction, a regional trend obscured the ENSO-like signal (not shown).The only case in which it is certain that using all available proxies is preferred is when each grid box contains a (non-climatic) noise-free proxy such that the network gives complete spatial and temporal coverage.
Pseudo-proxies are model or instrumental data series that are degraded by applying transformations and/or adding noise to simulate the behaviour of a real proxy (Mann, 2002).Proxy system models (Evans et al., 2013;Dee et al., 2015) attempt to characterise the mechanistic and forward processes connecting the response of a proxy archive to an environmental signal and the subsequent observation of this response; this includes accounting for nonlinearities, multivariate responses, and measurement limitations.Many PSMs require an isotope-enabled GCM or other data not available in the original sources of the proxy data collected for this study, but the VS-Lite model is an easily implemented PSM for simulating tree ring widths (TRW) that needs minimal input (Tolwinski-Ward et al., 2011;Breitenmoser et al., 2014;Evans et al., 2014).Using the R package VS-LiteR available from GitHub (https://github.com/suztolwinskiward/VSLiteR), GCM precipitation and temperature data were combined to create pseudo-TRW datasets, with which TRW proxies were represented in the pseudoproxy experiment here.The rest of the proxies were represented by the GCM raw temperature or precipitation series.
These precipitation, temperature, and TRW series (taken from the real proxy locations) were firstly degraded by adding white Gaussian noise.Information on true signal-tonoise ratios (SNRs) of the various proxies is sparse in the literature.However, where data were available on the relationship between the proxy and the target climate variable in the form of R values, SNR can be calculated via the following equation (Smerdon, 2012) Several original publications on individual publications provide such information (see Tables A1 and A2), and Tierney et al. (2015) have conducted a systematic comparison of most coral records included here against the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST) 1900-1990 dataset.Where no concrete information was available, a value of SNR = 0.4 was prescribed.While this has been shown to be a realistic average (Smerdon, 2012), comparison with the calculated SNRs in this study suggests that it is in fact a relatively conservative value.After the addition of noise, the pseudo-proxies were degraded further to reflect the length, time span, and resolution of the real proxies.Real proxy temporal resolution between 0 and 1000 yr BP was applied; this period was chosen as it is the focus of this study.The resolution was recreated by assuming that each data point represents an average of the previous unsampled years; for example, in a proxy with a 10-year resolution, data point n was recreated by taking the average over points (n − 9) : n.

B4 Calibration: network creation
After screening proxy records, selecting GCMs, and creating pseudo-proxies, a pseudo-proxy experiment was conducted to create an optimal network ensemble.The calibration stage builds proxy networks in a stepwise manner by incrementally adding the proxy that maximally improves the quality of the network reconstruction.All proxies were first subjected to inverse transform sampling (ITS) using the MATLAB script (translated to R) available at https://github.com/CommonClimate/common-climate/blob/master/gaussianize.m.ITS converts the distribution of a time series to a standard normal distribution, thus reducing the bias introduced by nonlinearities inherent to many proxy records (Emile-Geay and Tingley, 2016).Lastly, the pseudo-proxies were averaged to 30 years (using a simple average) to prevent high-resolution records from dominating the signal.This averaging period was chosen since our focus is on long-term ENSO state rather than inter-annual variability, and it reflects the resolution of many individual low-frequency proxy reconstructions (see Yan et al., 2011a, b;Anderson, 2012;Rodysill et al., 2012;Makou et al., 2010).Then at each stage (for each interim network), the pseudo-proxies in the network (derived from GCM cal past1000) were first normalised to a base period (n) where all proxies overlap via the equation where z i is the normalised proxy value at time t, c t is the original value at time t, µ n is the mean of series c over the base period n, and σ n is the standard deviation over the same base period.The length of n was ensured to be at least 100 years to reduce the influence of noise on µ and σ .In some cases, this requirement led to the rejection of one or more proxies because their length or position in time were not compatible with the other records.This process of normalisation is similar to the method used in Wilson et al. (2010).Proxies that fall in the same grid box were averaged after normalisation.This prevents over-representation of those locations and improves their SNR by cancelling out some of the stochastic noise and amplifying their signal (Wang et al., 2014a).Each normalised pseudo-proxy z at location [x, y] was multiplied by a scaled version of the EOF value at location [x, y].The scaling was such that at each time step the absolute sum of the weights was 1, which accounts for the fact that the number of locations with proxy data varies over the reconstructed period (especially the beginning and end).The weight w at proxy location i can thus be considered an indication of the relative sensitivity of location i to the ENSOlike EOF pattern, which changes depending on the number of proxies n available at each time step t and their associated EOF values d: The resultant vector w i is thus the weighting series for proxy location i.The n proxy series are multiplied by their corresponding w i and are finally summed to create a single reconstruction series (ENSO r ): This is essentially a sparse reconstruction of the PC, which is the dot product of the full raw dataset and EOF.The quality of a network was assessed by comparing ENSO r with the PC using the Pearson rank correlation R.
The entire calibration process was repeated 1000 times, each time using pseudo-proxies with newly generated noise iterations, resulting in 1000 proxy networks.

B5 Validation: network evaluation
The validation stage tests the robustness of the networks using independent data.Of each calibration network, 1000 pseudo-proxy versions were created using GCM val pseudoproxies and were tested for their ability to reconstruct the target (the PC created using the GCM cal EOF and the GCM val past1000 dataset) by calculating R.There are several measures of quality for validation statistics, the most common of which are the coefficient of determination R 2 , the reduction of error (RE) and the coefficient of efficiency (CE; Cook et al., 1994).Discussions on the relative merits and pitfalls of these measures can be found in the literature (see Cook et al., 1995;Emile-Geay et al., 2013a).Although CE is generally regarded as the appropriate indicator for low-frequency reconstructions (Bürger, 2007;Emile-Geay et al., 2013a), the nature of the method described here reduces its effectiveness as a measure of quality.The past1000 and piControl runs have little (or no) external forcing driving the simulation; hence, they have very little low-frequency variability or very few trends.Moreover, the data are z-normalised at various stages, removing any differences in means.As the CE effectively tracks changes in the mean, this removal of the mean renders the CE sensitive to spurious results.The diagnostic R was instead chosen as comparison of the three indicators, showing that it was more effective at picking high-quality reconstructions than the CE and RMSE, though generally a high R value did correspond to high CE and low RMSE (not shown).R is essentially equivalent to using R 2 , but it retains the ability to distinguish between positive and negative correlations.
The R values from the 1000 × 1000 validation reconstructions (R val ) were then compared to critical values R crit calculated for each network.Again, 1000 reconstructions were made for each network, but using GCM val piControl data instead of past1000.The piControl run contains no external forcing and is therefore essentially noise, but it retains the inherent climatological spatial correlations.From these reconstructions, R crit for each network was determined by taking the 95th percentile of the R values of the corresponding 1000 reconstructions.Where R val > R crit , the network was retained; where R val < R crit , the network was deemed unfit and was discarded.Networks sensitive to the choice of dataset are thus weeded out.
The combination of using pseudo-proxies, the add-one-in approach and R crit simultaneously accounts for proxy temporal resolution, spatial distribution, and temporal coverage (i.e.proxy start and end dates) and gives an estimate of the uncertainty due to proxy noise.However, an important assumption is that the signal in all proxies is solely temperature or precipitation, and it is thus still a best-case estimate.

B6 Proxy ENSO ensemble
The optimal networks that passed the R crit test were used to create an ensemble of real proxy reconstructions of ENSOlike climate change.The remaining networks may not all be unique, further reducing the effective number of networks.Presumably, networks that occur multiple times are more effective proxy combinations; retaining the duplicates accordingly upweights these networks in the final reconstruction.Error estimates for the reconstructions were made using the 1000 × 1000 RMSE values calculated at the validation stage.This was translated into ensemble member uncertainty limits by adding and subtracting the 1000 error series from the reconstruction time series (to get the maximum and minimum error limits respectively) and taking the 5-95th percentile over their full range.This error estimation explicitly takes into account the impact of network choice as well as random error affecting the proxies.
Once proxy reconstructions and associated uncertainty estimates were calculated for all ensemble members, they were renormalised to 100-650 yr BP to make the trends and amplitudes comparable within and between temperature and precipitation.The period 100-650 yr BP was chosen because it was common to all reconstruction time series and only covers one of the two periods of interest (the LIA).Although calibration to the instrumental period would potentially allow us to quantify the absolute amplitude, this was not done for two reasons.Firstly, the proxy data coverage during the instrumental period and the preceding century was relatively low, reducing the confidence in the reconstruction during that period; calibrating to this period would thus increase the uncertainty on the rest of the reconstruction.Secondly, any calibration to the instrumental data is necessarily biased towards high-frequency trends (Mann et al., 2008).Within a 30-yearaveraged series, the number of comparison points with the instrumental period is extremely low.The final proxy reconstruction was calculated as the ensemble mean.The corresponding error estimate is a combination of the reconstruction ensemble range and the error ranges for individual ensemble members.

Figure 1 .
Figure 1.Network creation process diagram.Overview of the network creation process.(a) A new network is created from the base network (base) plus each pseudo-proxy individually and is tested for its reconstructive power.This results in n test scores; (b) the highest score (max n) is selected and the associated proxy is moved from the test proxies to base.This is repeated until all test proxies are incorporated into base.(c) The optimal network is selected by cutting base where max n was highest.The entire process is repeated 1000 times, with new noise realisations being added to the pseudo-proxies at the start of each run.

Figure 2 .
Figure 2. Precipitation ENSO ensembles.(a) Precipitation reconstruction of ENSO-like climate change that has been 30-year averaged (black line).Individual network solutions are shown as orange lines, with the uncertainty envelope in orange shading.(b) Number of proxies included in the ensemble over time, with the median in black and the range in blue.The pink and purple shaded periods are the MCA and LIA respectively.

Figure 4 .
Figure 4. Within-ensemble correlations.Correlations between 1000 pairs randomly chosen from the precipitation (blue) and temperature (red) ensembles.Box plots encapsulate the space between the first and third quartile with the median shown as a black line; whiskers indicate the 95 % confidence interval of the median; points are values outside this confidence interval (outliers).Statistical significance of the median value is indicated at the bottom: * * * p < 0.001, * * p < 0.01, * p < 0.05, ^p < 0.1, p > 0.1.

Figure 5 .
Figure 5. Precipitation EOF with proxy locations.Background colours are scaled EOF values.Bubbles are individual proxies; size is indicative of how often the proxy is included in the network ensemble, and shading indicates relative weighting such that darker colours are more strongly positive or negative.

Figure 6 .
Figure 6.Temperature EOF with proxy locations.As Fig. 5, but for temperature.

Figure 7 .
Figure 7. Precipitation-temperature correlations.Correlations between the temperature and precipitation ensemble members based on 1000 randomly chosen pairs for the period 100-1500 yr BP ("All") and the MCA and LIA individually.See Fig. 4 for explanation of the box plots and significance.
Figure7shows the correlation between 1000 randomly chosen combinations of temperature and precipitation ensemble members as an indication of the agreement between the two climate variables.There is no correlation -positive or negative -apparent between the temperature and the precipitation reconstructions, neither over the entire 1500 years nor over the MCA or LIA individually.Whether this is a true physical phenomenon or simply a reflection of the high uncertainty in the reconstructions is difficult to separate.Therefore, it is not possible to categorically determine a systematic difference between the ENSO signals in temperature and precipitation proxies.The definitions of the MCA and the LIA used here are based on those given byYan et al. (2011b); there are many alternative definitions, however(Jones and Mann, 2004).To test the sensitivity of the results to the definition of these periods, we recalculated the LIA-MCA difference using two widely used alternative definitions: fromMann et al. (2009) (MCA = AD 950-1250, LIA = 1400-1700) and the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (MCA = AD 950-1250, LIA = 1450-1850;Masson-Delmotte et al., 2013).Figure8shows the range of LIA-MCA differences for the individual members within the precipitation and temperature ensembles, using the three definitions.Using theYan et al. (2011b) definition, the precipitation interquartile range indicates that the LIA is more El Niño-like than the MCA, though the difference is statistically insignificant (p = 0.22).For temperature, there is no evidence of any difference between the ENSO-like state of

Figure 8 .
Figure 8. Difference between the means of the MCA and LIA for three definitions.Difference calculated by subtracting µ MCA from µ LIA for each ensemble member, using the three MCA and LIA definitions listed in Sect. 4. A positive value indicates LIA is more El Niño-like than MCA.Precipitation is on the left in blue, temperature is on the right in red."Yan" refers to the definition used in Yan et al. (2011b), "Mann" in Mann et al. (2009), and "IPCC" in Masson-Delmotte et al. (2013).See Fig. 4 for explanation of the box plots and significance.
Felis et al. (2009) Sr / Ca,Kuhnert et al. (2005) Sr / Ca,Quinn et al. (1996),Quinn et al. (2006) Sr / Ca, and Kilbourne et al. (2008) Sr / Ca were excluded as their correlations with SST reported inTierney et al. (2015) were of opposite sign to what is expected (i.e.positive when the physical processes should lead to a negative correlation).
Henke et al.: ENSO reconstructions -MIROC-ESM (Atmosphere and Ocean Research Institute (The University of Tokyo), National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology).