Revisiting ENSO teleconnections in the boreal winter

1 Introduction

The El Niño-Southern Oscillation (ENSO) has been of significant interest to researchers due to its dominant influence on regional climate variability across the globe (Rasmusson and Carpenter, 1982; Ropelewski and Halpert, 1987; Trenberth et al., 1998). Despite its prominence, ENSO exhibits considerable complexity that introduces uncertainty in both its evolution and impacts. Its irregular periodicity (2–7 years), asymmetry between El Niño and La Niña in amplitude, duration, and spatial structure, and sensitivity to stochastic atmospheric variability all contribute to diversity in ENSO behavior (Cane et al., 1986; Capatondi et al., 2015; Chen et al., 2019a,b; Lau, 1985; Misra et al., 2007, 2008; Newman et al., 2011; Okumura and Deser, 2010; Penland and Sardeshmukh, 1995; Song et al., 2026; Thompson and Battisti, 2001; Tziperman et al., 1997; Wei et al., 2026). In addition, low-frequency climate modes such as the Pacific Decadal Oscillation, Interdecadal Pacific Oscillation, and Atlantic Multidecadal Oscillation modulate the background state of the ocean–atmosphere system, influencing ENSO characteristics and the strength and structure of its teleconnections (Cai et al., 2010; Dieppois et al., 2021; Imada and Kimoto, 2009; Sun and Okumura, 2020; Verdon and Franks, 2006; Wang et al., 2014). More recently, Deser et al. (2026) indicated that SST trends in the tropical oceans generate teleconnections that could counteract and weaken ENSO teleconnections. These interactions across timescales complicate attribution and reduce the robustness of ENSO-related climate signals.

Many studies have evaluated the ability of climate models to reproduce observed ENSO teleconnections (Joseph and Nigam, 2006; King et al., 2015; Sutton et al., 2024; Weare, 2013; Williams et al., 2023). Implicit in such evaluations is the assumption that the observational and reanalysis datasets used as benchmarks provide a robust and internally consistent characterization of ENSO teleconnections. However, few studies have quantitatively addressed how diagnosed teleconnection patterns themselves depend on the choice of observational dataset, ENSO index definition, temporal record length, and temporal filtering methods. Consequently, uncertainties associated with the observational benchmark remain poorly characterized.

This study addresses this gap through a systematic uncertainty assessment of ENSO teleconnections. Rather than treating observational estimates as a single reference, we examine ENSO teleconnections as a diagnostic quantity whose robustness depends on multiple methodological choices. The principal contribution of this study is not simply a comparison of datasets or indices individually, but the development of an integrated diagnostic framework for evaluating the robustness of ENSO teleconnections across observational products and methodological choices. By identifying the relative contributions of dataset selection, ENSO index definition, temporal coverage, and scale separation to teleconnection uncertainty, this work provides a more rigorous basis for interpreting observed ENSO impacts and for evaluating climate model simulations against observational benchmarks.

ENSO teleconnections in the higher latitudes are sensitive to the strength and location of tropical convection (or heating) anomalies, which often results in large diversity of teleconnections (Hoerling et al., 1997; Hoerling and Kumar, 1997). These teleconnection patterns can be understood as stationary Rossby Waves excited by diabatic ENSO-related heating anomalies in the equatorial Pacific, influencing regions such as North America and the Pacific-South American pattern (Hoskins and Karoly, 1981; Horel and Wallace, 1981; Mo and Ghil, 1987; Mo and Higgins, 1998). Modeling studies demonstrate that both barotropic and baroclinic processes are important, with zonal jet-sourced vertical wind shear modulating the amplification and structure of the response (Branstator, 1983; Lee et al., 2009). As a result, they find that ENSO teleconnections are enhanced in the winter hemisphere when the subtropical jet is close to the anomalous tropical heat source in the tropical Pacific. In other regions, such as the Indian monsoon region or the tropical Atlantic, ENSO teleconnections are manifested by the zonal shift in the Walker Circulation that modulates the large-scale atmospheric stability and vertical wind shear (Krishnamurthy and Goswami, 2000; Webster et al., 1998). Teleconnections can also produce extreme regional impacts, such as over East Asia (Chen et al., 2023), while weak responses in some extratropical regions may reflect compensating non-linear and linear atmospheric processes (Frauen et al., 2014). Thus, it may be a reasonable assumption that different reanalyses and observational datasets have variations in how they model these mechanisms, leading to uncertainties in the teleconnections.

The ENSO teleconnections discussed in the paper will be shown via contemporaneous and lagged correlations between the forcing and the response regions (Nigam and Baxter, 2014). While such correlation maps are not discriminatory between the forcing and the response regions, their physical interpretation is guided by the well-established framework of ENSO-forced atmospheric responses (e.g., Branstator, 1983; Lee et al., 2009; Mo and Ghil, 1987; Ropelewski and Halpert, 1987). This established understanding motivates the application of a more rigorous statistical test that is under-utilized in evaluating ENSO teleconnections to assess the robustness of these canonical teleconnection patterns. The objectives of this study therefore are to investigate the following:

• Are upper-air and surface ENSO teleconnection patterns consistent across a variety of atmospheric reanalyses and surface observations?
• Are ENSO teleconnections in the datasets sensitive to the choice of ENSO index?
• To what extent do other temporal scales besides interannual scales influence ENSO teleconnections?
• To what extent do the parametric choices in the chosen methodology affect the diagnosis of ENSO teleconnection?
• Which of the above factors contributes most to the uncertainty of ENSO teleconnections?

