Stability of ENSO and its tropical Paciﬁc teleconnections over the Last Millennium

Determining past changes in the amplitude, frequency and teleconnections of the El Niño–Southern Oscillation (ENSO) is important for understanding its potential sensitivity to future anthropogenic climate change. Palaeo-reconstructions from proxy records provide long-term information of ENSO interactions with the background climatic state 5 through time. However, it remains unclear how ENSO characteristics have changed through time, and precisely which signals proxies record. Proxy interpretations are un-derpinned by the assumption of stationarity in relationships between local and remote climates, and often utilise archives from single locations located in the Paciﬁc Ocean to reconstruct ENSO histories. Here, we investigate the stationarity of ENSO teleconnec-10 tions using the Last Millennium experiment of CMIP5 (Coupled Model Intercomparison Project phase 5) (Taylor et al., 2012). We show that modelled ENSO characteristics vary on decadal-to centennial-scales, resulting from internal variability and external forcings, such as tropical volcanic eruptions. Furthermore, the relationship between ENSO conditions and local climates across the Paciﬁc basin varies throughout the 15 Last Millennium. Results show the stability of teleconnections is regionally dependent and proxies may reveal complex changes in teleconnected patterns, rather than large-scale changes in base ENSO characteristics. As such, proxy insights into ENSO likely require evidence to be synthesised over large spatial areas in order to deconvolve changes occurring in the NINO3.4 region from those pertaining to proxy-relevant local 20 climatic variables. To obtain robust histories of the ENSO and its remote impacts, we recommend interpretations of proxy records should be considered in conjunction with palaeo-reconstructions from within the Central Paciﬁc.


Introduction
The El Niño-Southern Oscillation (ENSO) is an important determinant of climate variability, altering global rainfall patterns and modulating global temperatures. Understanding the long-term characteristics of ENSO variability and its sensitivity to external forcings, such as greenhouse gases, represents a fundamental climate modelling and 5 data challenge. Changes in ENSO behaviour may occur under future global warming (Power et al., 2013). However, there is a is a large dispersion in global climate model (GCM) projections of changes in ENSO characteristics (e.g. Collins et al., 2010;Vecchi and Wittenberg, 2010), and hence the sensitivity of the coupled ocean-atmosphere system to future changing boundary conditions remains uncertain (DiNezio et al., 10 2012). In addition, investigations of the sensitivity of ENSO to anthropogenic climate change are restricted by the relatively short instrumental record, which provides us with limited guidance for understanding the range of ENSO behaviours. For example, the observed changes in the character of ENSO since the mid-1970s towards a dominance of El Niño, rather than La Nina, episodes are difficult to evaluate in terms of a forced 15 response or unforced variability given the limited observational record almost certainly does capture the full range of internal climate dynamics.
High resolution palaeo-reconstructions, including from tree rings, sediment cores, corals and speleothems, have the potential to provide long-term information about changes in modes of climatic variability and their sensitivity to different boundary con-20 ditions. Some tropical proxy records reveal ENSO interactions with the background mean climatic state. For example, data from long-lived fossil corals are often interpreted quantitatively as estimates of ENSO changes through time that show a range of ENSO frequencies and amplitudes through time. Central Pacific coral reconstructions generally reveal a weakened ENSO during the early Holocene (McGregor et al., 25 2013) and highly variable ENSO activity throughout the Holocene (Cobb et al., 2013), which may have arisen from internal ocean-atmosphere variability (Cobb et al., 2003). Developing robust estimates of natural ENSO variability over a period longer than per-1581 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | mitted through the instrumental record is a useful research avenue, with the potential for informing meaningful adaptive strategies for future climate change.
Palaeo-ENSO proxy records of the Last Millennium (1000 years) are sparsely populated temporally and spatially, and reconstructions remain uncertain (Cobb et al., 2003;Khider et al., 2011). It also remains unclear as to precisely which climatic signals as- 5 sociated with ENSO are being recorded in these individual proxy records and whether these provide the necessary resolution to reconstruct ENSO changes. The assumption of stationarity of relationships between local and remote climates (teleconnections) underpins the interpretation of many palaeoclimate reconstructions, although stationarity should not necessarily be assumed in terms of ENSO variability (Gallant et al., 10 2013). Are palaeo-reconstructions from the tropical Pacific recording base changes in the ENSO system or rather changes in teleconnected patterns? Previous modelbased studies have identified sensitivity in the relationship between ENSO and the background climate state, and urged caution in the reconstruction of ENSO from proxy records under the assumption of stationarity of observed teleconnections (Coats et al.,15 2013; Gallant et al., 2013).
However, these studies have not comprehensively addressed the degree to which uncertainty about the non-stationarity of ENSO teleconnections can be assessed for particular locations and for particular mean climatic states. Furthermore, although we previously investigated the potential non-stationarity of hydrologic responses to  like conditions under disparate boundary conditions in idealised model simulations, we did not provide specific guidance for interpreting tropical proxy records (Lewis et al., 2014), which currently comprise our dominant source of information about ENSO characteristics beyond the instrumental record. In addition, while previous studies have utilised proxy records, together with simulations using global climate models (GCMs) to evaluate the representation of ENSO in the current generation of GCMs (Cobb et al., 2013), these approaches focused on using palaeo-ENSO reconstructions to test the performance of GCMs for the purpose of constraining uncertainty in future projections of ENSO behaviour under climate change.
1582 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | As such, precisely which expressions of ENSO are being recorded in proxy archives under differing climatic boundary conditions have not been comprehensively interrogated. Climate models, in addition to observational and proxy climate evidence, allow an understanding of long-term ENSO changes through time to be obtained (Schmidt, 2010). A new generation of climate models and experiments has recently become 5 available (Taylor et al., 2012), providing an opportunity for the first time to investigate ∼ 1200 years of ENSO variability and establish a framework for understanding ENSO changes through time, using more models than previously possible. Hence in this current study, we investigate changes in ENSO characteristics (frequency and amplitude) in model experiments of the Last Millennium ("past1000"). Focusing on three 10 key climatic regions (East, Central and West Pacific), where explicit palaeo-ENSO reconstructions have been made, teleconnected patterns (the relationship between local and remote climates) throughout the Last Millennium are examined for surface temperatures and precipitation. We ultimately aim to determine whether proxy archives in the tropical Pacific are likely to be recording alterations in ENSO base frequencies or 15 local-scale teleconnections under differing boundary conditions.

