Stability of ENSO teleconnections during the last millennium in CESM

The El Niño-Southern Oscillation (ENSO) has a significant impact on the global climate through atmospheric teleconnections. It is important to understand the stability of ENSO teleconnections, not only for future weather forecasting and climate projection, but also for ENSO reconstructions based on paleo-proxies. In this study, we investigate the decadal variations of ENSO teleconnections in global land surface temperature (LST) from 850 to 2005AD using 13 ensemble members of the Community Earth System Model-Last Millennium Ensemble (CESM-LME). The CESM can simulate the main Eurasian cooling and Arctic warming, known as the warm Arctic-cold Eurasia (WACE) pattern, during the boreal winter of an El Niño. Furthermore, it can also capture the western Antarctic warming during the developing and decaying summers of an El Niño. There is a dominant decadal variation in the ENSO-LST teleconnections, expressed as anomalous LST patterns that closely resemble those seen in the WACE pattern during boreal winter and the western Antarctic warming pattern during summer. This decadal variation of ENSO-LST teleconnections is primarily due to the varying positions of Rossby wave sources associated with distinct ENSO patterns, which are located either to the west or to the east of Hawaii. The LST response to ENSO over South Siberia, as well as the associated precipitation response over North Eurasia, even show opposite patterns at different phases of the decadal variation. The decadal variation in CESM is found to be related to the interdecadal Pacific oscillation (IPO) and is likely attributed to internal variability rather than external forcing. Our findings suggest that the decadal variation in ENSO teleconnections should be considered when using proxies from Eurasian regions to reconstruct ENSO variability.