Unlike previous studies that primarily assess ENSO teleconnections within individual datasets or evaluate model fidelity against a single observational reference, this study quantifies the observational uncertainty of ENSO teleconnections themselves. By jointly evaluating sensitivities to dataset choice, ENSO index definition, record length, and temporal scale separation, we identify the dominant sources of uncertainty in diagnosed teleconnection patterns and assess the robustness of commonly accepted ENSO impacts. The results provide guidance for both observational analyses and climate-model validation studies that rely on teleconnection diagnostics.

2 Datasets

Several indices are used to monitor ENSO signals. The ENSO indices used in this study are the ocean-based Oceanic Niño Index (ONI), Multi-variate ENSO Index (MEI), and atmospheric-based Southern Oscillation Index (SOI), which are summarized in Table 1.

Table 1. . Summary of ENSO indices used in this study.

The ONI is calculated using SST anomalies in the east central equatorial Pacific, also known as the Niño-3.4 region (5°N−5°S, 170–120°W; Barnston et al., 1997). It is reported as a 3-month running mean using ERSSTv5 SST data (Huang et al., 2017). Studies have shown that SSTs in the Niño-3.4 region have been found to statistically capture the state of ENSO and various ENSO teleconnections (Giese and Ray, 2011; Hanley et al., 2003; Trenberth et al., 1998).

The Multivariate ENSO Index (MEI) measures ENSO signal as expressed by the first principal component of six observed meteorological variables over the tropical Pacific: sea level pressure (SLP), surface zonal wind component, surface meridional wind component, SSTs, surface air temperature, and cloudiness (Wolter and Timlin, 1993, 1998). It is re-computed every month from the tropical Pacific Comprehensive Ocean-Atmosphere Data Set (COADS) records to monitor the strength of ENSO conditions for the preceding 2 months. Previous work has found that MEI correlates well with pressure-based and SST-based indices in terms of identifying ENSO phases (Hanley et al., 2003). However, its strongest benefit is its ability to capture a more complete picture of the atmospheric and oceanic anomalies caused by ENSO variations. A limitation of the MEI is its reliance on homogeneous surface wind data, which is unavailable prior to 1950.

The Southern Oscillation Index (SOI) is a pressure-based index based on the observed SLP differences between Tahiti, French Polynesia and Darwin, Australia (Ropelewski and Halpert, 1987; Trenberth, 1984; Walker, 1932, 1937). Unlike the previously discussed indices, the SOI is inversely related to ENSO teleconnections, with a negative (positive) phase of the SOI corresponding to an El Nino (La Nina) event. The SOI tends to exhibit greater short-term variability (i.e., greater noise) than SST-based indices (Hanley et al., 2003; Trenberth, 1976, 1984), which can hinder the clear identification and correlation of ENSO teleconnections.

The other datasets used in this study are listed and summarized in Table 2. The oldest, and spatially coarsest datasets amongst the reanalyses used in this study is the NCEP-DOE AMIP-II Reanalysis (NCEP-R2) which has a grid size of 2.5° by 2.5° and is available from 1979 to 2018 (Kanamitsu et al., 2002). While it incorporates improved physical parameterizations relative to its predecessor, the NCEP-NCAR Reanalysis (NCEP-R1; Kalnay et al., 1996), it lacks more recent advancements such as finer spatial resolutions and the assimilation of radiance, SSM/I, and rainfall data. JRA55 is more contemporary (Japan Meteorological Agency/Japan, 2013) and is the longest third-generation reanalysis to use a full observing system with a 4D-var assimilation system. It extends back to 1958 and has a finer spatial resolution (1.25° x 1.25°) than NCEP-R2. Contemporary datasets have an improved representation of the atmosphere through finer spatial grid scales, improved data assimilation, and better model parameterization. ERA5 and MERRA2 incorporate a larger volume of observations (Gelaro et al., 2017; Hersbach et al., 2020) and operate at spatial resolutions finer than 1°. Although MERRA-2 only slightly exceeds NCEP-R2 temporally, ERA5 extends much further back, reaching 1940. Despite these shared advancements, they do still hold some differences in their assimilation approaches. Most notably, ERA5 utilizes a 4D-var data assimilation system, whereas MERRA2 employs a 3D-var system. This distinction may add uncertainty in teleconnections between datasets, as 4D-var systems incorporate the temporal evolution of the forecast model, which allows for an improved representation of the time distribution of observations (Rabier and Liu, 2003).

Table 2. A summary of datasets used in the study.

The Climate Hazards group Infrared Precipitation with Station data version 2 (CHIRPSv2) and Climate Research Unit (CRU) precipitation datasets are surface observational products (Funk et al., 2015; Harris et al., 2020; Xie et al., 2007). CHIRPSv2 integrates station measurements and satellite data, using the latter to fill gaps in regions with insufficient station coverage. CRU, however, relies solely on weather station records and uses a climatological value in grid points with insufficient reporting stations. Both datasets have spatial resolutions finer than 1° and are thus regridded to a common 1° resolution as is done with precipitation re-analyses from ERA5 and MERRA2 in this study.

For land surface temperature, CRU was the sole observational dataset used in this study (University of East Anglia Climatic Research Unit et al., 2023). Very few multi-decadal datasets of land surface temperature exist, primarily because satellite measurements typically represent skin temperature instead of the 2-m air temperature. However, unlike rainfall data, the spatial and temporal gradients of surface temperature are less heterogenous and hence far less uncertain. But the density of stations is highly variable, which makes construction of global temperature analysis a challenge. Therefore, land surface temperature analysis in this study relies on a single source of land-based meteorological station data.