Definitions
The study is primarily focused on palaeo-ENSO variability from the tropical Pacific. Model data were investigated in three regions that have been identified as sensitive to 20 modern ENSO variability and have also been used explicitly to reconstruct past ENSO changes (e.g. Cobb et al., 2013;McGregor et al., 2013). Area-mean anomalies for precipitation and surface temperature were calculated for the West ( These regions are not intended to provide exhaustive coverage of ENSO impacts, but are large enough to provide useful comparisons with model-based data. El Niño episodes were defined based on surface temperature anomalies in the NINO3.4 region, with events defined in the models when NINO3.4 temperature anomalies were > 0.5 K for at least six consecutive months (Trenberth, 1997). Conversely, La 5 Niña episodes were defined when NINO3.4 temperature anomalies were < −0.5 K for at least six consecutive months. Spatial patterns are examined by compositing temperature and rainfall anomalies into positive (El Niño) and negative (La Niña) phases using these definitions for all CMIP5 models analysed. We utilise the NINO3.4 region as an index to classify ENSO conditions. Although the NINO3.4 region is commonly used to 10 categorised ENSO episodes, it should be noted that there are other indices of ENSO that may also provide useful information beyond the Central tropical Pacific conditions described by the NINO3.4.

Models and experiments
CMIP5 data (Taylor et al., 2012) were downloaded from the Project for Model Diagno-15 sis and Intercomparison (PCMDI) through the Earth System Grid (ESG). Simulations were used of the historical (1850-2005 CE) experiment, which is forced using changing atmospheric compositions due to observed anthropogenic and volcanic influences, solar forcings and emissions of short-lived species from natural and anthropogenic aerosols. In addition, simulations were used of the Last Millennium (past 1000) (850-20 1850 CE), in which reconstructed time-evolving exogenous forcings are imposed, including changes in volcanic aerosols, well-mixed greenhouse gases, land use, orbital parameters and solar changes. Each model's pre-industrial control simulation (piControl) with non-evolving pre-industrial forcings was analysed. Models were analysed where past 1000 simulations were provided (Table 1). One model (bcc-csm1-1) was excluded from analysis because its representation of ENSO spectra is too short, producing a biennial oscillation rather than a broad spectral peak in the 2-8 year band. The MIROC-ESM model was also excluded, as it exhibits large drift related error in the 1584 form of long-term trends that cannot be attributed to natural variability (Gupta et al., 2013). Further details are provided in Sect. 3 on Model Evaluation. Data (precipitation (pr) and surface temperature (ts)) for six remaining models were regridded onto a common 1.5 • latitude by 1.5 • longitude grid. For the piControl and past 1000 experiments, monthly anomalies were calculated by subtracting the mean sea-5 sonal cycle for each model. For the historical experiment, only model years 1976-2005 are considered, as this experiment necessarily produces a non-stationary climate due to the time-evolving anthropogenic greenhouse gas forcings imposed. Additional experiments were analysed for CMIP5-participating models, where available. For GISS-E2-R (Schmidt et al., 2014) and IPSL-CM5A-LR (Dufresne et al., 2013) models, extended control simulations of > 500 years in duration were analysed and compared to forced, past1000 experiments.

Diagnosing ENSO changes and teleconnections
The location of ENSO activity in the historical and Last Millennium experiments in was first explored using the leading empirical orthogonal function (EOF) of the tropical Pa-15 cific surface temperature. These spatial patterns were compared to the NINO3.4 index to determine possible non-stationarities in the site of ENSO activity through time (Li et al., 2011). This EOF analysis (Fig. S1) demonstrates that in both experiments, the surface temperature patterns are loaded in the NINO3.4 region. Although there are some differences in the spatial patterns of the leading EOF mode across the equato-20 rial Pacific, the similarity in model experiments in this particular region indicates that areal-average NINO3.4 temperatures provide a useful metric of ENSO activity in both experiments.
Hence, a wavelet analysis was next used to examine the frequency and amplitude of NINO3.4 surface temperature variability in each model for statistically significant 25 changes. Wavelet analysis is useful for examining non-stationary signal and provides time and frequency localisation. A Morlet mother wavelet (Torrence and Compo, 1998) with degree 6 was used to calculate the wavelet power spectra and identify large-1585 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | scale changes in variance. Wavelet spectral estimates were tested against red noise, represented as a first order autoregressive process. First, the wavelet power spectrum of NINO3.4 surface temperatures was considered, with variance in the wavelet power spectrum plotted as function of time and period (Fig. 2a). Next, the NINO3.4 mean wavelet power spectrum, generated using a Morlet wavelet of 6 • , was used as a metric 5 for ENSO amplitude. The spectral power was calculated for the historical simulation and compared to the range of spectral power displayed in the past1000 experiment, calculated using ten 100 year epochs (Fig. 2b). The relationship between ENSO variability and teleconnected patterns in the tropical Pacific regions (East, Central and West) was diagnosed through several complemen-10 tary approaches. First, an ordinary least squares regression between monthly NINO3.4 mean surface temperature and remote area-mean surface temperature, and between monthly NINO3.4 mean surface temperature and remote area-mean precipitation was compared for the historical and Last Millennium experiments, for each region. Second, the relationship between local and NINO3.4 climates was considered using the corre-15 lation between variables (Corr(Local, Remote)), analogous to considering land-surface coupling strength (Lorenz et al., 2012). Correlations coefficients were calculated for monthly timeseries in ten 100 year epochs comprising the Last Millennium. Values were determined at each model gridbox and an area-weighted mean calculated for each region. The significance of correlations was assessed at the 95 % confidence 20 level for each coefficient using a t test. Third, the significance of identified changes in local-remote relationships during the Last Millennium was investigated.
For each 100 year epoch comprising the Last Millennium, the El Niño-and La Niñaassociated local temperature and precipitation anomalies were selected for each region. A two-sided Kolmogorov-Smirnov (KS) test was used to investigate whether the 25 distribution of local climate variables in 100 year epochs within the Last Millennium could statistically have been drawn from the same population (at the 5 % significance level). A two-sided KS test was applied to each ENSO phase for each variable (surface temperature, precipitation) in each region (East, Central, West) comparing every permutation of epochs sequentially (e.g. comparing El Niño-associated Central Pacific temperatures during 850-949 with 950-1049, then 1050-1149, then 1150-1249 etc.). A KS test was used for detecting changes in ENSO-remote climate relationships in Last Millennium timeseries as it is non-parametric and requires no assumptions to be made regarding the distribution of the data. A change is detected where the null hypothesis 5 (that the distributions considered were drawn from the same population) is rejected at the 5 % significance level.

Model evaluation
The basic properties of El Niño-Southern Oscillation (ENSO) simulated in Coupled Model Intercomparison Project phase 5 (CMIP5) models (Taylor et al., 2012) have been 10 comprehensively examined, relative to observations, in previous studies (Bellenger et al., 2013;Guilyardi et al., 2012). Here, ENSO was examined through 6 metrics −1) ENSO amplitude (Niño3 sea surface temperature (SST) SD), (2) structure (Niño3 vs. Niño4 amplitude), (3) frequency (root mean square error of Niño3 SST anomaly spectra), (4) heating source (Niño4 precipitation SD), (5) the amplitude of the ENSO bien-15 nial component (the ratio of the Niño3 SST anomaly timeseries power in the 3-8 years and 1-3 years bands) and (6) seasonality of ENSO (ratio between winter November-January over spring March-May average Niño3 SST anomalies SDs). Notably, some models have a tendency to represent ENSO with a too short period of about 2 years (bcc-csm1-1, Fig. S2) or to have a spectral peak displaced towards longer periods of 20 between 3 and 8 years (CSIRO-Mk-3.5 and MIROC5). Models that have contributed Last Millennium simulations to CMIP5, but are identified as erroneously biennial in terms of ENSO spectral power are excluded from this analysis. In addition, long-term drift in long model simulations must be considered. The MIROC-ESM model is excluded from this analysis as it exhibits large drift related error in the form of long-term 25 trends that cannot be attributed to natural variability, but instead relate to deficiencies in model physics and numerics (Gupta et al., 2013) (Fig. 3).