Introduction
The El Niño-Southern Oscillation (ENSO) is considered one of the most important interannual variabilities of Earth's climate system, causing global impacts on precipitation and temperature over both ocean and land (Mcphaden et al. 2006;Guilyardi et al. 2009;Collins et al. 2010;Deser et al. 2010;Andréa et al. 2021).The sea-surface temperature (SST) anomalies in the central-to-eastern equatorial Pacific associated with ENSO events can change tropical climate by modifying the Walker circulation (Bjerknes 1969;Gill 1980;Ropelewski and Halpert 1987;Klein et al. 1999;Alexander et al. 2002;Meyers et al. 2007).The ENSO also modifies extratropical temperature and precipitation through atmospheric teleconnections, i.e., via the well-known Pacific-North American (PNA) and Pacific-South American (PSA) teleconnections.The El Niño will cause cooling over Eurasia and North America against the background of global warming (Hoskins and Karoly 1981;Wallace and Gutzler 1981;Mo and Ghil 1987;Lau and Nath 1994;Hoerling et al. 1997;Bladé et al. 2008).Additionally, the ENSO also has substantial impacts on polar regions, and an El Niño event can induce warming in the Arctic and Antarctic by weakening the polar vortex (Domeisen et al. 2019;Li et al. 2021).
Understanding the stability of ENSO teleconnections is essential for ENSO reconstructions based on paleoclimate records.Paleoclimate data, such as tree-ring chronologies, corals, and sediment records, can provide insights into past ENSO variability (Gergis et al. 2006;Christie et al. 2009;Carr et al. 2013).Li et al. (2011) reconstructed a long-term ENSO index for 1100 years using tree rings from Southwest America, which was based on the PNA teleconnections of ENSO.Undoubtedly, these paleoclimate reconstructions rely on the assumption that ENSO teleconnections remain stable over the past few centuries.However, recent studies have revealed that ENSO and its regional teleconnections can vary on interdecadal to centennial timescales (Diaz et al. 2001;Ault et al. 2013aAult et al. , 2013b;;Tejavath et al. 2019), and some studies have cautioned that stationarity should not be assumed when reconstructing ENSO variability from proxy records (Coats et al. 2013;Gallant et al. 2013; Lewis and LeGrande 2015;Hope et al. 2017).
Previous studies have found that the relationship between ENSO and tropical oceans exhibits interdecadal variability (López-Parages et al. 2016;Park and Li 2019;Liu et al. 2021).The impact of El Niño on the warming of North Tropical Atlantic SST is dependent on the phases of the Atlantic multi-decadal oscillation (AMO), and is typically intensified during the negative phase of the AMO (Park and Li 2019;Wang et al. 2017).The relationship between ENSO and Indian Ocean SST also presents a decadal variation, depending on the strength of the linkage between decadal variabilities in the Pacific Ocean and Indian Ocean (Xie et al. 2010;Du et al. 2013;Tao et al. 2015;Ham et al. 2017;Kim and An 2019;Liu et al. 2021).Moreover, the ENSO-North Pacific teleconnection is significantly stronger in late-20th century than that in mid-20th century (Christopher 2018).In addition to the oceans, a non-stationary signature of ENSO in the atmosphere has also been detected.Recent studies indicate that there has been a substantial interdecadal change in the relationship between ENSO and the North Atlantic Oscillation (NAO) over the past century (Knippertz et al. 2003;Greatbatch et al. 2004;Rodríguez-Fonseca et al. 2016).Furthermore, the once robust relationship between ENSO and the Asian summer monsoon has broken down in recent decades (Kumar et al. 1999;Zhou et al. 2007;He and Wang 2013;Geng et al. 2018), and the impact of ENSO on rainfall in Europe has become more pronounced in the 20th century than in previous periods (López-Parages and Rodríguez-Fonseca 2012).
The causes of decadal-to-interdecadal changes in ENSO teleconnections are still under ongoing discussion.Some works suggest that the diversity of the ENSO, such as the amplitude and spatial structure of SST anomalies, could potentially affect the teleconnection patterns (Murphy et al. 2014;Ashok et al. 2007;Rodrigues et al. 2015;Zhang et al. 2015).For instance, the ENSO-Indian Ocean teleconnection has become stronger since 1970s, primarily due to the increased ENSO intensity which enhances the local thermocline feedback (Chowdary et al. 2012;Xie et al. 2016).Low-frequency multidecadal variability of the background state, such as the AMO and interdecadal Pacific oscillation (IPO), can also modulate the decadal change of ENSO teleconnections (Wang et al. 2008;Chen et al. 2010;Watanabe and Yamazaki 2014;Geng et al. 2017;Feng et al. 2021;Li et al. 2021;Liu et al. 2021).According to Geng et al. (2020), ENSO events are considerably weaker during the positive AMO phase than during the negative phase.Power et al. (1999) also reported that the IPO in its negative phase tends to increase the frequency of La Niña events, which show more significant correlations with Australian precipitation as compared to El Niño events.
Besides the internal variability in the climate system, volcanic eruptions as the strongest external natural forcing were found to not only increase the likelihood of El Niños (Adams et al. 2003;Emile-Geay et al. 2008;Timmreck 2012;Stevenson et al. 2016;Liu et al. 2022), but also affect the decadal phase change of the AMO (Knudsen et al. 2014) and the IPO (Maher et al. 2014) Hence, the modulations of decadal variations of ENSO-teleconnection may be attributed to internal variability, external forcing, or both.
The non-stationarity of ENSO teleconnections poses as a challenging problem for reconstruction of past climate in both tropical and extratropical regions.Therefore, in order to determine whether ENSO teleconnections based on recent observations can be employed to reconstruct the ENSO index during the last millennium, it is essential to conduct an investigation into the stability of ENSO teleconnections.The long-term outputs from ocean-atmosphere coupled numerical models have become valuable tools for studying the decadal changes of ENSO teleconnections, since the limited observation datasets are too short to analyze the decadal variations (Christopher 2018;Yu et al. 2018;King et al. 2020).
This study aims to detect non-stationary behaviors of atmospheric teleconnections associated with ENSO during a millennium period and to explore their decadal-to-multidecadal modulators.The following two questions will be addressed in this paper: (1) What are the main features of non-stationary teleconnections during the last millennium?(2) What factors modulate the variation of teleconnections on the decadal-to-multidecadal timescale?The paper is organized as follows.The datasets and analysis methods used in this study are described in Sect. 2. In Sect.3, we present the main results in detail.Summary and discussion will be given in Sect. 4.