We make use of the full extent of the overlapping dataset between the ENSO index (Table 1) and the upper-air/surface variable (Table 2) in our analysis as indicated in last column in Table 2. However, while intercomparing reanalyses (ERA5, MERRA-2, JRA55, and NCEP-R2), we also use a limited period of 39 years (1980–2018) when all reanalyses datasets in Table 2 overlap to make a fair comparison. For the comparison of precipitation observational datasets (CHIRPSv2 and CRU), this overlapping period is 42 years (1981–2022). The datasets in Table 2 used for the diagnosis for ENSO teleconnections were linearly interpolated to a common 1° grid spacing for fair comparison across datasets. This was also done to reduce the computational burden of temporal filtering of these datasets (explained in the next section). Furthermore, since we are examining the large-scale ENSO teleconnection patterns, we find that degrading the resolution of some of the datasets does not affect the results of this study.

3 Methodology

3.1 Empirical ensemble mode decomposition

We make use of Empirical Ensemble Mode Decomposition (EEMD; Wu and Huang, 2009) to conduct a temporal analysis by separating the scales in the various variables analyzed in this study to synoptic, subseasonal, interannual, decadal, and secular timescales (Table 3). EEMD is a data-adaptive technique for analyzing time series without relying on predefined basis functions such as the sine and cosine waves used in Fourier methods. Instead, it extracts intrinsic oscillatory modes directly from the data based on local time-scale separation. EEMD is extensively used in climate studies (e.g., Feng et al., 2014; Franzke, 2009; Huang et al., 2012a,b; Misra et al., 2012; Wu et al., 2009).

Table 3. Temporal thresholds used for scale separation.

EEMD is a noise-assisted data analysis technique. It is an enhancement of the original Empirical Mode Decomposition (EMD), which is well-suited for analyzing non-linear and non-stationary signals (Huang et al., 1998). EEMD “robustly” isolates the oscillations in the datasets using the EMD approach by introducing ensembles of white noise realizations into the decomposition process. Throughout the trials, the white noise is averaged out, leaving behind the decomposed oscillatory modes or the Intrinsic Mode Functions (IMFs).

The specific parameters used for EEMD in this paper are stated in Table 4. The amplitude of the noise (ϵ) added to signal is 0.2 × standard deviation of the timeseries following Wu and Huang (2009). The number of sifting iterations is 10 per IMF, and the tolerance level to end the sifting process is ϵN, where N is the number of ensembles (Wu and Huang, 2009). In the process of sifting, EEMD isolates the secular time scale as the residual (last) IMF, wherein, no further iterations are likely to result in oscillatory modes. In this study, we find that this residual mode appeared as a secular variation only in the surface temperature field, while the other variables (e.g., geopotential height fields, precipitation) displayed the last IMF with an incomplete oscillatory mode with a single sign change over the period of the dataset. Hence, these incomplete oscillatory modes were not regarded as secular variation.

Table 4. EOF-EEMD defined parameters.

The datasets (Y(t, s)), used in this study vary in space (s) and time (t). Instead of conducting EEMD analysis at every grid point, which can be computationally intensive, we first decompose the datasets using EOF analysis. Such datasets can be decomposed into m Principal Components (PCs) of Empirical orthogonal functions (EOFs) as:

Y(t,s)=∑mKem(s)pm(t)(1)

where e_*m* and p_*m*are EOFs and PCs, respectively. In this study, EOF analysis is first performed, and the first 10 (K = 10) principal components (PCs) are retained for subsequent EEMD decomposition. The sensitivity of the results to the choice of number of retained PCs is discussed in the results. EEMD decomposes the time series of the PCs into complete sets of near-orthogonal components called Intrinsic Mode Functions (IMFs). These IMFs of the PCs are then used to obtain the anomalies (Y′(t, s)) for each frequency as:

Y′(t,s)=∑m=1Kem(s)IMFm(t)(2)

The reconstructed components, Y′(t, s), were then grouped based on their frequency following Table 3. It may be noted that Misra and Jayasankar (2024) find that truncating the EOFs with 10 PCs is a conservative choice, since even for a discontinuous variable like precipitation, the percentage of variance explained by individual PCs begin to asymptote after the first few PCs. We will show the sensitivity of our results to this truncation of PC later in this study. Owing to the substantial computational cost of applying the EOF–EEMD procedure to global datasets at 1° spatial resolution, the temporal decomposition of ENSO teleconnections in this study is restricted to ERA5. Therefore, the purpose of this analysis is to understand the impact of the various temporal scales of variations on the interannual ENSO teleconnections in ERA5.

3.2 Benjamini-Hochberg test

Since most meteorological variables are spatially correlated, a field significance test is necessary to determine whether a group of spatially distributed hypothesis (like grid point correlations passing a two tailed t-test) show a significant overall signal after accounting for spatial dependence among the individual tests with each grid point (Livezey and Chen, 1983; Ventura et al., 2004; Wilks, 2020). The spatial correlation of the variables itself inflates the significance of the correlations with ENSO indices (i.e., the teleconnection) at the individual grid points. Therefore, testing the significance of the correlation at each grid point would thus increase the “False Discovery Rate” (FDR), which is the expected proportion of false positives among all significant results when conducting multiple comparisons.

In this study, we utilize the Benjamini-Hochberg (BH) test for field significance (Benjamini and Hochberg, 1995). The BH test adjusts the calculated p-value at all grid point locations to control the FDR, thereby minimizing the likelihood of incorrectly rejecting the null hypothesis across multiple comparisons. By managing the FDR, the reliability of the significance of the results is enhanced. A grid point is deemed significant if its calculated p-value is equal to or less than the critical p-value, determined by the following formula:

Pk≤kmα(3)

In this equation, k represents the rank order of the p-value at a specific grid point and m denotes the total number of grid points. The significance level (α ) used in this study was 0.05. It may be noted that when a temporal correlation is considered robust or significant hereafter, then it refers to BH test being passed with 5% significance level.