1587
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The remaining six models that contributed CMIP5 Last Millennium simulations (Table 1) were compared to twentieth century reanalysis data (20CR) (Compo and Whitaker, 2011) as a proxy for the observed climate that has been widely used to investigate observed ENSO-rainfall relationships (King et al., 2015;Klingaman and Woolnough, 2013). The CMIP5 models were compared to reanalysis precipitation and 5 surface temperature for several ENSO-related characteristics over the common period, derived from the CMIP5 "historical" simulation. First, to investigate the model representation of ENSO spatial patterns, the first empirical orthogonal function of the tropical Pacific surface temperature anomalies was calculated for 20CR reanalysis and CMIP5 multi-model mean (MMM) EOF ( Fig. 4a and b) and similarly for precipitation anoma-10 lies ( Fig. 4c and d). This demonstrated the surface temperature and precipitation patterns are qualitatively similar for reanalysis and models; surface temperature patterns are generally of the same sign, although the meridional width of tropical temperature anomalies is narrower than in the reanalysis estimates, and simulated precipitation patterns are similar to the reanalysis estimate in the Central Pacific, although posi-15 tive anomalies are located too far westward in the CMIP5 MMM, compared with observed. In addition, the relationship between NINO3.4 surface temperature anomalies and global precipitation fields in reanalysis was compared to the CMIP5 MMM ( Fig. 4e and f). The correlation coefficients between NINO3.4 temperature anomalies and local precipitation are generally of the same sign in simulated and reanalysis fields, includ- It should be noted, however, that EOF analysis does not necessarily reveal modes that can be readily interpreted physically. In addition, only the MMM patterns are for each experiment, possible biases in ENSO representations in CMIP5 models are not considered prohibitive to investigating changes in the stability of teleconnections through time. 5 Models demonstrate significant variance in the ENSO-relevant band (2-8 years) throughout the Last Millennium simulation (Fig. 2a). Several models also display significant low-frequency variability beyond the 8 year threshold; for example, the FGOALS-s2, GISS-E2-R and MRI-CGC3 models exhibit multi-decadal scale variability in NINO3.4 temperatures during the Last Millennium simulations. These decadal-to 10 centennial-scale El Niño-and La Niña-like episodes during the Last Millennium simulations are evident in all models analysed here (Fig. 3). Models also demonstrate a range of ENSO amplitudes (Fig. 2b), with low spectral power shown, for example, in the GISS-E2-R and MRI-CGMC3 (Stevenson, 2012) models. There are differences of ENSO amplitude between the historical and Last Millennium experiments within 15 individual models. In the historical, the amplitude is generally weaker at all periods for MRI-CGMC3 and GISS-E2-R. The amplitude of higher ENSO-relevant periods (5-8 years) in the historical is largely outside the range exhibited in the Last Millennium (Fig. 2b). Previous model-based studies have also revealed strong inter-decadal to inter-centennial modulation of ENSO behaviour (Coats et al., 2013;Wittenberg, 2009) 20 and warn that such modulation may not be fully revealed by the comparatively short instrumental climate record available and hence could attach large uncertainties to ENSO metrics diagnosed from short records.