Data and method
The Community Earth System Model-Last Millennium Ensemble (CESM-LME) Project is an ensemble of global Earth System Model experiments spanning the period of 850-2005AD (Otto-Bliesner et al. 2016).These outputs have been widely used to analyze characteristics of paleoclimate (Abram et al. 2016;Shi et al. 2018;Stevenson et al. 2019;King et al. 2020;Liu et al. 2022).The monthly outputs used in this study come from 13 ensemble members under a full set of transient forcing (solar variability, volcanic eruption, land use, greenhouse gases, orbital change, and ozone aerosol) at a resolution of 2°×2° for the atmosphere and land, and 1°×1° for the ocean and sea ice.
The averaged SST in the Niño3.4region (5° S-5° N; 120° W-170° W) in the boreal winter (December to February) is defined as a proxy of the ENSO.The other fields such as land surface temperature (LST), precipitation, geopotential height (GPH) at 200 hPa and global SST were seasonally averaged over May to September (MJJAS) and November to March (NDJFM), providing the boreal-summer and borealwinter fields.
To extract the decadal variability of ENSO teleconnections, we performed an empirical orthogonal function (EOF) analysis of 11-year-sliding regression maps for all 13 ensemble members spanning the period of 850-2005AD.First, we calculated the regression maps of LST onto the ENSO for each 11-sliding years with an interval of five years, resulting in 229 sets of 11-year-sliding regression maps for each of the 13 ensembles.The climatological regressions of LST onto ENSO were then obtained by averaging all of the 229 × 13 maps.The anomalous fields were calculated by subtracting the climatological regression from each regression map.Afterward, we applied the EOF analysis to these 229 × 13 anomalous maps to obtain the decadal changes.To study the seasonal evolution from ENSO developing summer to decaying summer, a seasonal EOF analysis (Wang and An 2005) was performed by combining the three regression maps for developing summer, peaking winter, and decaying summer of ENSO into one regression map.The principal components (PCs) were normalized by their standard deviations, which were then projected onto the corresponding EOF patterns.
To verify our findings based on the CESM-LME, we also compared them to the observations from the NASA Goddard's Global Surface Temperature Analysis (GIS-TEMP; Hansen et al. 2010) and the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST; Rayner et al. 2003) for the period from 1880 to 2021.A Student's t test was used to determine the significance level at 95% confidence for both regression and composite analyses.The CESM-LME has been found to provide a fairly good simulation of the surface temperature (Otto-Bliesner et al. 2016).Additionally, it simulates ENSO evolution well in terms of seasonality, frequency, and associated teleconnections (Bellenger et al. 2014;Capotondi et al. 2020), as well as the ability to simulate El Niño response to large tropical eruptions (Liu et al. 2018(Liu et al. , 2022;;Stevenson et al. 2016Stevenson et al. , 2017;;Yang et al. 2022).As shown in Fig. S1, despite an overestimation of the intensity of El Niño and La Niña, all 13 ensembles are capable of simulating the seasonal phase-locking of the ENSO evolution in the winter time, consistent with previous works (Stevenson et al. 2016).
The ENSO teleconnections are considered to be the stationary waves that originate from a source region (Hoskins and Karoly 1981;Hoskins and Ambrizzi 1993).Therefore, the Rossby wave theory (Rossby 1945;Yeh 1949), which is based on dispersion and group velocity, is commonly applied to understand the variation of these teleconnections (Hoskins and Karlay 1981;Karlay 1983;Li and Li 2012;Li et al. 2015;Zhao et al. 2015).Thus, the wave ray based on the wave theory linearized on the zonal mean flow (Hoskins and Karoly 1981), was used to depict the propagation of Rossby waves.
A wave ray is a curve that depicts the wave energy dispersing at its group velocity.Its displacements along the longitude X and latitude Y follow the formulas where c g = (c gx , c gy ) denotes group velocity.The group velocity is in a form of whereas k , l are zonal and meridional wavenumbers, and − q y is the meridional gradient of the basic-state absolute vorticity.For a stationary wave, k is the zonally mean zonal wind in the Mercator Projection − u ∕cosϕ , and ϕ is the latitude.For a given flow and zonal wavenumber, we integrated Eq. ( 2) to derive its wave ray with the meridional wavenumber varying equation as following In order to examine the effect of the zonally variation of the zonal wind on the propagation, we used the local quantities of basic state (the bar terms in the above formulas) instead of the zonal-mean quantities (Hoskins and Ambrizzi 1993).Although wave theories based on more complicated flows have been developed (Karoly 1983;Li and Li 2012;Li et al. 2015;Zhao et al. 2015), this simple theory still works well on the analysis of ENSO teleconnections in mid-to-high latitudes of the Northern Hemisphere.To calculate the wave ray, the fourth-order Runger-Kutta method was applied to integrate Eqs. ( 2) and (3).The time interval was set as 1 h and the total integration time is 30 days.

Climatic ENSO effects on global LST
To examine overall impact of ENSO on global LST, we conducted a diagnosis of the climatological regression of global LST onto boreal-winter ENSO, i.e., the average of 229 × 13 regression maps during 850-2005AD (Fig. 1a, c).
In this work, our analyses only involved the regression on the Niño3.4index and did not consider the asymmetry between the El Niño and the La Niña.Therefore, the teleconnection maps we presented are related to an El Niño, while those related to a La Niña would present similar patterns but with the sign reversed.
The teleconnection patterns during the developing summer (MJJAS − 1), boreal winter (NDJFM), and decaying summer (MJJAS + 1) exhibit distinct differences in term of both the sign and intensity.In MJJAS − 1 (Fig. 1a), the impact of El Niño on LST manifests as cooling in middle latitude of Asia, central Asia, East Asia, the Maritime Continent, northern Australia, eastern and southwestern North America, while it triggers warming in Africa, Alaska, South America, western Antarctic near the Amundsen Sea, and Dronning Maud Land in Antarctica.
During NJDFM when the El Niño peaks (Fig. 1b), there is a noticeable intensification of the ENSO-LST teleconnections.An El Niño is associated with the strongest cooling in the northern Eurasian continent, while cooling in the Tibetan Plateau and mid-latitude North America is also strengthened.Other regions in the world exhibit notable warming, especially West Asia, high-latitude North America, South America, Greenland, Australia, and the Antarctic.In the mid-to-high latitudes of the Northern Hemisphere, the significant characteristic of the ENSO-LST teleconnections is the considerable cooling in Eurasia, while warming is observed in northern Canada and Greenland.This phenomenon is known as the warm Arctic-cold Eurasia (WACE) pattern (Cohen et al. 2012;Mori et al. 2014;Overland et al. 2015;Sun et al. 2016;Feng et al. 2021).
During the El Niño decaying summer (Fig. 1c), the El Niño-LST teleconnections are not as robust as those in boreal winter, but they are still stronger than in previous summer.There is still robust warming in South America and West Antarctic near the Amundsen Sea.Additionally, the signs of the regression coefficients change from negative to positive over Europe and India, indicating an opposite teleconnection between the boreal winter and the decaying summer phase of an El Niño.

Decadal variations of ENSO-LST teleconnections
To investigate the decadal-to-multidecadal variations of ENSO-LST teleconnections, we employed a seasonal EOF analysis on the anomalous regression maps.EOF1 explains about 7.7% of the total variance, which can be significantly separated from the other higher modes.The significance for EOF1 is determined by whether the differences between PC1 > + 1.0σ and PC1 < − 1.0σ (where σ is the standard deviation of PC1) are significant at the 95% confidence level from Student's t test.It can be seen that the strongest decadal variability of ENSO-LST teleconnections occurs in boreal winter (Fig. 1e), which is consistent with the strongest LST responses to ENSO in winter (Fig. 1b).For an El Niño, EOF1 mainly exhibits significant cooling in the northern Eurasian continent during boreal winter (Fig. 1e), with the temperature dropping below − 0.4 ℃.On the other hand, the warming of over 0.3 °C can be observed in western North America from Alaska to Idaho and Greenland.While decadal variabilities of other regions that pass the significance test are less than 0.05 ℃.This result suggests that the anomalous pattern of the dominant decadal variation of El Niño-LST teleconnection resembles the climatological structure, i.e., the WACE pattern.
During both MJJAS − 1 (Fig. 1d) and MJJAS + 1 (Fig. 1f), the decadal variability does not show as much strength as that during the winter time.However, it is worth noting that despite the smaller magnitude of changes, there is still a significant cooling trend of more than − 0.2 ℃ in the West Antarctic near the Amundsen Sea during both MJJAS − 1 and MJJAS + 1.The results suggest that the most significant decadal variation of ENSO teleconnections on global LST occurs in the mid-to-high latitudes of the Northern Hemisphere during the boreal winter.During both developing and decaying summers of the ENSO, the decadal variation is notable in the West Antarctic region.Meanwhile, we also conducted a similar analysis using a longer sliding window, i.e., a 31-year sliding window (Fig. S2), and found that it produced similar results to those obtained with the 11-year sliding window.
To identify the decadal non-stationary teleconnections, we composited the two opposite phases of the dominant decadal variation, selected based on the amplitudes of PC1 being greater than + 1.0σ and less than − 1.0σ, respectively.When the amplitude is greater than + 1.0σ (Fig. 2b), El Niño events are mainly associated with the extremely cold temperature in the northern Eurasia that extends to East Asia during boreal winter, while the Greenland experiences extreme warming.On the contrary, when the amplitude is smaller than − 1.0σ, the cold temperature is weakened over the high latitudes of Eurasia, and LST response to the El Niño over South Siberia is even reversed to show warming.At the same time, warming in Greenland is also weakened, and the eastern part of this island even shows cooling (Fig. 2e).Additionally, it is observed that during both developing and decaying summers of an El Niño, the intensity and extent of warming over the Antarctic region for PC1 > + 1.0σ are smaller than for PC1 < − 1.0σ.In the following text, we will use the terms "EA-cold phase" and "EA-warm phase" to refer to the two phases with PC1 > + 1.0σ and PC1 < Each ensemble is a segment of the whole sequence before the EOF analysis.Stippling indicates statistical significance at the 95% confidence level for the difference between PC1 > + 1.0σ and PC1 < − 1.0σ, where σ is the standard deviation of PC1 − 1.0σ, respectively, as the most significant distinctions in ENSO-LST teleconnections are observed in Eurasia during the boreal winter.
Figure 3 exhibits the composite ENSO-related precipitation anomalies during the above two opposite phases.In MJJAS − 1, the differences of precipitation anomalies associated with two different phases are found in the middle and lower reaches of the Yangtze River in East China, with less rainfall in the EA-cold phase and more rainfall in the EA-warm phase (Fig. 3a, d).Meanwhile, the opposite eastwest distribution of precipitation anomalies in Greenland is reversed between different phases.During the boreal winter, northern Eurasia tends to experience less rainfall in the EAcold phase and more rainfall in the EA-warm phase (Fig. 3b,  e), and opposite rainfall responses also occur in the eastern Greenland between these two phases.During the El Niño decaying summer, the differences of precipitation responses to different phases are relatively weak.In western Antarctic near the Amundsen Sea, more rainfall is found in the EAwarm phase compared to that in the EA-cold phase during both developing and decaying summers of the El Niño.In southwest North America, precipitation responses to an El Niño for the EA-cold and warm phases are similar from the El Niño developing summer to decaying summer, indicating that there is a stable wet response to El Niño during the last millennium in CESM.
To understand why there is significant decadal variation in ENSO-LST teleconnections, we examined the composite ENSO-related SST structure and its 200-hPa geopotential height teleconnections for the EA-warm and EA-cold phases (Fig. 4).The SST anomalies exhibit the well-known evolution of El Niño, Indian Ocean dipole (IOD) (Klein et al. 1999;Meyers et al. 2007;Du et al. 2009;Chowdary et al. 2014;Zhang et al. 2015;Stuecker et al. 2017) and North Atlantic Tripole (NAT) anomalies (Sutton and Hodson 2003;Rodríguez-Fonseca et al. 2016; Jiménez-Esteve and Domeisen 2018) from MJJAS − 1 to MJJAS + 1.The most substantial difference between the EA-cold and warm phases is the location of the maximum SST warming in the tropical Pacific during the boreal winter.In the EA-cold phase, the maximum SST warming associated with the El Niño is located closer to the dateline, more like the Central-Pacific (CP) El Niño (Ashok et al. 2007).This westward shift of the maximum SST warming is also observed during the boreal summer.In the EA-warm phase, the maximum SST warming is located in the eastern Pacific, as the Eastern-Pacific (EP) El Niño.4b, e).In the Southern Hemisphere, a PSA-like wave train, represented by southward high + low + high pressure anomalies from the equatorial Pacific to the Antarctic along 120° W, is also observed from the El Niño developing summer (Fig. 4a, d) to decaying summer (Fig. 4c, f).
Notably, atmospheric teleconnections associated with these different types of El Niño in boreal winter are quite different.During the EA-cold phase associated with a CP-like El Niño, the PNA teleconnection exhibits a zonaltype teleconnection pattern characterized by the 200 hPa geopotential height anomaly centers.This pattern shows a "+ − + − +" distribution, spanning from the tropical Pacific, across the North Pacific, to North America, Greenland and the North Atlantic (Fig. 4b).As a result, the Arctic polar vortex is strengthened over Greenland and weakened over northern Eurasia, leading to strong northerly wind and cooling over northern Eurasia.During the EA-warm phase (Fig. 4e), however, the PNA teleconnection exhibits a meridional-type teleconnection pattern, stretching from From MJJAS − 1 to MJJAS + 1, a robust PSA-like wave train pattern, characterized by a low-pressure system over the middle latitudes of the South Pacific, a high-pressure system near the Amundsen Sea of the western Antarctic and a low-pressure system over South America and the South Atlantic, is emanated from the tropical Pacific in both phases.The weakening of the Amundsen Sea low due to western Antarctic high-pressure anomalies results in a substantial warming in the western Antarctic.The strength of the PSA-like wave train varies depending on the type of El Niño.The PSA-like wave train is stronger for an EP-like El Niño (Fig. 4d, f) compared to a CP-like El Niño (Fig. 4a,  c), which leads to a greater weakening of the Amundsen Sea low and more warming in western Antarctic during the EP-like El Niño.
In a linear theory proposed by Hoskins and Karoly (1981) and Liu and Wang (2013), the teleconnections present in the upper-troposphere are primarily determined by the background wind and Rossby wave source that originates from the associated convection.By utilizing the wave ray theory, we sought to understand why the westward shift of CP-like El Niño convection can produce a zonal-type teleconnection, whereas an EP-like El Niño favors a meridional one.In the CESM-LME, background zonal wind differences between EA-cold and warm phases are quite small (Fig. 6).Under different background flows between two phases, the wave rays starting from the same location (at 162.5°W or Different wave rays (thick black curve) from starting points (black dot) to the west of Hawaii for a EA-cold phase (PC1 > + 1.0σ) and to the east of Hawaii for b EA-warm phase (PC1 < − 1.0σ), respectively.Distributions of 200-hPa geopotential height anomalies (shading; m) associated with ENSO and background zonal wind (black dashed contours of 20, 30, 40, 50, 60 m/s) are also shown 145° W) do not reveal any significant differences in the trajectory (not shown).Under the same background flow (average of two phases), the wave ray originating from 162.5° W in EA-cold phase shows different propagation trajectory with that from 145° W in EA-warm phase, although they both exhibit a trajectory resembling a "great circle" (Fig. 6).In EA-cold phase, the ray reaches the Greenland in a more zonally propagating pattern, whereas in EA-warm phase, it displays a more meridionally propagating pattern towards the Arctic, consistent with the CESM simulations (Fig. 4b,  e).This result implies that the shift of the wave source location, rather than the background flow differences, determines the distinct ENSO teleconnections between the EA-cold and warm phases.