4 Results

4.1 A comparison of ENSO indices

Figure 1a shows the cross-correlations between the three ENSO indices, which indicates that most indices are very strongly correlated to each other except the SOI. The unfiltered SOI time series shows the lowest (and statistically insignificant) correlations with the rest of the ENSO indices. However, in Figure 1b all the ENSO indices including the SOI are filtered (using EEMD) and only the interannual scales are retained. We now find that in Figure 1b, most of the correlations resemble Figure 1a except for the correlations of the other ENSO indices with the SOI, which now displays uniformly a statistically significant negative correlation. As noted earlier, the high contribution of noise in the SOI timeseries is a significant issue (Hanley et al., 2003; Trenberth, 1976, 1984). But by filtering this noise, the SOI timeseries could be made comparable to other ENSO indices (Figure 1b). For consistency, we will use the interannually filtered values of all ENSO indices for examining the teleconnections and refer hereafter to these as filtered ENSO indices.

4.2 ENSO teleconnections with 200 hPa geopotential height field

The filtered ONI- and MEI-based correlations with 200 hPa geopotential heights exhibit remarkably consistent large-scale teleconnection patterns across all reanalyses, characterized by widespread positive correlations in the tropical Pacific and a prominent positive height anomaly over the North Pacific–North American sector (Figure 2). Significant correlations are concentrated over the tropical Pacific, the North Pacific, and parts of the Southern Ocean, indicating a robust ENSO teleconnection (Figure 2). The MEI generally produces slightly stronger and more extensive correlation patterns than ONI.

In contrast, correlations based on the filtered SOI display weaker amplitudes, reduced spatial coherence, and fewer statistically significant regions (Figure 2). While some features of the canonical ENSO teleconnection remain evident, particularly over the North Pacific and portions of the Southern Hemisphere, the spatial structure varies more among the reanalyses. Differences are especially apparent at mid- and high latitudes, where the magnitude and extent of significant correlations differ among datasets (Figures 2c, f, i, and l). The corresponding lagged correlations with the ENSO indices leading by 1 month (Supplementary Figure S1) show similar patterns as in Figure 2.

4.3 ENSO teleconnections with 500 hPa geopotential height field

Figure 3 presents the correlations between DJF-mean 500-hPa geopotential height anomalies from ERA5, JRA55, MERRA2, and NCEP-R2 and the three filtered ENSO indices. Compared with the corresponding 200-hPa correlations shown in Figure 2, the 500-hPa teleconnection patterns are generally weaker, less spatially coherent, and considerably more zonally asymmetric in the tropics. Significant correlations are primarily confined to the Pacific–North American sector and portions of the Southern Hemisphere, with much less extensive coverage than at upper-tropospheric levels. These results are consistent with previous studies showing that ENSO teleconnections are strongest and most zonally symmetric in the upper troposphere, where weaker damping and more rapid propagation of equatorial wave disturbances facilitate a broader atmospheric response (Barnett, 1985; Branstator, 1983; Dima et al., 2005; Newell and Wu, 1992; Wallace et al., 1998; Wu et al., 2000; Yulaeva and Wallace, 1994).

Figure 3 also reveals substantially greater sensitivity of the teleconnection patterns to the choice of reanalysis dataset than was evident at 200 hPa. While the broad structure of the ENSO-related circulation anomalies is generally preserved, notable differences emerge in the amplitude, spatial extent, and statistical significance of the correlations, particularly over the North Pacific, North America, and high-latitude regions. In contrast, the differences among the three ENSO indices are relatively modest within a given reanalysis. Correlation patterns obtained using filtered ONI and MEI are highly similar, while those based on the filtered SOI exhibit somewhat weaker amplitudes but largely retain the same large-scale features. It may be noted that the corresponding lagged correlations with the ENSO indices leading by 1 month (Supplementary Figure S2) show nearly the same structure as in Figure 3.

4.4 ENSO teleconnections with precipitation

The precipitation teleconnection patterns in Figures 4 and 5 exhibit a high degree of spatial complexity and regional variability compared with the geopotential height variations shown in Figures 2 and 3. Significant ENSO-related precipitation anomalies are concentrated primarily over the tropical Pacific, the Maritime Continent, Australia, southern Africa, and portions of North and South America, reflecting the strong influence of ENSO on tropical convection and associated regional rainfall variability. The overall structure of the teleconnections is generally consistent with canonical ENSO impacts, including enhanced precipitation over the central and eastern equatorial Pacific during El Niño conditions and reduced rainfall over the western Pacific and Maritime Continent. Figure 4 also demonstrates substantial agreement between CHIRPS2 and CRU in the location of the major precipitation teleconnection centers, although differences are evident in the magnitude and spatial extent of the correlations, particularly over data-sparse regions such as northern Africa and central Asia. Similarly, the one season lagged correlations (with ENSO index leading by a season) show similar patterns of correlation as the contemporaneous correlations (Figure 4).

The reanalysis-based precipitation correlations in Figure 5 display considerably larger differences than those found among the observational datasets. While all four reanalyses reproduce the principal ENSO-related precipitation anomalies as noted in Figure 4, notable discrepancies occur in the strength, spatial coherence, and statistical significance of the teleconnections, especially over tropical Africa, and the extratropics. In contrast, differences among the ENSO indices are generally smaller within a given dataset. Filtered ONI and MEI produce nearly identical patterns, whereas the filtered SOI yields slightly weaker correlations but preserves the primary teleconnection features. Overall, the results indicate that uncertainties associated with precipitation datasets and reanalysis products exceed those arising from the choice of ENSO index, highlighting the greater sensitivity of precipitation teleconnections to observational and model-related differences than to the metric used to characterize ENSO.