ENSO characteristics
Low frequency modulation of Last Millennium ENSO characteristics may result from internal variability, or they may be externally driven. External forcings, specifically large 25 tropical volcanic eruptions occurring between 1250 and 1600 CE (Fig. S3), may pro-1589 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | duce decadal-to centennial-scale ENSO responses. Across the model ensemble investigated here, La Niña-like post-volcanic responses appear to be simulated on multidecadal timescales (Fig. 3). However, the response of ENSO to explosive volcanism is uncertain and may be dependent on the strength of the forcing (i.e. the size of the eruption). On seasonal to annual timescales, model evidence suggests the radiative forcing 5 due to volcanic stratospheric aerosols induces a La Niña episode that is followed by an El Niño episode after the peak of the forcing (McGregor and Timmermann, 2011). Proxy reconstructions derived from tree rings across the Pacific reveal similar ENSO responses to those simulated, with anomalous cooling reconstructed in the East-Central tropical Pacific in the year of volcanic eruption, followed by anomalous warming occurring one year after . This pattern of response is likely sensitive to the specific set of eruptions investigated (Wahl et al., 2014). Previously, the association of eruptions and subsequent El Niño episodes has been demonstrated for forcings larger than that observed during the historical period for Mt Pinatubo (Emile-Geay et al., 2008). By comparison, when considering the strongest eruption of the Last Mil-15 lennium in 1258 CE, the eruption likely produced a moderate or strong El Niño episode that interrupted the prevailing La Niña-like conditions induced by increased solar activity during the Medieval Climate Anomaly (Fig. S4). Low solar activity during the Little Ice Age (1450-1850 CE) has been linked to generally cooler conditions, in contrast to a possible solar forcing of warmth during the Medieval Climate Anomaly (950-1250 CE) 20 (Schurer et al., 2014). Analysis of the frequency characteristics of various solar forcings imposed on the models during the Last Millennium simulations (Fig. S5) shows some decadal-and centennial-scale power, but this is not qualitatively consistent with variance in simulated ENSO conditions (Fig. 2).
Alternatively, decadal-to centennial-scale modulation of ENSO behaviour may result from internal ocean-atmosphere dynamics , rather than a response to exogenous forcings. Previous studies suggest low-frequency ENSO characteristics may not be produced in long control simulations (Landrum et al., 2013), which exclude external climate forcing. When the properties of ENSO simulated in the control simulations where compared to Last Millennium CMIP5 experiments, the smaller subset of extended control simulations available and analysed here exhibit qualitatively similar variability to that shown in the externally-forced Last Millennium experiment . For the GISS-E2-R (Schmidt et al., 2014) and IPSL-CM5A-LR (Dufresne et al., 2013) models, extended unforced pre-industrial control ("piControl") simulations 5 of > 500 years in duration were provided as part of CMIP5. The smoothed mean annual timeseries of NINO3.4 region temperature anomalies during the control simulations (Fig. 5) exhibit qualitatively similar variability to the forced Last Millennium simulations (Fig. 3), including multi-decadal to centennial-scale El Niño-and La Niña-like phases. Using wavelet analysis (Torrence and Compo, 1998), the amplitude and frequency of ENSO variability is investigated throughout the piControl simulations. The models reveal some differences in low-frequency variability beyond the 2-8 ENSO-relevant band, when comparing the control (Fig. 6) and the Last Millennium (Fig. 2). In GISS-E2-R, for example, there is weaker spectral power exhibited on multi-decadal timescales in the control experiment, compared with the Last 15 Millennium, although some low frequency modulation is apparent in both models due to internal variability. This limited suite of CMIP5 models hints that ENSO variability on timescales that are poorly captured in the short instrumental record results from internal climate variability, and hence occurs in both piControl and Last Millennium experiments.