IPO-determined decadal variation in ENSO teleconnections
To understand the factors contributing to the decadal variation of ENSO teleconnections in the CESM-LME, a regression analysis of global 11-year smoothing SST anomalies against PC1 was conducted (Fig. 7a). Figure 7a exhibits a typical negative phase of IPO, which is characterized by significant negative SST anomalies in the central and eastern equatorial Pacific and positive anomalies in the western Pacific.Along with this simulated negative IPO, the North Atlantic displays negative SST anomalies.The SST patterns associated with the IPO usually persist for around 20-30 years.This is evidenced in the spectrum analysis of PC1, which indicates prominent peaks of oscillation at roughly 20 years, 30 and 34 years (Fig. 7b).ENSO teleconnections are usually triggered by its associated convection anomaly over the tropical Pacific.During the positive phase of IPO, the EP-like El Niño, which is associated with the EA-warm phase, tends to trigger a positive convection anomaly that remains constrained to the eastern tropical Pacific.However, during the negative phase of IPO and EA-cold phase, the positive convection anomaly associated with CP-like El Niño moves westward (Fig. 5).During positive and negative IPO phases, the change of the El Niño convection location contributes to different ENSO teleconnections.
Decadal climate variability can arise from both internal variability and external forcing.Internal variability is often regards as the "noise", while the external forcing due to solar variability, volcanic activity, and greenhouse gas emission is generally considered as the "signal".To investigate the extent which internal variability and external forcing contribute to the decadal variations of ENSO teleconnections in CESM-LME, three approaches were employed in this study.
We first employed the average of 13 PC1s to filter out the effects of internal variability (Fig. 8a).The remaining signals can be attributed to external forcing (Zanchettin et al. 2013;Si and Hu 2017;Jia et al. 2020).Throughout the last millennium from 850 to 2005AD, the global temperature underwent several notable periods, including the Medieval Warm period around 1000AD, the cold Little Ice Age from the 14th century to 19th century, and the current Modern Warm Period (Chai et al. 2020).However, the average of  8a).Notably, these latter two periods of negative signals are not closely associated with the wellknown large eruptions of 1695 and 1982AD, respectively (Liu et al. 2022).During 850 to 2005AD, there are over 40 large tropical volcanic eruptions, yet the ensemble average only displays significant signals three times.This suggests that the small ensemble size of 13 may be responsible for the significance, and that averaging only 13 ensembles may not have entirely removed the interannual variability.To further check whether the long-term forcing from solar radiation and anthropogenic greenhouse gas concentration are responsible for any observed signals, we performed the same EOF analysis on the 50-and 100-year-sliding regression.Despite this, no signals can be simulated in the 13-ensemble average (Fig. 8b, c).
The second approach we used to quantitatively estimate the relative contributions of internal variability and external forcing is the analysis of variance (ANOVA) method (Harzallah and Sadourny 1995;Ting et al. 2009).The I 2 and F 2 represent the biased estimates of the internal and external forced variances, respectively.The r is the ratio of forced variance to total variance.For a relatively large ensemble sizes, the F 2 should be well approximated by the ensemble averaged variance.The small and negative F 2 and r in Fig. 8 indicates that the decadal variation of ENSO-LST teleconnections may be due to the dominant role of internal variability, although a relatively small ensemble size of 13 could also be a reason.
The third approach involves employing the EOF analysis on the averaged regression maps of these 13 ensembles, and the leading EOF mode of the ensemble mean presents the dominant forced signal (Allen and Smith 1997;Chang et al. 2000;Ting et al. 2009;Li and Ting 2015;Jia et al. 2020).It can be seen that the decadal variability of the forced signal (Fig. S3) is much smaller than that shown in Fig. 1.This finding further illustrates that the internal variability plays a more important role on the decadal variability than external forcing.
Based on these analyses, it can be concluded that the dominant decadal variation of ENSO teleconnections in CESM-LME is mainly attributed to internal variability rather than to external natural or anthropogenic forcing.