4.5 ENSO teleconnections with surface temperature

The surface temperature teleconnections identified in CRU in Supplementary Figure S3 exhibit weaker correlations than the precipitation and mid-tropospheric circulation teleconnections. In Supplementary Figure S3, all three filtered ENSO indices display nearly identical patterns in their correlations, characterized by positive (negative) correlations over tropical South America, Africa, and northern Australia and negative (positive) correlations over northeast Asia for ONI and MEI (SOI), albeit these correlations are largely statistically insignificant Similarly, the one season lagged correlations (with ENSO index leading by a season) show similar patterns of correlation (Supplementary Figure S3d-f) as the contemporaneous correlations in Supplementary Figure S3a-c.

Supplementary Figure S4 demonstrates that these ENSO-related temperature teleconnections are inconsistent across the four reanalysis datasets. For instance, the positive correlations over South America, south Africa that appear in correlations with filtered ONI and MEI in ERA5 are not observed in other reanalyses. The correlations over North America however appear to be consistent across reanalyses and ENSO indices. The correlation patterns displayed by ERA5 match the closest with the correlations shown in Supplementary Figure S3, suggesting that ERA5 displays more fidelity in this ENSO teleconnection than other reanalyses. Overall, surface temperature displays weak, statistically insignificant correlations with ENSO index regardless of the choice of the dataset or ENSO index.

4.6 Impact of the discrete temporal scales on ENSO teleconnection

In this analysis, we filter out the various scales of temporal variations (synoptic, sub-seasonal, interannual, decadal, and secular) in the 200 hPa and 500 hPa heights from ERA5, surface precipitation from CHIRPSv2, and surface temperature from CRU using EOF-EEMD to assess their impact on the ENSO teleconnections.

In retaining only the interannual variations in the 200 hPa geopotential height field from ERA5 reanalysis, the correlations are made more similar across all three filtered ENSO indices with subtle differences (Supplementary Figure S5a-c). For example, all indices originally provide a uniform, statistically significant correlations in the global tropics in ERA5 (Figures 2a–c). Figures 6a–c indicates that the area of statistically significant correlations is further enhanced globally when 200 hPa is filtered for interannual scales, with the largest enhancement happening in the global tropics (Figure 6j–l) followed by that in the Northern Hemisphere mid-latitudes (Figure 6g–i). The removal of the synoptic and sub-seasonal scales had minimal impact on the correlations of the 200hPa height variations with the ENSO indices (Figures 6 and Supplementary Figure S5). In contrast, removal of the decadal variations in the 200 hPa height variations of ERA5 enhanced correlations in a comparable manner to the isolation of interannual variations, with most of the enhancement stemming from the tropics.

The influence of the different discrete temporal scales on the correlations of the 500 hPa height variations of ERA5 with the three filtered ENSO indices is far weaker and more non-uniform across ENSO indices (Supplementary Figures S6 and 7). For example, the correlation with ONI shows the most prominent changes in the area covered by statistically significant correlations, when only the interannual scales are retained while the correlations with MEI and SOI indicate some changes in the Southern Hemisphere mid-latitudes but globally, they are minimal.

The various discrete temporal scales of variations of precipitation from CHIRPS2 on its correlations with the filtered ENSO indices appears to be impactful (Supplementary Figures S7 and 8). We find that retaining only the interannual scales of variation in precipitation reduces the area of significant correlation globally, with the global tropical and mid-latitude areas contributing to this reduction. It is observed that the removal of the synoptic, sub-seasonal and decadal variations separately also cause comparable reduction in the areas of significant correlation with the ENSO indices (Figure 8). This feature is uniformly observed across all three three ENSO indices.

The areal coverage of the correlations of the CRU surface temperature with the filtered ENSO indices (viz., ONI and SOI) uniformly get enhanced (and become statistically significant) when only the interannual scales are retained (Supplementary Figures S8 and 9). Interestingly, the correlations with the filtered MEI index show the least sensitivity to separation of the temporal scales (Figure 9) for reasons we are currently uncertain. Furthermore, Figure 9 suggests that this mild global enhancement of the areal coverage of the correlations with filtered ONI and SOI appear largely in the northern polar latitude and tropical latitude bands and in the Southern Hemisphere mid-latitudes (only with SOI).

4.7 Sensitivity to methodological choices

There are a few methodological choices in this study that we check here to assess their sensitivity to our results. In testing these sensitivities, we choose to use the correlations with ONI as the ENSO index for brevity.

4.7.1 Sensitivity to coarsening the resolution

In this case, we examine the sensitivity of the correlations of the ENSO index with the atmospheric variables at their native and coarsened 1° grid spacing. Furthermore, we analyzed this sensitivity with datasets that were available at the highest resolution on their native grids, which implies we used ERA5 for geopotential heights, CHIRPS2 for precipitation, and CRU for surface temperature. This ensured that we were testing the sensitivity of the correlations to the largest coarsening of the grids to the uniform 1° grid spacing.

Supplementary Figure S9a and b shows the sensitivity of filtered ONI correlations with 200 hPa and 500 hPa geopotential heights from ERA5, respectively. This sensitivity is expressed as the difference in the ENSO correlations between the coarsened and the native grid. In both instances, the differences in the correlations are statistically insignificant. This indicates that coarsening the ERA5 geopotential heights grid resolution to 1° did not have an impact on the ENSO correlations with mid and upper-air geopotential heights. The sensitivity of the correlations to coarsening the CHIRPSv2 precipitation from its native 0.25° grid to 1° grid had a marginal, but uniform decrease in correlation ranging between 0.02 and 0.10 (Supplementary Figure S9c). However, this change in the correlations in Supplementary Figure S9a from coarsening the CHIRPSv2 precipitation to 1° grid will only marginally impact the ENSO correlations as they appear generally stronger ≥|0.5| in Figure 4, where the correlations pass the BH test. Likewise, coarsening the grid of the surface temperature in CRU from its native resolution of 0.5° to 1° grid had negligible impact on its teleconnection with ONI. Therefore, we can conclude that coarsening the native resolution of the datasets to a uniform 1° grid did not have a significant bearing on our results.