20
Overall, there may be important difference in ENSO characteristics in the Last Millennium experiment, compared to the historical and piControl simulations, relating to both external forcings from large tropical volcanic eruptions or internal variability revealed on longer than inter-annual timescales. Differences between the piControl and Last Millennium experiments may represent an artefact of comparing a limited number of model Rather, we suggest that these differences highlight the considerable natural variability in the ENSO system that cannot be determined entirely from centennial and shorter records (Wittenberg, 2009). Proxy reconstructions providing "snapshots" of climate his-5 tory throughout the Last Millennium are often too short to resolve unforced and forced changes and these brief snapshots may reflect sampling of specifics ENSO characteristics within a short window. Longer reconstructions of spliced coral records, however, suggest ENSO responds to both low-frequency internal variability and external forcings (Cobb et al., 2003;McGregor et al., 2013).

ENSO impacts and teleconnections
Models show broadly similar global impacts associated with NINO3.4 regional temperature anomalies in the Last Millennium and historical experiments (Fig. 1). Composited patterns of global surface air temperature anomalies associated with positive (El Niño) and negative (La Niña) ENSO phases derived from all analysed models are spatially 15 coherent across the experiments. Also, the Last Millennium patterns largely resembling those in the historical experiment. Although ENSO surface temperature anomalies across the Pacific are qualitatively similar, anomalies associated with the historical period  are generally of greater magnitude, particularly at remote locations outside the equatorial Pacific, including over North America and the south Pa-20 cific. These differences in magnitude between the Last Millennium and the historical may relate to the differing boundary conditions during the historical period associated with anthropogenic forcings, such as long-lived greenhouse gases, or simply from the greater diversity of ENSO episodes represented in the longer Last Millennium simulation. 25 Hence, we next explore the temporal stability of teleconnected ENSO patterns and ENSO impacts in different domains during the Last Millennium. The relationship between NINO3.4 regional temperature anomalies and the mean local climate is exam-1592 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ined in each analysed Pacific region (East, Central, West). Overall, the linear relationship between local and NINO3.4 surface temperature for the East Pacific and Central Pacific regions is consistent in the Last Millennium and historical experiments (Fig. 7). However, when considering data from the analysed models for the West Pacific region, this linear relationship is different for these experiments. Also, there is a statistically sig-5 nificant difference between experiments in the relationship between local precipitation and NINO3.4 conditions in the East Pacific. Complexity in the relationships between local (surface temperature and precipitation) and remote (NINO3.4 surface temperature) conditions during the Last Millennium simulation is further considered using the correlation between variables (Corr(Local, Remote)). This approach is analogous to considering land-surface coupling strength (Lorenz et al., 2012). We diagnose temporal stability using this correlation in ten 100 year epochs that comprise the Last Millennium (Fig. 8). The strength of the remotelocal relationship varies temporally and is regionally and climate variable dependent. In the West Pacific, particularly, this coupling is generally weak and varies significantly 15 between epochs. For local precipitation, the strongest coupling occurs for the Central Pacific.
We also investigate the significance of these identified Last Millennium changes in local-remote relationship. A Kolmogorov-Smirnov (KS) test was used to determine whether the distributions of El Niño-and La Niña-associated local temperature and 20 precipitation anomalies in each region in 100 year Last Millennium epochs could statistically have been drawn from the same population. There are detectable differences (at the 5 % significance level) in the distribution of ENSO-associated local climate variables in these 100 year epochs. West Pacific El Niño-and La Niña-associated temperatures, for example, significantly vary in character through the Last Millennium. Tempo-25 ral changes in local ENSO fingerprints (Corr(Local, Remote)) of the Last Millennium, also likely result from external forcings and/or internal ocean-atmosphere dynamics. However, these same relationships were not explored in the extended control simulations because of the small number of contributions available from different models.