Observed decadal variations of ENSO teleconnections
To explore whether the simulated decadal variation of ENSO teleconnections exists in the observations, we performed the same analyses on the GISTEMP and HadISST spanning the period 1880-2021.As shown in Fig. 9a, c, the climatological ENSO-LST teleconnections in the observations are similar to those simulated in CESM-LME (Fig. 1a, c).The observation also displays a WACE pattern, i.e., significant winter cooling in Eurasia and warming in Arctic over northern Canada and Greenland, related to an El Niño.The leading EOF of 11-year-sliding regression in the observations, which accounts for about 19.7% of the total variance, also exhibits prominent wintertime cooling over northern Eurasia, consistent with the simulations, although North America displays an opposite pattern compared to the simulations (Fig. 9e).Related to this observed decadal variation of ENSO-LST teleconnections, the cold SST anomalies are also observed in the eastern equatorial pacific, although not significant, whereas significant warming is observed in the subtropic western Pacific (Fig. 10).This anomalous SST gradient mimics that of IPO simulated by CESM-LME.The tropical North Atlantic, however, exhibits significant warming related to cold eastern equatorial Pacific in the observations, against the cooling in the simulations.

Summary and discussion
In this study, we investigated decadal variations of ENSO-LST teleconnections, from the developing summer to decaying summer of ENSO.During the peak phase of an El Niño, i.e., boreal winter, CESM-LME shows substantial climatic teleconnections, which results in significant cooling in Eurasia along with warming Arctic over northern Canada and Greenland, known as the WACE pattern.It also displays warming in the western Antarctic during both developing and decaying summers of an El Niño.The dominant decadal variation of ENSO-LST teleconnections exhibits significant perturbations in the mid-to-high latitudes of the Northern Hemisphere during the boreal winter, with a variability pattern similar to the WACE.Additionally, the decadal variation is also notable in the western Antarctic region during developing and decaying summers of ENSO.
The Eurasian winter LST response to an El Niño displays opposite variations during the positive and negative phases of the decadal variation, referred to as the "EA-cold phase" and "EA-warm phase".These different responses are attributed to the distinct patterns of anomalous Arctic polar vortex, forced by the diversity of El Niño.During the EA-cold phase, the El Niño exhibits a CP type and its Rossby wave source excited by precipitation is situated to the west of Hawaii.This Rossby wave forcing tends to excite a zonal-type teleconnection, which weakens the Arctic polar vortex near Eurasia, results in a cold and dry winter over the high latitudes of Eurasia and a warm and dry winter in the high latitudes of North America.During the EA-warm phase, the El Niño exhibits an EP type with its Rossby wave source located to the east of Hawaii.As a result, it triggers a meridional-type teleconnection to shift the Arctic polar vortex toward northern Eurasia.This leads to a warmer and wetter winter over Eurasia compared to the EA-cold phase.We further discovered that the decadal variations of ENSO teleconnections during the last millennium are mainly attributed to the internal IPO mode rather than external forcing in CESM-LME.This conclusion is based on a relatively small ensemble size of 13, and a large ensemble size should be considered to support this finding in future work.
The stability of ENSO teleconnections is a critical concern when using proxy records, such as coral, tree rings, and sediment cores, to reconstruct paleo-ENSO.Several studies have demonstrated the decadal-to-centennial variability of the ENSO (Diaz et al. 2001;Hope et al. 2017), and the non-stationarity relationships between ENSO and remote climates in specific oceans (Lewis and LeGrande 2015;Christopher 2018;Park and Li 2019;Liu et al. 2021), land areas (López-Parages et al. 2012, 2016;Coats et al. 2013), and even polar regions during the last millennium (Feng et al. 2021;Li et al. 2021).As a supplement and extension of previous studies, we investigated the nonstationary behaviors of atmospheric teleconnections associated with the ENSO during a millennium period by using the outputs of the CESM-LME.The significant perturbation simulated in the WACE pattern highlights the importance of considering the decadal variability when reconstructing ENSO using proxies from these regions.It is worth noting that the relatively stable precipitation response over southwestern North America, in the presence of this decadal variability, provides further confidence in the reliability of using tree rings from this region to reconstruct the ENSO variability (Li et al. 2011).