4.7.2 Sensitivity to truncation of the number of PCs

In this case of testing methodological choices, we examine the correlations of the ENSO index with the interannually filtered atmospheric variables truncated at different number of PCs. Note that in all instances (except in case of precipitation), the interannual filtered atmospheric variables showed the most enhanced ENSO teleconnections (Figures 6–9). Therefore, we choose to use the correlations of the filtered ONI with interannually filtered variables to test this sensitivity. The correlations of the filtered 200 hPa and 500 hPa geopotential heights with the filtered ONI show a significant change (an increase) from raising the truncation of the PCs from 5 to 8 (Supplementary Figure S10). However, the differences in the correlations as the truncation of the PC is raised to 10 and 12 from 8 appears minimal (Supplementary Figure S10). In contrast, the correlations of precipitation from CHIRPSv2 with ONI appear similar across truncation of PCs at 5, 8, 10, and 12 in Supplementary Figure S11a, b, c, and d, respectively. However, the correlations of the CRU surface temperature with ONI at the different PC truncations show some notable differences, with truncation of PC at 12 showing the most robust ENSO teleconnections in the Southern Hemisphere (Supplementary Figure S11e-h). But qualitatively, the correlations of ENSO teleconnections with ONI at truncation of PC at 10 (Supplementary Figure S11g) is like truncation of PC at 12 (Supplementary Figure S11h). Noting these features of sensitivity to PC truncation to the EOF-EEMD procedure, we choose to use a truncation of 10 PCs as it was an optimal compromise between accuracy and computational efficiency. Truncation at higher PCs requires more computational time, which seems unnecessary given the insignificant differences in the results from using a higher truncation.

4.7.3 Sensitivity to the length of the data record

We recomputed the 200 hPa geopotential height correlations shown in Figure 2 using the common overlapping period (1980–2018) across all reanalyses, and the results are presented in Supplementary Figure S12. The correlation patterns (shaded regions) in Supplementary Figure S12 closely resemble those in Figure 2. However, the spatial extent of statistically significant correlations (stippled regions) associated with the filtered ONI and MEI is substantially reduced in Supplementary Figure S12 relative to Figure 2 for all reanalyses except ERA5. This marked reduction in statistical significance indicates that the 39-year overlapping period is insufficient to robustly diagnose the relationship between ENSO and 200-hPa geopotential height variations in JRA55, MERRA2, and NCEP-R2.

Similarly, restricting the analysis to the 1980–2018 period yields correlation patterns between the filtered ENSO indices and 500-hPa geopotential height (Supplementary Figure S13) that closely resemble those shown in Figure 3. Nevertheless, the spatial extent of statistically significant correlations is considerably reduced relative to Figure 3, particularly in MERRA2 and NCEP-R2. This result further demonstrates that a 39-year record is inadequate for robustly characterizing ENSO teleconnections in 500-hPa geopotential height.

Supplementary Figure S14a–f shows the correlations between the filtered ENSO indices and the two two precipitation observational datasets over the common 1981–2022 period. For both CHIRPS2 and CRU, the correlation patterns (shaded regions) and their statistical significance (stippled regions) at both zero- and one-season lag closely resemble those in the corresponding panels of Figure 4. These results suggest that a 42-year record is, in fact, sufficient to robustly diagnose ENSO teleconnections in these precipitation datasets. Moreover, the differences in ENSO teleconnections between CHIRPS2 and CRU evident in Figure 4 are also reproduced in Supplementary Figure S14.

For surface temperature reanalyses, the correlations shown in Supplementary Figure S4 were similarly recomputed for the overlapping 1980–2018 period (Supplementary Figure S15g–l). As in Supplementary Figure S4, no statistically significant correlations are found, indicating that the shorter overlapping record is just as insufficient to diagnose ENSO teleconnections in surface temperature.

5 Discussion

5.1 Implications for ENSO teleconnections

This study evaluated the sensitivity of ENSO teleconnections to several sources of uncertainty, including the choice of reanalysis and observational datasets, ENSO index, temporal record length, spatial resolution, temporal filtering, and methodological decisions such as EOF truncation and field significance testing. Rather than acting independently, these sources of uncertainty interact throughout the analysis chain. Dataset characteristics influence the representation of climate variability, which subsequently affects the diagnosed teleconnection patterns and the statistical significance of those patterns. Consequently, the robustness of an ENSO teleconnection depends not only on the strength of the physical relationship but also on how uncertainties propagate through the observational datasets and analysis methodology.

The largest source of uncertainty arises from the selection of reanalysis and observational datasets, whereas the choice of ENSO index produces comparatively smaller differences across all variables examined. This result is consistent for 500-hPa geopotential height, surface precipitation, and land surface temperature. Differences among datasets likely originate from variations in observational coverage, data assimilation techniques, model physics, and temporal coverage. These dataset-dependent differences provide the background state upon which subsequent methodological choices operate. For example, the differences in 500-hPa geopotential height teleconnections between ERA5 and MERRA-2 may reflect the use of four-dimensional variational (4D-Var) data assimilation in ERA5 versus three-dimensional variational (3D-Var) assimilation in MERRA-2. Likewise, observational uncertainties associated with sparse observing networks and retrieval methodologies contribute directly to differences in precipitation teleconnections diagnosed from CHIRPS2 and CRU, particularly over Africa and Central Asia, where rainfall analyses remain highly uncertain (Hu et al., 2016, 2018; Sylla et al., 2012). Because these structural differences exist before any statistical analysis is performed, they propagate into all subsequent diagnostic steps, including temporal filtering, EOF decomposition, and significance testing.