1593
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Regardless of attributing the source of ENSO variability in the Last Millennium simulations to external forcings or internal dynamics, it is evident in the model experiments that differing teleconnections may result at different points in time and may differ from present-day relationships. In addition, Last Millennium variability in ENSO-local climate relationships across sites the tropical Pacific suggests that global ENSO changes do 5 not necessarily scale linearly to local scales and cannot be assumed to do so.

ENSO under differing boundary conditions
The CMIP5 archive also provides simulations of the mid-Holocene (midHolocene, circa 6000 years ago) from multiple participating climate models. The mid-Holocene provides a well-constrained target for model-based studies (Schmidt et al., 2004) with 10 substantially larger time-evolving forcings than those imposed during the Last Millennium, and this period has also been the target of palaeo-reconstructions. Hence, these simulations are also briefly investigated here, in addition to the information provided by the Last Millennium experiment. Mid-Holocene simulations are run for at least 100 years after reaching equilibrium and have changed orbital parameters and atmo-15 spheric concentrations of greenhouse gases imposed. Other boundary conditions such as aerosols, solar constant, vegetation and topography are prescribed as the same as in the pre-industrial control simulation.
By way of context, Cobb et al. (2013) report that Central Pacific corals record highly variable ENSO activity through the Holocene, although no systematic trend in ENSO 20 variance was demonstrated in this study. A complementary Central Pacific reconstruction from Kiritimati Island suggests that ENSO variance was persistently reduced by 79 %, compared with today at this location about 4300 years ago (McGregor et al., 2013). Central Pacific coral-based evidence of ENSO variability is substantially different from lower-resolution records from the eastern equatorial Pacific (Conroy et al.,25 2008; e.g. Moy et al., 2002). Collectively, East Pacific records suggest a systematic decrease in mid-Holocene ENSO variance. On the West Pacific side of the basin, corals from northern Papua New Guinea reveal a reduction in ENSO frequency and ampli-tude over the period of 7.6-5.4 ka (thousand years ago) compared with today, and also identifies large and protracted El Nino events for 2. 5-1.7 ka (McGregor and Gagan, 2004). These Mid-Holocene ENSO reconstructions do not necessarily provide contradictory information, but may instead reflect geographic complexities (Carre et al., 2014;Cobb et al., 2013). However, as proxy-based reconstructions from each of these 5 regions have been used to infer changes in the same coupled ocean-atmosphere system, we also examine teleconnected ENSO patterns under these significantly different boundary conditions that characterise the mid-Holocene.
In this study, we consider the subset of participating CMIP5 models with contributions of mid-Holocene simulations (MRI-CGCM3, IPSL-CM5A-LR, FGOALS-s2, CCSM4) 10 and find a general reduction in spectral power across ENSO-relevant frequencies (Fig. S6a) that has also been reported in model experiments of this period conducted prior to the release of CMIP5 (Chiang et al., 2009). This reduced spectral power in the ENSO band can be considered a metric for reduced ENSO amplitude (Stevenson, 2012). Previous model and proxy-based studies have also hinted at subdued ENSO 15 activity in the mid-Holocene. For example, early studies using simple numerical models of the coupled ocean-atmosphere system by Clement et al. (2000) demonstrate increasing ENSO variability throughout the Holocene in response to time varying orbital forcings. The impact of mid-Holocene orbital changes on ENSO variability has not been demonstrated comprehensively from proxy records. However, various fossil coral 20 reconstructions indicate that there may have been reductions in ENSO variability in the mid-Holocene (Cobb et al., 2013).
In addition, when CMIP5 midHolocene model data are composited into positive (El Niño) and negative (La Niña) phases, the magnitude of simulated mid-Holocene spatial patterns of ENSO impacts (Fig. S6b and c) are subdued, relative to the historical. 25 The relationship between NINO3.4 mean surface temperature anomalies and regional (East, Central, West Pacific) temperature and precipitation was examined (see Fig. S7 for West Pacific). This indicates that the relationship between West Pacific surface temperature anomalies and corresponding NINO3.4 temperature anomalies differs from 1595 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the midHolocene and historical simulations. The frequency of high and low local surface temperature anomalies in the West Pacific during El Niño defined conditions is also reduced in the midHolocene experiment compared with the historical (Fig. S7). The West Pacific surface temperature anomaly distributions for both El Niño and La Niña phases are notably narrower for the midHolocene, relative to the historical, indi-5 cating a reduced magnitude of ENSO impacts in this region. Conversely, the precipitation distributions for the West Pacific are for the most part similar for each experiment. The NINO3.4 impacts on East and Central Pacific regional temperatures are broadly similar for the historical and mid-Holocene. 10 Palaeoclimate studies typically utilise archives from single locations located in the Pacific Ocean to reconstruct generalised basin-scale histories of ENSO. However, we demonstrate that proxies in one location alone should not be considered regionally representative, or singularly insightful about robust ENSO reconstructions without explicit examination of the stability of ENSO teleconnections. Using standardised simula-15 tions performed with multiple participating climate models, we have demonstrated that there is model evidence for decadal-to centennial-scale ENSO variability during the Last Millennium that may not be fully revealed by the comparatively short instrumental climate record. The existence of varying ENSO characteristics throughout the Last Millennium is also supported by proxy-based climate reconstructions (Cobb et al., 2003), 20 which show variable ENSO characteristics include changing frequency and amplitude compared to modern during the Last Millennium.