Fig. 1
Fig. 1 ENSO teleconnections and its decadal change.Climatological regression of global land surface temperature (LST; °C) during a developing summer (MJJAS − 1), b boreal winter (NDJFM), and c decaying summer (MJJSA + 1) with respect to boreal winter Niño3.4 index for 13 CESM-LME ensembles from 850 to 2005AD.Stippling indicates regression relationship significant at the 95% confidence level from Student's t test.d-f Show the first mode of three-

Fig. 2
Fig. 2 Different global LTS responses to ENSO.Composite LTS regression (°C) in a, d developing summer, b, e boreal winter, and c, f decaying summer with respect to ENSO when PC1 > + 1.0σ (left

Fig. 3
Fig. 3 Different global precipitation responses to ENSO.Same as Fig. 2, except for precipitation (mm) over land

Fig. 4
Fig. 4 Different global SST and circulation responses to ENSO.Same as Fig. 2, except for SST (shading; °C) and 200-hPa geopotential height (Contour; m).The zonal mean component at each latitude Figure 5 presents the boreal winter zonal structures of tropical Pacific (5° S-5° N) SST anomalies and 200-hPa geopotential height anomalies along 20 °N related to an El Niño during the EA-cold and warm phases.The wave ray source, i.e., the center of geopotential height anomalies, originates at 162.5° W to the west of Hawaii in the EA-cold phase, related to a CP-like El Niño.The wave ray source, however, is shifted eastward and starts at 145° W to the east of Hawaii in the EA-warm phase, linked to an EP-like El Niño.

Fig. 5
Fig. 5 Different zonal structures in Pacific between EA-cold and warm phases.The average tropical Pacific (5° S-5° N) SST anomalies (red lines) and 200-hPa geopotential height anomalies (blue lines) along 20° N for EA-cold phase (solid lines; PC1 > + 1.0σ) and EAwarm phase (dashed lines; PC1 < − 1.0σ), respectively.The black dots denote the locations of wave source in the two phases

Fig. 8
Fig. 8 Internal variability determined changes of ENSO teleconnections.Plots of 13-ensmeble mean of PC1 for a 11-year-sliding, b 50-year-sliding, and c 100-year-sliding regression maps with an interval of five years.Shading shows the 95% confidence bounds based on two-tailed Student's t test.The I 2 and F 2 represent the biased estimates of the internal and external forced variances, respectively.The r is the ratio of forced variance to total variance