The interaction between dataset characteristics and temporal filtering is particularly important. The analysis demonstrates that low-frequency variability, including decadal and secular oscillations, modulates the diagnosed ENSO teleconnections in both 200- and 500-hPa geopotential height. Removing these longer temporal scales strengthens the interannual ENSO signal, indicating that unresolved low-frequency variability can partially mask teleconnection patterns when observational records are relatively short. Consequently, uncertainty associated with record length cannot be considered independently of temporal filtering. A shorter record inadequately samples low-frequency climate variability, reducing the effectiveness of separating interannual ENSO variability from background climate fluctuations. This interaction explains why filtering substantially improves the statistical significance of the CRU land surface temperature teleconnections while having comparatively little influence on precipitation, suggesting that precipitation teleconnections are less sensitive to low-frequency modulation than temperature and geopotential height.

These findings are consistent with previous studies showing that ENSO teleconnections vary with the evolving background climate state. Van Oldenborgh and Burgers (2005) argued that observational records remain too short to robustly quantify multiplicative decadal modulation of ENSO teleconnections. Likewise, several studies have demonstrated that low-frequency modes such as the Pacific Decadal Oscillation influence the amplitude and spatial structure of ENSO teleconnections (Gershunov and Barnett, 1998; McCabe and Dettinger, 1999; Misra and DiNapoli, 2013). Hoerling and Kumar (1997) further showed that individual El Niño events produce different teleconnection patterns because the associated SST anomalies differ in spatial structure and amplitude, while Chiang and Sobel (2002) proposed that tropical tropospheric temperature adjustments transmit ENSO signals throughout the tropics. Differences in the representation of these processes among datasets therefore amplify the uncertainty introduced by limited record length and temporal filtering.

The comparatively robust and spatially coherent teleconnections diagnosed in the tropical upper troposphere are consistent with the tropospheric temperature mechanism of Chiang and Sobel (2002), whereby ENSO-induced diabatic heating anomalies generate large-scale tropical temperature adjustments that extend well beyond the Pacific basin. Accordingly, the 200-hPa geopotential height teleconnections exhibit stronger zonal symmetry and greater robustness than those at 500 hPa, reflecting the more efficient propagation of equatorial wave disturbances in the upper troposphere (Walker, 1937; Wu et al., 2000). The persistence of these patterns across datasets suggests that physically robust teleconnections are less sensitive to methodological choices than weaker regional responses.

In contrast, methodological choices such as EOF truncation, spatial regridding, and field significance testing introduce relatively modest uncertainties. The sensitivity experiments demonstrate that retaining 10 principal components and regridding all datasets to a common 1° resolution have negligible effects on the diagnosed teleconnection patterns. These methodological decisions therefore do not substantially alter the large-scale ENSO signal once the underlying dataset characteristics have been established. Similarly, the Benjamini–Hochberg field significance test primarily affects the confidence assigned to the diagnosed correlations rather than the spatial structure of the teleconnections themselves. Thus, methodological uncertainties generally influence the statistical interpretation of teleconnections, whereas dataset selection and temporal sampling primarily determine the diagnosed patterns.

Collectively, these results indicate a hierarchy of uncertainty. Dataset selection and observational characteristics constitute the dominant source of uncertainty, while temporal sampling and low-frequency variability modulate the diagnosed teleconnections through their interaction with temporal filtering. Methodological choices, including EOF truncation, regridding, and field significance testing, contribute comparatively little to the diagnosed spatial patterns. Consequently, differences in ENSO teleconnections reported across studies may often reflect differences in the observational datasets and temporal characteristics employed rather than fundamental inconsistencies in the underlying physical teleconnection mechanisms. Future efforts to quantify ENSO teleconnections would therefore benefit from longer observational records, improved representation of low-frequency climate variability, and continued advances in data assimilation and observational coverage.

5.2 Limitations

Several limitations should be considered when interpreting the results of this study. First, the diagnosed ENSO teleconnections are based on linear correlation analysis, which quantifies statistical associations rather than causal relationships. Consequently, the identified teleconnection patterns should not be interpreted as direct evidence of physical causation. Instead, their physical basis is inferred from the extensive body of literature demonstrating that these statistical relationships arise from the large-scale atmospheric response to tropical Pacific SST anomalies associated with ENSO (e.g., Branstator, 1983; Chiang and Sobel, 2002; Hoerling et al., 1997; Wallace et al., 1998). Likewise, the EOF–EEMD framework represents a linear decomposition of climate variability and therefore cannot fully capture the nonlinear interactions among ENSO, decadal variability, and regional climate processes. Nevertheless, the linear framework provides a transparent and computationally efficient means of isolating the influence of different temporal scales on ENSO teleconnections, consistent with the methodology commonly adopted in teleconnection studies (e.g., Mo, 2010).

Second, the EOF–EEMD methodology itself has inherent limitations. These include mode mixing, endpoint effects arising from finite record lengths, and sensitivity to algorithmic parameters such as the number of ensemble realizations, the amplitude of the added white noise, the number of sifting iterations, and the stopping criteria used to extract intrinsic mode functions. In this study, mode mixing was mitigated by employing a relatively large ensemble size of 1,000 realizations. Endpoint effects, however, cannot be completely eliminated and remain dependent on the length of the available observations. This limitation motivated the use of the full available observational record, in addition to analyses restricted to the common overlapping period among datasets, to better resolve the low-frequency modes of variability.