Towards reconstructing robust ENSO histories
Furthermore, we find ENSO-local climate relationships vary significantly throughout the Last Millennium in the various tropical Pacific regions. In addition to (1) the sensitivity of ENSO to the background climate state and (2) the decadal-to centennial-scale 25 modulation of ENSO behaviour from internal variability and external forcings such as volcanic eruptions, this suggests that (3) changes in ENSO do not necessarily scale linearly to local scale impacts. Models suggest it may be inherently difficult to deconvolve variability in the NINO3.4 region and local-scale, teleconnected climatic change in the West and East Pacific regions. The West Pacific Warm Pool is likely sensitive to subtle shifts in the western extent of the warm (cool) tongue characterising positive (El Niño) and negative (La Niña) episodes. Hence, subtle changes in the Pacific basin 5 may impact this region through several ocean-atmosphere mechanisms.
We argue that proxy insights into change and variability in ENSO system are likely to be most robust when evidence is be synthesised over large spatial areas. That is, only incomplete information about temporal changes in a large-scale climate system can be provided by considering changes at a singular location (i.e. a time series of a climatic variable). Rather, considering multi-dimensional information in the form of spatial patterns of change through time is likely to yield more robust insights in large-scale systems. This provides a framework for enhanced interpretations of the invaluable information of palaeoclimatic change provided by proxy records. For example, combined evidence from the West and Central Pacific is more likely to reveal the potentially sub-15 tle changes in ENSO-associated spatial patterns of temperature and precipitation perturbations across the Pacific. For remote regions outside the equatorial Pacific, the non-stationarity of ENSO teleconnections is likely to be more problematic. These sites should be considered in conjunction with palaeo-reconstructions from within the Central Pacific basin, the so-called "centre of action" of ENSO (Cobb et al., 2013). We do 20 not expect that this recommendation is dependent on the NINO3.4 metric used to define ENSO utilised here. Indeed, such recommended spatially integrated approaches have already been undertaken and provide valuable information over the recent past (e.g. Li et al., 2013). It is also worth noting that under boundary conditions significantly different from present, such as the mid-Holocene, the stability of ENSO teleconnections 25 is likely to be more variable, and hence potential non-stationarities in local-remote relationships require explicit consideration in proxy interpretations.
This study has used palaeoclimate simulations conducted using a suite of CMIP5participating models with various forcing to investigate changes in ENSO and its 1597 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | teleconnections under differing boundary conditions (the Last Millennium and mid-Holocene). Such palaeo-modelling approaches do not attempt to replace proxy reconstructions, but rather combining multiple approaches can provide enhanced interpretations of reconstruction of past climate guiding our understanding of the most consistent physical explanations for change (Schmidt, 2010). In this study, we considered 5 a subset of CMIP5 models that contributed palaeo-simulations and these necessarily contain systematic biases in ENSO representations (Power et al., 2013). In their study focused on understanding ENSO responses to volcanic forcings, Emile-Geay et al. (2008) suggested further forcing/response insights could be provided by GCMs with realistic ENSO cycles and asked whether the current generation of models were 10 up to the task. Deficiencies in our theoretical knowledge of ENSO and the difficulties in representing physically realistic ENSO cycles in GCMs (Guilyardi et al., 2012) are a limit on providing robust quantitative understanding of forced and unforced changes in the ENSO system. Nonetheless, existing model simulations remain useful for examining palaeoclimates, despite their biases. These reveal spatially and temporally 15 complex changes in ENSO and its teleconnected patterns during the Last Millennium. These complex spatial and temporal changes should be considered when developing robust proxy interpretations and ENSO histories in order that these are most useful for constraining future ENSO behaviour under greenhouse forcings.
This study has identified several avenues for further research. First, additional tar-20 geted experiments within a single climate model that well represents ENSO spatial dynamics, particularly on the western extent of the warm/cold tongue, would provide further insight into the apparent complexity of ENSO impacts through time. In particular, explicit investigation of the sensitivity of changes in the location of maximum base ENSO activity across the Pacific to varying boundary conditions would provide use- 25 ful information about the range of variability of this coupled system. Furthermore, our present study did not comprehensively investigate the relative influences on various external forcings (solar and volcanics) and internal variability on ENSO characteristic, which would provide useful information for comparison with proxy records.
More direct comparisons between model output and proxy reconstructions can be provided by using a larger set of hydrological diagnostics and by employing pseudoproxy techniques. In the first approach, the relationships between local climate parameters and circulation, beyond the instrumental period, are examined using global climate models. Incorporating water isotope tracers into global climate models as a di-5 agnostic tool brings output one step closer to direct model-proxy comparisons. Ongoing research using water isotope tracers in GISS ModelE-R in experiments of the historical and Last Millennium period will provide further location-specific information about water isotope proxies and ENSO variability. Future studies focusing on water isotope variability may provide useful insights into changing ENSO impacts across the 10 Pacific as these hydrological diagnostics have been shown to integrate regionally coherent climatic information, compared with precipitation, which is highly variable (Lewis et al., 2010). Furthermore, in the second approach, a simulated time series intended to mimic actual proxy records ("pseudo-proxy") is generated from a climate model simulation (Anchukaitis and Tierney, 2013). The pseudo-proxy approach can be used to 15 interrogate the necessary proxy density required for producing skilful regional climate field reconstructions and provide guidance on interpretations of reconstructions from particular locations (Smerdon, 2011;Wahl et al., 2014).  Historical past1000 Historical past1000 Historical past1000 Historical past1000 Historical past1000 West Pacific Figure 7. Scatter plots of mean surface temperature anomalies (K) in the NINO3.4 region against area-mean anomalies in surface temperature (left) and precipitation (right) for the East (a, d), Central (b, e) and West (c, f) Pacific. Data points are shown for the historical (grey) and past1000 (blue). Lines of best fit are calculated for each experiment using ordinary least squares regression.