Third, computational limitations constrained the scope of the scale-separation analysis. Because EOF–EEMD analysis over a global domain is computationally demanding, the geopotential height analysis was performed only for ERA5, the most recent reanalysis considered in this study. Similarly, the precipitation scale-separation analysis was restricted to CHIRPS2. Consequently, the diagnosed effects of temporal scale separation should be regarded as dataset-specific and should not be generalized to the remaining reanalysis or observational datasets without comparable analyses. While the sensitivity experiments indicate that several methodological choices, such as EOF truncation and spatial regridding, have only a minor influence on the diagnosed teleconnections, the extent to which temporal scale separation affects other datasets remains an open question.

Finally, the present study focuses primarily on quantifying uncertainties in linear ENSO teleconnections across datasets, ENSO indices, and temporal scales. It does not investigate the underlying dynamical mechanisms responsible for the diagnosed differences. Complementary analyses using nonlinear statistical techniques or dynamical diagnostics, such as Rossby wave source analyses (e.g., Nie et al., 2019), could provide additional insight into how dataset characteristics and low-frequency climate variability influence ENSO teleconnections. Such analyses represent an important avenue for future research and would help determine the extent to which the conclusions presented here can be generalized beyond the datasets and methodologies examined.

6 Summary and conclusions

ENSO-driven climate variability is one of the most important sources of seasonal climate variability and associated societal impacts worldwide. Consequently, ENSO teleconnections are widely used to understand regional climate variability and are frequently employed as benchmarks for climate model evaluation. The primary objective of this study was to systematically quantify the sensitivity of boreal winter ENSO teleconnection diagnostics to the choice of observational and reanalysis datasets, meteorological variables, ENSO indices, temporal record length, and methodological choices. In contrast to many previous studies, the statistical significance of the diagnosed teleconnections was evaluated using the Benjamini–Hochberg field significance test, thereby accounting for false discoveries associated with multiple spatially dependent hypothesis tests.

This study provides one of the first systematic assessments of how these factors influence ENSO teleconnection diagnostics across commonly used validation datasets. Specifically, we (1) compared teleconnection patterns among upper-air and surface variables, (2) evaluated differences across multiple reanalyses and observational datasets, (3) examined the sensitivity to commonly used ENSO indices, and (4) assessed the influence of temporal-scale separation using an EOF–EEMD framework, although the latter analysis was necessarily restricted to ERA5 and CHIRPS2 because of computational constraints.

The principal findings are summarized as follows:

• Although ONI, MEI, and SOI are strongly correlated with one another, they produce measurably different teleconnection patterns, illustrating that high correlations among ENSO indices do not imply equivalent relationships with atmospheric variables.
• The dominant source of uncertainty in ENSO teleconnection diagnostics is the choice of reanalysis or observational dataset. Differences among datasets consistently exceed those arising from the choice of ENSO index, particularly for 500-hPa geopotential height, land surface temperature, and precipitation. These differences likely reflect variations in observational coverage, data assimilation methodologies, and model formulations.
• The second most important source of uncertainty is the interaction between interannual ENSO variability and climate variability on other temporal scales. Removing decadal and secular variability substantially strengthens the interannual ENSO teleconnections in ERA5 geopotential height and CRU land surface temperature, demonstrating that unresolved low-frequency variability can mask ENSO signals, particularly when observational records are relatively short.
• Record length therefore represents another important source of uncertainty because it determines the extent to which low-frequency variability can be adequately sampled and separated from the interannual ENSO signal. The common observational period was found to be insufficient for robustly characterizing ENSO teleconnections in most atmospheric variables, although precipitation appeared comparatively insensitive to this limitation.
• In contrast, methodological choices such as EOF truncation and regridding to a common 1° resolution exert only minor influences on the diagnosed teleconnections, indicating that these processing decisions are considerably less consequential than dataset selection or temporal sampling.
• Among the atmospheric variables examined, tropical upper-tropospheric (200-hPa) geopotential height exhibits the most robust and spatially coherent ENSO teleconnections across datasets and ENSO indices, whereas mid-tropospheric circulation and surface variables display substantially larger structural uncertainty.

Taken together, these results establish a clear hierarchy of uncertainty in observational ENSO teleconnection diagnostics. Dataset selection is the largest contributor to uncertainty, followed by temporal sampling and the interaction of ENSO with lower-frequency climate variability, whereas ENSO index selection exerts a comparatively smaller influence, and methodological choices such as EOF truncation and spatial regridding contribute the least. Consequently, differences in ENSO teleconnections reported across studies are more likely to arise from differences in the underlying datasets and temporal characteristics than from the particular ENSO index employed.

These findings have important implications for both observational analyses and climate model evaluation. Because observational teleconnections themselves are subject to substantial structural uncertainty, evaluating climate models against a single observational or reanalysis product may lead to misleading conclusions regarding model performance. Multi-dataset observational benchmarks should therefore become standard practice, particularly for surface variables and mid-tropospheric circulation fields. Likewise, future observational analyses would benefit from longer, more homogeneous climate records that better resolve low-frequency variability, continued improvements in data assimilation and observational coverage, and the incorporation of complementary dynamical and nonlinear diagnostics to better understand the physical mechanisms responsible for differences among datasets. Such advances will improve the robustness of observational benchmarks and provide a stronger foundation for evaluating the representation of ENSO teleconnections in present and future climate models.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fclim.2026.1830449/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Matthew Collins, University of Exeter, United Kingdom

Reviewed by: Lin Chen, Nanjing University of Information Science and Technology, China

Paulo Miguel De Bodas Terassi, University of São Paulo, Brazil