Identifying quasi-periodic variability using multivariate empirical mode decomposition: a case of the tropical Pacific

A variety of statistical tools have been used in climate science to gain a better understanding of the climate system's variability on various temporal and spatial scales. However, these tools are mostly linear, stationary, or both. In this study, we use a recently developed nonlinear and nonstationary multivariate time series analysis tool – multivariate empirical mode decomposition (MEMD). MEMD is a powerful tool for objectively identifying (intrinsic) timescales of variability within a given spatio-temporal system without any timescale pre-selection. Additionally, a red noise significance test is developed to robustly extract quasi-periodic modes of variability. We apply these tools to reanalysis and observational data of the tropical Pacific. This reveals a quasi-periodic variability in the tropical Pacific on timescales ∼ 1.5–4.5 years, which is consistent with El Niño–Southern Oscillation (ENSO) – one of the most prominent quasi-periodic modes of variability in the Earth’s climate system. The approach successfully confirms the well-known out-of-phase relationship of the tropical Pacific mean thermocline depth with sea surface temperature in the eastern tropical Pacific (recharge–discharge process). Furthermore, we find a co-variability between zonal wind stress in the western tropical Pacific and the tropical Pacific mean thermocline depth, which only occurs on the quasi-periodic timescale. MEMD coupled with a red noise test can therefore successfully extract (nonstationary) quasi-periodic variability from the spatio-temporal data and could be used in the future for identifying potential (new) relationships between different variables in the climate system.

1 Introduction

The climate system is a highly complex system consisting of variability across many different timescales (e.g. Baede et al., 2001). Among these, quasi-periodic patterns of variability are important as they can be leveraged for medium-to-long-range predictions (e.g. Mariotti et al., 2018, 2020; L'Heureux et al., 2020). Statistical tools have been particularly helpful in identifying quasi-periodic variability.

The following are typical statistical tools used for exploring the patterns of variability on different temporal and/or spatial scales (Broomhead et al., 1987; Hasselmann, 1988; Penland and Sardeshmukh, 1995; Ghil et al., 2002; Huang et al., 1998; Ghil et al., 2002; Froyland et al., 2021): Fourier transform (FFT), wavelet transform, principal component analysis (PCA), (multi-channel) singular spectrum analysis ((M)SSA), principal oscillation patterns (POPs), linear inverse model (LIM), or even nonlinear Laplacian spectral analysis. However, these methods do not necessarily have a multivariate extension and/or are either stationary, linear, or both. Additionally, some methods, such as (M)SSA, require specification of a window over which quasi-periodic signals are sought, making them somewhat subjective. These can be drawbacks when studying complex, inherently nonlinear, and nonstationary systems, such as the climate system.

Multivariate empirical mode decomposition (MEMD) addresses these drawbacks as it is an analysis tool that is entirely data adaptive and is designed to extract nonlinear and nonstationary signals. MEMD is a generalisation of the empirical mode decomposition (EMD; Huang et al., 1998) to multivariate datasets of more than two time series (Rehman and Mandic, 2010). EMD is a 1-D time series analysis tool that is based on Hilbert transform and takes advantage of the instantaneous frequency, allowing a local extraction of modes of variability. Each mode that the EMD extracts consists of two elements: (1) typical timescale of the mode, i.e. average instantaneous frequency of the mode; and (2) the time series of the mode. Beyond this, MEMD extracts timescales common to all input time series (i.e. synchronises signals) and provides modes (patterns) of variability according to these timescales.

Despite their appeal, MEMD and EMD have hardly been used in climate research. EMD and its 1-D extension ensemble EMD (EEMD; Wu and Huang, 2009) have been used for smoothing, filtering, extracting trends, variability, and testing for red noise distribution of climate data (e.g. Duffy, 2004; Wu et al., 2007; Franzke, 2009; Lee and Ouarda, 2011; Qian et al., 2011; Franzke and Woollings, 2011; Franzke, 2012; Ezer and Corlett, 2012; Ezer et al., 2013; Wang and Ren, 2020). MEMD has only been used for an idealised analysis of atmosphere–ocean coupling strength (Alberti et al., 2021). Moreover, neither MEMD nor EMD has been used for extracting quasi-periodic modes of climate variability.

A major challenge in applying MEMD in climate analysis is that no statistical null hypothesis test for red noise has been developed. When applied to climate data, MEMD can reveal many modes that are consistent with red (or white) noise. In particular, sea surface temperature (SST) exhibits a red spectrum because it represents the integral response of the ocean to stochastic higher-frequency atmospheric (e.g. weather or white noise) variability (e.g. Hasselmann, 1976; Frankignoul and Hasselmann, 1977). However, there are also patterns of quasi-periodic variability (e.g. repeating every 2–8 years) that reflect more complex climate dynamics. To identify such quasi-periodic variability in the climate system, any analysis would first require a significance test that could distinguish this variability from the red noise. In other words, we seek spectral peaks (modes of variability) that pass the red noise threshold (Gilman et al., 1963; Madden and Julian, 1971; Bretherton et al., 1999; see also Sect. 4, Appendix B2).

In this study, we combine the MEMD method with a red noise test (Sects. 3 and 4) to robustly detect quasi-periodic modes of variability. Thus, MEMD becomes well suited for analysing nonlinear and nonstationary climate data. It also has the advantage of objectively detecting intrinsic timescales without pre-selecting and/or filtering for a frequency band (e.g. as in MSSA – multichannel singular spectrum analysis) and without any periodic signal or basis of functions specifications (e.g. as in Fourier transform and wavelet analysis). Since MEMD and (especially) its combination with a red noise test is a new tool in climate science, we test it on a known example of quasi-periodic variability. Namely, we analyse the tropical Pacific atmosphere–ocean variability to extract the El Niño–Southern Oscillation (ENSO). In the future, this tool may be used in other spatio-temporal applications where quasi-periodic variability has not yet been identified.

ENSO is a quasi-periodic phenomenon occurring on (interannual) timescales of 2–8 years (e.g. Philander, 1989; Wang and Fiedler, 2006; Timmermann et al., 2018). ENSO events are characterised by warming tropical Pacific SST during the development of El Niño (warm phase) and cooling of SST afterwards leading into La Niña (cold phase). These events are typically characterised by ocean–atmosphere interaction, whereby atmospheric changes in winds can lead to changes in the distribution of warm and cold waters in the ocean that in turn impact the atmosphere. ENSO exhibits significant nonlinearities, with marked skewness and phase locking to the seasonal cycle (e.g. Stein et al., 2010; Dommenget et al., 2013). It is also nonstationary (e.g. Crespo et al., 2022; Fedorov and Philander, 2000). This complex coupled dynamics is different from red noise and is therefore the focus of the present study.

The paper is structured as follows. Section 2 provides description of data used; Sect. 3 describes the EMD and MEMD implementation (Appendix A provides further details); Sect. 4 provides description of the red noise test (Appendix B provides further details); Sects. 5 and 6 identify modes of variability in the tropical Pacific and explore the physical mechanisms relevant on different timescales. Conclusions are given in Sect. 7.

2 Data

We focus on the intrinsic variability of the tropical Pacific and analyse monthly mean data of three different variables relevant for atmosphere–ocean exchange: sea surface temperature (SST) from the HadISST observational dataset (Rayner et al., 2003), as well as surface zonal wind stress (τ_x) and thermocline depth (i.e. the depth of the 20 ^∘C isotherm) – both from SODA2 ocean-reanalysis dataset (Carton and Giese, 2008). These three variables are used to assess the ENSO dynamics as these quantities typically play a leading role in the onset and decay of ENSO events and are typically used in the oscillator models that explain ENSO dynamics (see Sect. 6 and e.g. Wang, 2018). The data are analysed in the tropical Pacific (25^∘ S–25^∘ N, 110^∘ E–65^∘ W) over the period 1871–2010 for which all datasets are available. Note that surface wind stress and the ocean subsurface data reconstructions in the 19th and early 20th century are less reliable than in the late 20th century due to sparser data coverage; thus, the results presented here are only as accurate as these reconstructions can be (Wittenberg, 2004; Crespo et al., 2022). We use the early record data here to show the nonstationarity and nonlinearity of the tropical Pacific variability over the last ∼ 140 years. The identification of the two main quasi-oscillatory modes (Sect. 5) is also not affected by the inclusion of these data.

The MEMD analysis (Sect. 3) is performed on all fields simultaneously with SST at the highest resolution (1^∘ in latitude and longitude), whereas the thermocline depth and τ_x (both 5^∘ resolution in latitude, 9^∘ resolution in longitude) have much lower resolution. This gives greater weight to SST data in the analysis and less towards the other variables, such that the mode does not change significantly by adding other variables in the analysis (i.e. results below for the SST are similar whether we use SST alone or together with other fields). This choice is made because we are primarily interested in the quasi-periodic behaviour in SST, since SSTs are used for defining ENSO (e.g. in the Niño3 region – see Table 1). SSTs are also typically smoother than other fields. This is especially true for wind stress, which is strongly affected by the noisy atmospheric variability. This ensures that modes that emerge from MEMD analysis are representative of quasi-periodic variability in SST (and thus ENSO), while the rest of the variables are enslaved to SST variability. The other variables are added to the MEMD input data to describe the climate dynamics involved in the quasi-periodic SST variability, e.g. ENSO (Sect. 6). Note that MEMD can be sensitive to input data; thus, we must carefully consider the input data structure (relevant to a specific study).

Table 1 Tropical Pacific regions used for computing time series (see text for details). The rightmost column lists the variables that are averaged over specified regions.

While the SST, thermocline depth, and τ_x play an important role in the ENSO dynamics, it is specific regions (see Table 1) that have historically been analysed in more detail, for example, in conceptual oscillator models (e.g. Jin, 1997a; Wang, 2018). The time series in the specific regions are thus used to assess if MEMD modes on ENSO timescales are consistent with physics described by conceptual oscillator models (Sect. 6). Additionally, analysing such time series helps with a simpler visualisation of temporal evolution of different variables. Thus, we average τ_x over Niño4 and Niño5 regions separately, the thermocline depth over the Niño6 region (off-equatorial thermocline depth) and over the tropical Pacific (Pacific mean), and the SST over the Niño3 region (again, see Table 1; see also Fig. 3 in Wang et al., 1999). Note that regional averages (Table 1) are computed after the MEMD analysis is performed (for details see Sect. 3).

Before performing the analysis, we detrend the data and remove their seasonal cycle, which is done the following way (see also, e.g., de la Cámara et al., 2019). First, we calculate 30-year means centred on every 10th year for each individual month. This yields one value for each individual month every 10 years. Then we interpolate between these values (of every 10th year) to obtain yearly time series, again for each individual month. This yields a smooth seasonal cycle that includes a trend and seasonal cycle for every month in the record. Detrended and deseasonalised data are then computed as the difference between the original monthly time series and the smooth trend/seasonal cycle. We do this at every grid point and for every variable separately.

This is done to avoid domination of the seasonal cycle or trend in the statistical analysis below, even though the MEMD can generally extract nonlinear trends by itself. Note that this means that we cannot assess the impact of long-term variability or seasonal cycle on ENSO variability in this study, but the latter may still be present indirectly, as ENSO is phase-locked to the seasonal cycle (e.g. Stein et al., 2010; Wengel et al., 2018).

3 (Multivariate) empirical mode decomposition

3.1 MEMD description

We employ MEMD to objectively identify intrinsic modes of variability in nonlinear and nonstationary spatio-temporal data. To understand MEMD, it is easier to first consider the simpler implementation of the 1-D version, empirical mode decomposition (EMD), as outlined by Huang et al. (1998):

• i.

  Local minima and maxima of the input (1-D) time series (e.g. see black solid line in Fig. 1) are identified.

• ii.

  Envelopes are created by interpolating between subsequent maxima (upper envelope, shown as grey dotted line in Fig. 1) and between subsequent minima (lower envelope, represented by grey dashed line in Fig. 1).

• iii.

  An average envelope is obtained by taking the mean of the upper and lower envelopes (depicted by the red solid line in Fig. 1).

• iv.

  The average envelope is subtracted from the original time series data.

• v.

  The subtracted data, i.e. the original data minus the average envelope, represent the first mode of variability and typically correspond to the highest-frequency signal in the dataset. However, the average envelope can be further analysed.

• vi.

  Steps i–vi are repeated on the average envelope until only a trend or residual component remains. This occurs when we can no longer find at least two extrema in the dataset, which is a condition that needs to be satisfied by EMD's modes of variability.

The resulting modes of variability (i.e. the output of the EMD analysis), represented by their respective time series, are known as intrinsic mode functions (IMFs). The timescale of each IMF is characterised by the time lapse between two subsequent extrema, and the mean timescale of an IMF represents the average over all time lapses within its time series.

Figure 1 Schematic for obtaining the average envelope during the EMD process. The black line shows a simple input signal that is a sum of two sine waves (a low-frequency and a high-frequency wave). Grey dotted and grey dashed lines show the upper and lower envelope, respectively. The red line shows the average envelope, which represents the low frequency of the input signal. If we remove the average envelope from our input data, we obtain the high-frequency signal (e.g. the first mode of EMD).

Additionally, each IMF has to satisfy two criteria: (a) the number of extrema and the number of zero crossings differs at most by one and (b) the mean value of the envelope of the IMF is zero. Note that the procedure from i to iv does not necessarily satisfy (a)–(b) immediately; thus, an additional sifting process (typically iterative) is used, which requires a stopping criteria to ensure the physical meaning of the IMFs. The stopping criteria can be based on the standard deviation of each IMF and on the maximum number of iterations, among others, which set tolerance and confidence limits for the IMF (for details see Huang et al., 1998; Rilling et al., 2003; Huang et al., 2003). Ultimately, this process extracts modes of variability that consist of (1) the typical timescale of the mode and (2) its time series. The modes are automatically ordered from the highest (shortest) in the first mode to lowest (longest) frequency (period) in the last mode. Note that modes that emerge are largely independent with only small correlations between them.

MEMD (Rehman and Mandic, 2010) is a generalisation of the EMD to multivariate datasets of more than two time series (for bivariate and trivariate data separate methods exist; Rilling et al., 2007; Rehman and Mandic, 2010). The method solves a similar problem as in i–vi, but the mean envelope is computed as (Rehman and Mandic, 2010) “an integral of all the envelopes along multiple directions in an N-dimensional space” (i.e. on an N sphere). This is much more complex, but the basic idea remains similar to the 1-D method, and the method retains similar stopping criteria for the sifting process. For further details on and visualisations of the method the reader is referred to Rehman and Mandic (2010) and Alberti et al. (2021).

The MEMD ultimately extracts timescales common to all input time series (i.e. synchronises signals; Rehman and Mandic, 2010) and provides multivariate IMFs (i.e. the outputs of MEMD method) describing those timescales. The timescales of IMFs are then consistent across the input time series – to visualise this, see Supplement Figs. S6–S10, Tables S1–S2, and Sect. S.3 for an idealised example.

As with all statistical methods, it is important to be aware of the drawbacks associated with the (M)EMD (e.g. Stallone et al., 2020). Similar to other time-series-filtering methods, (M)EMD can encounter issues at the edges of the time series, which can result in travelling waves and thus unrealistic peaks in the time series (e.g. Stallone et al., 2020). Another common challenge with (M)EMD is the mixing of modes (see below) where a single genuine mode may be split into multiple modes if inappropriate parameters are chosen (e.g. Huang et al., 1999, 2003; Stallone et al., 2020). However, this issue can also arise when there is no clear timescale separation. To address these challenges, it is crucial to test the physical relevance of the modes of variability identified using (M)EMD and to ensure convergence and stability of the modes through different parameter sweeps that are related to the stopping criteria (see Sect. 3.2). Note that different parameter sweeps may be relevant for different applications.

3.2 MEMD parameter sweep

The code for the MEMD is freely available on GitHub (https://github.com/mariogrune/MEMD-Python- (last access: 1 December 2020); similarly for the EMD discussed above: https://github.com/laszukdawid/PyEMD, last access: 4 December 2023), and the user ultimately only decides about the stopping criteria. These are set to the “fix_h” parameter following Huang et al. (2003), who suggest that limiting iterations yields better-behaved IMFs than other stopping criteria. We limit the number of iterations to 15 (parameter “n_iter” is 15), though other values were tested and a range for “n_iter” around 5–30 yielded similar results, suggesting some convergence for the significant modes of variability (see below).

Note that at higher/lower frequencies we find mode mixing in our MEMD analysis where timescales are not clear (here, this occurs on timescales shorter than about 8 months and longer than about 700 months), as well as with a larger number of iterations. These modes are not detected as different from red noise (see Sects. 4 and 5 for more details). However, the significant modes of variability on interannual timescales that are of interest here are largely unaffected by this (Sect. 5).

3.3 Obtaining modes of variability via MEMD

As mentioned in Sect. 2, we use 3-D data relevant for ENSO dynamics, i.e. $\mathbf{A}(L_t,L_y,L_x)$ (with L_t the length of the time dimension, L_y the length of the latitude dimension, L_x the length of the longitude dimension, and A the selected variable(s)/field(s); bold letters represent two- or three-dimensional arrays) as the basis for the MEMD analysis, which is done the following way.

• i.

  First, we remove the smooth trend/seasonal cycle (described above; see also de la Cámara et al., 2019) from A to get a temporal anomaly, A^′.

• ii.

  We divide data by their standard deviations (σ; done separately for each variable and at each grid point).

• iii.

  Since we use more than one variable (i.e. SST, surface wind stress, thermocline depth) in the analysis, we concatenate (denoted with ⊕) the different variables in their spatial dimensions; i.e. we get ${\mathbf{A}}^{\prime }(L_t,L_y=L_{y\mathrm{1}}\oplus L_{y\mathrm{2}}\oplus \mathrm{\dots },L_x=L_{x\mathrm{1}}\oplus L_{x\mathrm{2}}\oplus \mathrm{\dots })$ (with subscripts 1, 2, etc., representing spatial dimensions of the different variables).

• iv.

  Then, we reduce the dimensionality by computing spatial patterns (empirical orthogonal functions, EOFs) and their time series (principal components, PCs) via singular value decomposition (SVD).

• v.

  Finally, we only retain the first 20 PCs that explain the majority of the variance in the field A^′. The PCs are then used as input data for the MEMD algorithm (for further details on the whole procedure see Appendix A).

MEMD analysis identifies 21 IMFs that are ordered by frequency from the highest (IMF 1) to the lowest (IMF 21) with the last 21st mode typically representing a trend, which in our case was already removed (see above). Namely, we find 21 potential intrinsic timescales within the tropical Pacific, i.e. common to all input PCs. This means that we obtain 21 IMFs for each PC time series, i.e. PC${}_m\left(L_t\right)={\sum }_s{\mathrm{IMF}}_s\left({\mathrm{PC}}_m\right(L_t\left)\right)$, with s the IMF number and m the PC number, where the sth IMF of each PC_m has the same timescale (see Table S2). Since we initially obtained PCs via SVD analysis and we thus have corresponding EOF patterns, we can then reconstruct spatio-temporal patterns of variability for each field (i.e. SST, wind stress, and thermocline depth) from PCs and EOFs. The interest here is in the variability on different timescales (i.e. for each (sth) IMF separately, we can perform this reconstruction for each IMF across all 20 PCs/EOFs):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathrm{IMF}}_s^{\mathrm{spatial}}(L_t,L_y,L_x)\\ & \simeq \sum _m{\mathrm{IMF}}_s\left({\mathrm{PC}}_m\right(L_t\left)\right){\mathrm{EOF}}_m(L_x,L_y).\end{array}\end{array}$$

Please recall that input data for MEMD analysis were divided by σ; thus, for the variability of a field in its original units we need to multiply spatio-temporal IMFs by σ.

To compute an index, such as the eastern Pacific SST (Niño3), we can average over a latitude–longitude region (Table 1) from spatio-temporal IMFs (Eq. 1) to obtain time series of, for example, the eastern Pacific SST (Niño3) for each IMF separately – see Figs. S6–S10 for time series of different variables from Table 1 and Table S1 for their timescales. This yields an index for each IMF (e.g. IMF_s(SST (Niño3)), with s the IMF number) corresponding to an equivalent index computed from input data (A^′). The latter is approximately the same as the sum of indices computed from all IMFs, e.g. $\mathrm{SST}\left(\mathrm{Ni}\text{ñ}\mathrm{o}\mathrm{3}\right)\left(L_t\right)\simeq {\sum }_s{\mathrm{IMF}}_s\left(\mathrm{SST}\right(\mathrm{Ni}\text{ñ}\mathrm{o}\mathrm{3}\left)\right(L_t\left)\right)$.

4 Statistical significance test for climate

Once we have computed the IMFs, we need to test if they are statistically significant. The importance of each IMF can be assessed by computing the explained variance of each IMF relative to the input field (e.g. retaining those IMFs that explain more than 0.1 % variance) or through other significance tests (e.g. white noise test; Appendix B1; Wu and Huang, 2004). We develop a test for variability that can be distinguished from red noise, as appropriate for studying climate variability (see Introduction).

The red noise test can be performed in different ways (see also Sect. 5.1): (i) we can choose an index of interest and perform a red noise test on its IMFs (see below; Fig. 2) – this can be used when we are interested in the quasi-oscillatory behaviour of a specific index (e.g. areal average) we constructed or a known climate index we wish to test; (ii) more objectively, we can use IMFs of specific PC time series (e.g. PC1; Fig. S2) – this can be used when we are interested in quasi-oscillatory behaviour of an index that explains a specific amount of variance in a field or region (e.g. tropical Pacific); and (iii) we can test for quasi-periodic variability at each grid point (Sect. 5.1, Fig. 5) – this can be used to identify potential regions of quasi-oscillatory behaviour that have not been identified before. In general, all methods can be used together (as done here) for robustness, but their use will ultimately depend on the specific scientific question.

Figure 2 The red noise significance test for spatio-temporal IMFs (from MEMD) via SST averaged over the Niño3 region (for details see text and Appendices A and B). Each blue dot represents the average squared amplitude (E_s; Eq. B1) and average period (T_s; in months; for periods of IMFs see Table S1) of each IMF that we have identified within our time series/data. The average period is computed from the instantaneous frequency obtained via Hilbert transform (see text under Eq. B4). For visualisation purposes we obtain a natural logarithm of both the average squared amplitude (log _e(E_s); y axis) and the average period (log _e(T_s); x axis) of each IMF and plot them as a scatterplot (blue dots). Note that the logarithms of periods (log _e(T_s)) are ordered from the shortest period (highest frequency; IMF 1; leftmost blue dot) to the longest period (trend; IMF 21; rightmost blue dot). The red solid line represents the theoretical red spectrum fit (Eqs. 2 and B7–B9), and the red dotted line represents the one-tailed 95th percentile confidence level (via χ² test). For further descriptions of the figure see text.

Here, we first test for potential quasi-periodic variability using SST time series that are relevant for (eastern Pacific) ENSO variability, i.e. the eastern Pacific SST (Niño3) from input data and the corresponding eastern Pacific SST (Niño3) from spatio-temporal IMFs (as described in Sect. 3; IMF_s (SST (Niño3))). Then, we compute the power spectrum for each timescale/period (i.e. each IMF) by obtaining the average squared amplitude of each (sth) IMF (E_s), i.e. $E_s=\left({\sum }_j\right[{\mathrm{IMF}}_s\left(j\right){]}^{\mathrm{2}})/L_t$ (with L_t the length of the time series and j each time step). The average timescale (T_s) of each IMF is computed from the instantaneous frequency using Hilbert transform (see text around Eq. B4 in Appendix B1), which yields similar results to computing the time lapse between two extrema in the time series. E_s is then plotted against T_s (here we use their logarithmic values) to yield the power spectrum plot (blue dots in Fig. 2).

The shape of the red noise fit (red solid line in Fig. 2), or rather red noise spectrum for every sth IMF, can be computed from the lag-one autocorrelation of the input data (e.g. eastern Pacific SST (Niño3)) as (see also Gilman et al., 1963; Kolotkov et al., 2016; for full derivation and further discussion see Appendix B)

$$\begin{array}{}\text{(2)}& E_s^{\mathrm{red}}\left({\mathit{\nu }}_s\right)={\displaystyle \frac{\mathrm{1}-r^{\mathrm{2}}}{\mathrm{1}-\mathrm{2}r\cos \mathrm{2}\mathit{\pi }{\mathit{\nu }}_s+r^{\mathrm{2}}}}({\mathit{\nu }}_s{\mathit{\alpha }}_s-{\mathit{\nu }}_s/{\mathit{\beta }}_s).\end{array}$$

$E_s^{\mathrm{red}}$ is a theoretical estimate of the (mean) energy of the red noise of the sth IMF, $\mathit{\nu }=\mathrm{1}/t$ is frequency, t is time, r is the lag-1 auto-correlation of input data of the eastern Pacific SST (Niño3), ${\mathit{\beta }}_s=\sqrt{{\mathit{\nu }}_s/{\mathit{\nu }}_{s+\mathrm{1}}}$ and ${\mathit{\alpha }}_s=\sqrt{{\mathit{\nu }}_{s-\mathrm{1}}/{\mathit{\nu }}_s}$, and subscript s represents the sth IMF of frequency ν_s (ordered from highest to lowest frequency).

$E_s^{\mathrm{red}}$ is then scaled so that the variance of the red noise fit is identical to the variance of the spectra computed from the spatio-temporal IMFs (E_s; Eq. B9; see also Madden and Julian, 1971; Bretherton et al., 1999). The 95 % confidence curve (red dashed line in Fig. 2) is computed by multiplying the red noise fit by the χ² statistic for a one-sided p value of 0.05. The number of degrees of freedom for the sth IMF is computed as ${\mathrm{DoF}}_s=E_s{L}_t(\mathrm{1}-r^{\mathrm{2}})/(\mathrm{1}+r^{\mathrm{2}})$ (for further details see Appendix B2; see also Bretherton et al., 1999; Wu and Huang, 2004; Kolotkov et al., 2016). Note that the same procedure (as described above) can be applied to any other time series/IMFs.

5 Tropical Pacific modes of variability

As MEMD in conjunction with a red noise test has not been applied in climate science before, we perform an extensive analysis of the method itself (in addition to an analysis of ENSO dynamics – see below). Then, we compare it to the basic band-pass filtering (fifth-order Butterworth filter) and to Fourier transform analysis to ensure consistency with other methods. Please recall that while other spectral methods often require prior knowledge about the spatial/temporal patterns of interest in order to construct appropriate indices, the MEMD method allows for the objective extraction of significant patterns and modes of variability from data without pre-existing knowledge (see also below and Appendices A and B).

5.1 Significant modes of variability

We identify two significant modes of variability in the eastern Pacific: the 11th IMF (IMF 11) and 12th IMF (IMF 12) (shown with blue arrows in Fig. 2) with average timescales of ∼ 23 and ∼ 39 months (∼ 2–3 years; see also Table S1), respectively. Their timescales fall well within the typical ENSO timescale range (2–8 years), as can be identified via a typical power spectrum analysis (Appendix B2, Fig. B3). These timescales (with their uncertainty ranges; see below) are consistent with previously identified quasi-biennial (QB; ∼ 2 years) and low-frequency/quasi-quadrennial (LF/QQ; ∼ 4 years) ENSO modes (e.g. Jiang et al., 1995; Allan, 2000; Kim et al., 2003; Keenlyside et al., 2007; Bejarano and Jin, 2008; Jajcay et al., 2018; Froyland et al., 2021).

The quasi-periodic mode of variability with a 39-month average period (IMF 12) is more statistically significant (than IMF 11). This can be seen in different ways. First, considering the Niño3 index (i.e. the eastern Pacific SST in the Niño3 region), IMF 11 (blue dot) lies only slightly above the 95 % red noise threshold (red dashed line), whereas IMF 12 lies well above the threshold (Fig. 2). Second, considering the PC1's IMFs (i.e. before we do spatio-temporal reconstruction of IMF data; Fig. S2), only IMF 12 is found to be significant, despite the generally good agreement between the Niño3 index and PC1 of the tropical Pacific SST (e.g. Ashok et al., 2007). Finally, we can perform a red noise test on each grid point of the tropical Pacific (similarly as for specific time series shown above) and plot the spatial structure of SST together with the significance test (Sect. 6, Figs. 5, S3). This shows that none of the grid points of spatio-temporal IMF 11 pass the red noise threshold (Fig. S3), but many grid points of the spatio-temporal IMF 12 are significant (not shaded in grey in Fig. 5).

Nevertheless, even though IMF 11 is marginally statistically significant, the MEMD analysis suggests that there may be two quasi-periodic modes of interannual SST variability in the tropical Pacific region. This agrees with previous results that have identified two statistically significant modes of variability with specific frequencies and distinguished by their dynamics. Note that we also find that the two modes of variability are well separated (in terms of timescale) from the other modes and from each other, as shown by Fig. 2. This ensures no mode mixing on the quasi-periodic timescales.

On longer timescales, we do not find any behaviour that would be discernible from red noise, suggesting that the lower-frequency range of ENSO (timescales longer than ∼ 4.5 years) (e.g. Allan, 2000; Jajcay et al., 2018) is better represented by red noise. Recall, however, that any potential oscillatory behaviour on timescales longer than 30 years was removed via detrending. Similarly, we do not find any quasi-periodic modes of variability on shorter timescales. Also, on these longer (> 700 months) and shorter (< 8 months) timescales we find some mode mixing. This is likely a consequence of inputting several different variables (i.e. SST, τ_x, thermocline depth) into the MEMD algorithm, where each variable can have different timescales represented. For example, wind stress can be much noisier than SST (or different parts of the tropical Pacific have different variabilities) and can thus lead to identification of several high-frequency modes of variability, resulting in mode mixing (see somewhat overlapping blue dots in Fig. 2).

5.2 Time series of significant modes of variability

Above we have established that the Niño3 index in the tropical Pacific exhibits two quasi-periodic modes of variability with average periods ∼ 2–3 years (via MEMD analysis). We now analyse the associated time series of the two IMFs (IMF 11 and IMF 12) and compare them with a band-pass-filtered (16–53 months) Niño3 index (Fig. 3).

Figure 3 Time series (1871–2010) of the Niño3 index from IMF 11 (red dotted line) and IMF 12 (red dashed line) obtained via MEMD. The sum of the two IMF indices (black dashed line) and the band-pass-filtered (16–53 months) Niño3 index (black solid line) are also shown. All data are standardised.

The period/frequency of the two modes is not constant (i.e. varies in time). Thus, we also specify a range of periods/frequencies for the two modes. The mean periods of IMF 11 and 12 with their uncertainty ranges (in square brackets) are: 23 [16, 33] (IMF 11), and 39 [29, 53] (IMF 12) months. These ranges are defined based on the 6.7th and 93.3rd percentiles of IMF 11 and IMF 12's instantaneous period/frequency values, which roughly correspond to ±1.5σ(instantaneous periods). Figure S1 shows power spectra of IMF 11 and IMF 12, which visualise the range of periods the two modes have. Recall that the instantaneous period/frequency can be obtained via Hilbert transform (see Eq. B4 in Appendix B1). This range was chosen as it captures most of the variability in a given mode. Additionally, this range yields a good agreement with other methods (e.g. band-pass filtering), although other (reasonable) percentile ranges give qualitatively similar results. We then use these period/frequency ranges to perform a band-pass filter (via fifth-order Butterworth filter) of the Niño3 index that is consistent with the individual modes (e.g. 29–53-month band of IMF 12 is used to construct band-pass-filtered Niño3 index composites in Fig. 6b; more below). We also use a band-pass filter of the Niño3 index consistent with the sum of the two modes where the band-pass range encompasses periods of both significant modes, i.e. 16–53 months (e.g. used to construct black solid line in Fig. 3).

Time series of the sum of IMFs (black dashed line) and the band-pass index (black solid line) largely agree (i.e. their correlation is 0.95; Fig. 3). Also, the individual modes show very good agreement with the band-passed index, with correlations of 0.68 (IMF 11) and 0.83 (IMF 12), which can be increased further if we consider only the specific IMF's timescale range (specified above) when band-pass filtering the Niño3 index. This merely confirms that MEMD results are consistent with other filtering methods.

Additionally, time series of the Niño3 index extracted from IMF 11 and 12 (Fig. 3) are largely consistent with the modes identified in Jiang et al. (1995) (their Fig. 9a), who used MSSA analysis, and with Wang and Ren (2020) (their Fig. 3), who used EEMD on the Niño3.4 index. IMF 12 is also similar to the Froyland et al. (2021) 4-year mode (their Fig. 10), which used an operator-theoretic approach. The average periods of the significant modes of variability (IMF 11 and 12) in this study are typically lower than in other studies; however, please recall that we have used much longer time series and that the time series of IMFs are nonstationary. Thus, within the overlapping time periods (e.g. recent decades), the timescale (and corresponding time series) is generally consistent across different studies. These similarities provide further confidence in the results from the MEMD. Note that these studies have focused on different ENSO timescales and associated different spatial patterns during QB and LF/QQ ENSO events. However, IMFs are nonstationary and can thus capture changing patterns over time (Sect. 5.3). Also, previous studies have not necessarily considered relationships between different variables that are relevant for the ENSO dynamics (Sect. 6) and related conceptual oscillator models on specific timescales.

Overall, MEMD is consistent with other filtering methods (i.e. it acts like an effective band-pass filter) and can thus be used for further analysis of ENSO dynamics and its nonstationarity to determine if IMFs hold any physical significance. If IMFs have physical meaning, the information inferred from them can be valuable for enhancing climate models, for long-term predictions, for understanding teleconnections, and for exploring the underlying physics and variability of specific fields of interest, such as ENSO. In the following sections we address some of these aspects.

5.3 Nonstationarity of ENSO

The period ranges provided above (see also Fig. S1) suggest that there is some overlap between the period/frequency bands of IMF 11 and IMF 12, which is a result of nonlinear and nonstationary evolution of the modes (period is not constant as seen in Fig. 3). Indeed, there is a low (yet statistically significant) correlation between the two modes (∼ 0.28). Thus, in some time windows the two IMFs can describe the variability of a similar timescale (e.g. similar temporal evolution in Fig. 3 around years 1923–1927, 1951–1955, and 1997), but in other time windows they describe variability on different timescales. Despite changing periods, the two IMFs pass the red noise threshold in the Niño3 region and are therefore still considered significant and thus quasi-periodic (Sects. 4 and 5.1).

In some decades the band-pass-filtered Niño3 index (black solid line) is more consistent with the lower-frequency IMF 12 (red dashed line; approximately 1870–1917 and 1968–2000), and in other periods it is more consistent with the higher-frequency IMF 11 (red dotted line; approximately 1917–1968 and 2000–2010) (Fig. 3). This is consistent with the interdecadal shifts in the frequency of ENSO (Hu et al., 2017, 2020) that have occurred around years 1970 (from higher frequency to lower frequency) and 2000 (from lower frequency to higher frequency). Similar behaviour can also be seen earlier in the record. We can also see interdecadal changes in the amplitude of the ENSO modes, i.e. middle panel (period 1920–1965) in Fig. 3 compared with top and bottom panels (periods 1870–1920 and 1965–2010). This is somewhat consistent with Crespo et al. (2022), who found reduced amplitude of ENSO during the 1901–1931 and 1935–1965 periods relative to the post-1970 period.

Similarly, IMFs can capture different propagation directions of SST anomalies (Fig. 4). Previous studies have highlighted that SST anomaly propagation changed from westward to stationary or eastward around 1970 (Fedorov and Philander, 2000, 2001; Wang and An, 2001). Indeed, the two IMFs show westward propagation of SST anomalies prior to 1970. However, after 1970 it is slightly more complex. SST anomaly propagation becomes stationary or eastward in IMF 11 (except for a period just before 1990; Fig. 4a), but IMF 12 (Fig. 4b) still indicates either stationary or westward propagation with some rare (eastward) exceptions. This is then reflected in the sum of the two modes (Fig. 4c), which generally shows variability in propagation of SST anomalies post-1970. Furthermore, some of the characteristics of SST anomaly propagation from the sum of IMF 11 and IMF 12 carry over to 1-year low-passed SST anomalies (e.g. in the 1960s and 1980s; Fig. 4d).

ENSO events have also been characterised as east Pacific (EP) or central Pacific (CP) depending on the longitude where SST anomalies maximise (e.g. Kao and Yu, 2009; Singh and Delcroix, 2013; Zhang et al., 2019). The two IMFs are less able to reproduce this feature of ENSO diversity (Fig. 4). Most EP events are captured by IMF 11 and 12. However, some CP events (e.g. early 1990s; Fig. 4d) are not captured by quasi-periodic modes of variability (Fig. 4a, b, c). This suggests that other processes are likely causing them, such as sub-ENSO variability (Keenlyside et al., 2007). Similarly, some persistent La Niña events (e.g. mid 1970s, mid-1980s, and around year 2000; Fig. 4d) are not necessarily captured by the quasi-periodic modes of variability alone (potentially related to ENSO asymmetry; e.g. Choi et al., 2013; An et al., 2020). These results could also reflect that some aspects of ENSO nonlinearity are not represented by these two IMFs (e.g. Dommenget et al., 2013).

Figure 4 Time–longitude Hovmöller diagram of tropical Pacific SST anomalies (in K) averaged between 5^∘ S and 5^∘ N from (a) IMF 11, (b) IMF 12, (c) IMF 11 + IMF 12, and (d) 1-year low-passed (via fifth-order Butterworth filter) SST data. Black arrows help visualise the direction of propagation of SST anomalies in some periods (other periods are more stationary). We show the time period around year 1970 (1950–2000) where a change in propagation direction has been identified in previous work. Note that the colour scale in panels (a)–(b) is smaller than in panels (c)–(d).

Note that the magnitude of ENSO ultimately depends on all underlying modes of variability in the tropical Pacific (not just on the IMFs discussed here) – as also seen in Fig. 4. In fact, we find ∼ 5 modes (i.e. IMFs 10–14 here) with timescales ranging from ∼ 1 to ∼ 11 years (assessed via the 6.7th and 93.3rd percentiles as above) that can reproduce the majority of the ENSO variability (not shown), but the rest of these modes (i.e. IMF 10, 13, 14 here) are consistent with red noise (Fig. 2).

Figures 3 and 4 also show that weak Niño3 events have either small amplitudes (e.g. 1930–1940) of both IMFs or opposite amplitudes (e.g. 1985–1987). However, strong Niño3 events generally show a constructive interference or mode combination (e.g. 1997, a super-El Niño event), which is consistent with Slawinska and Giannakis (2017), Jajcay et al. (2018), Wang and Ren (2020), and Froyland et al. (2021).

This section has clearly demonstrated that MEMD together with a red noise test is suitable for identifying nonstationary quasi-periodic multivariate signals. This is a clear advantage of MEMD over other multivariate signal processing methods (e.g. PCA, MSSA). Below we now turn to ENSO dynamics to show that IMFs are also physical.

6 ENSO dynamics

The dynamics of ENSO typically involves positive (e.g. Bjerknes) and negative feedbacks between the atmosphere and ocean. The Bjerknes feedback (Bjerknes, 1969) refers to any decrease (increase) of trade winds that leads to reduced (enhanced) ocean upwelling (downwelling) and thus warming (cooling) in the eastern tropical Pacific, leading to reduced (enhanced) zonal SST and pressure gradients, which in turn reinforce the initial increase (decrease) of the trade winds. The most prominent and also the simplest negative feedback in the tropical Pacific involves recharge/discharge of ocean heat content via Sverdrup transport (e.g. Jin, 1997a, b; Burgers et al., 2005). Other negative feedbacks involve propagation and reflection of ocean Rossby and Kelvin waves (e.g. Suarez and Schopf, 1988; Battisti and Hirst, 1989; Picaut et al., 1997; Weisberg and Wang, 1997; Wang et al., 1999; Wang, 2001, 2018), where the latter can also be wind-forced. These processes generally involve changes in the thermocline depth, surface wind stress, and SST anomalies.

6.1 Quasi-oscillatory timescale

As ENSO dynamics is primarily related to the evolution of the ocean surface zonal wind stress (τ_x), thermocline depth, and the SST in the tropical Pacific (e.g. Wang, 2018), we demonstrate their relationship with phase composite maps of IMF 12 (Fig. 5). Shading in Fig. 5 represents SST anomalies, contours represent the thermocline depth anomalies (solid contours represent positive values and deeper thermocline), yellow arrows represent τ_x anomalies, and grey shading represents grid points of the SST that do not pass the red noise threshold. All values are standardised (i.e. divided by standard deviation).

The phase composites are computed using the instantaneous phase of the IMF 12's eastern Pacific SST (Niño3) time series (i.e. IMF 12 (SST (Niño3))) that we can obtain via Hilbert transform (Appendix B1, Eq. B4). This assigns every point in the eastern Pacific SST (Niño3) time series a phase between 0 and 360^∘, which can then be split into 12 phases (marked phase 0 through 11; e.g. phase 0 is 0–30^∘ and phase 11 is 330–360^∘), and all points in the time series (of 1-D or 3-D fields) belonging to a specific phase are then averaged to form map phase composites shown in Fig. 5 (see also Fig. S3 for IMF 11). Line phase composites in Fig. 6 (see also Fig. S4 for IMF 11) are constructed similarly. The eastern Pacific SST (Niño3) is chosen as an index here to focus on east Pacific (EP) ENSO events and related dynamics.

Figure 5 Latitude–longitude phase composite (phases 0 to 11 as labelled) of IMF 12: shading for SST, contours for thermocline depth (contour interval is the same as in the colour bar with solid contours representing positive values, and dashed contours represent negative values), and arrows for τ_x (the scale is shown in the bottom left corner of panels for phases  10 and 11). All data are standardised and all fields were composited based on the phase of the eastern Pacific SST (Niño3) index.

Figure 5 shows a typical cycle of EP ENSO in the tropical Pacific (on a ∼ 3-year timescale) as composited over the IMF 12's Niño3 index. This cycle can also be summarised with line phase composites (Fig. 6a, c) by averaging over specific regions of the tropical Pacific (as labelled; see also Table 1). Here, we analyse the eastern Pacific SST (Niño3), central Pacific τ_x (Niño4), western Pacific τ_x (Niño5), off-equatorial western Pacific thermocline depth (Niño6), and Pacific mean thermocline depth as they have historically been considered important for ENSO dynamics and have been used in ENSO conceptual oscillator models (e.g. Jin, 1997a; Burgers et al., 2005; Wang, 2018).

Figure 6 Phase composites of the eastern Pacific SST (Niño3) (black solid line), off-equatorial western Pacific isotherm/thermocline depth (Niño6) (black dashed line), Pacific mean isotherm/thermocline depth (blue dotted line), central Pacific τ_x (Niño4) (grey dashed line), and western Pacific τ_x (Niño5) (grey solid line). All fields are composited over the phases of the eastern Pacific SST (Niño3), such that they fit the phases in Fig. 5. (a) Composites of IMF 12 for data divided by the standard deviation of the corresponding time series (e.g. IMF 12 (thermocline)/σ (thermocline)); (b) composites of band-pass-filtered (29–53 months) standardised time series; (c) as in panel (a) but IMF time series are divided by IMF's standard deviation (i.e. IMF 12 (thermocline)/σ (IMF 12 of thermocline)).

Together the two figures (Figs. 5, 6a, c) suggest the following sequence of events:

1. during La Niña (phases 5–7) we have negative SST anomalies and a shallower thermocline in the Niño3 region, stronger easterly wind stress in the Niño4 region, and deeper thermocline in the western Pacific (including the Niño6 region);

2. as La Niña weakens (phases 8–10), the westerly wind stress in the Niño5 region and thermocline depth averaged across the tropical Pacific peak, starting the El Niño cycle;

3. as the SST warms, the eastern Pacific thermocline (Niño3) becomes deeper, the central Pacific wind stress (Niño4) becomes westerly, and the thermocline in the western Pacific (including the Niño6 region) becomes shallower (phases 11,0,1);

4. El Niño weakens (phases 2–4), the western Pacific τ_x (Niño5) becomes easterly, and the thermocline averaged across the Pacific becomes shallower, starting a La Niña event (phases 5–7);

5. the cycle repeats.

The evolution described above is also seen in the band-pass-filtered (29–53 months; 2.5–4.5 years) data (Fig. 6b). Note that the values in Fig. 6b are slightly larger than in Fig. 6a, because slightly different frequency ranges are ultimately represented in the two panels, but they remain qualitatively similar.

While Fig. 6 describes an average evolution of the variables considered, not all individual events have this exact behaviour. This is partly reflected in the reduced amplitudes of τ_x and thermocline depth in Fig. 6a compared with the SST amplitude. The full time series of these parameters from IMF 12 (Fig. 7; line types/colours are the same as in Fig. 6) show that the relationships from Fig. 6 occur often in the analysed period, especially for strong events. However, for weak events (middle panel in Fig. 7) the relationships are harder to establish – every event seems to be slightly different. This is somewhat consistent with Crespo et al. (2022), who noted that the dynamics of ENSO was different prior to 1970 relative to after 1970, with a dominant recharge–discharge oscillator (Jin, 1997a, b; Burgers et al., 2005) in the latter period.

Figure 7 IMF 12's standardised time series of the eastern Pacific SST (Niño3) (black solid line), off-equatorial western Pacific thermocline depth (Niño6) (black dashed line), Pacific mean thermocline depth (blue dotted line), central Pacific τ_x (Niño4) (grey dashed line), and western Pacific τ_x (Niño5) (grey solid line).

Similar results (phase composites) can also be obtained for the 16–33-month band-passed data and IMF 11 (Figs. S3, S4). This suggests that on average the QB and LF/QQ ENSO events have similar evolution and associated dynamics. However, the frequency of events (see full time series in Fig. S5) that follow the dynamics identified in the phase composites (Fig. S4) is lower in the IMF 11 case than in the IMF 12 case. This may be indicative of other processes that could be relevant for the QB ENSO.

The results presented here and summarised in the line phase composites (Figs. 6, S4) have some interesting implications for the dynamics of ENSO and conceptual oscillator models. The out-of-phase relationship between the Pacific mean thermocline depth (blue dotted line) and eastern Pacific SST (Niño3; black solid line) is a typical feature of the recharge–discharge oscillator (Jin, 1997a, b; Burgers et al., 2005). This suggests that IMFs that emerge from MEMD analysis can capture physical processes in the tropical Pacific. Furthermore, co-variability of the Pacific mean thermocline depth and the western Pacific wind stress (Niño5; grey solid line) suggests that wind forcing in the western Pacific may play a role in the recharge–discharge process or in the onset of ENSO in general as it precedes ENSO events (see also Roulston and Neelin, 2000; Capotondi et al., 2018; Lopez et al., 2013; Lopez and Kirtman, 2014; Timmermann et al., 2018).

On the other hand, the co-variability between the central Pacific wind stress (Niño4; grey dashed line), off-equatorial western Pacific thermocline depth (Niño6; black dashed line), and eastern Pacific SST (Niño3; black solid line) suggests that the unified oscillator proposed by Wang (2001) may need to be revised (see also Graham et al., 2015). This is because (i) in the unified oscillator model these three quantities should all be somewhat out of phase, but here we show that only the western Pacific wind stress and Pacific mean thermocline depth exhibit an out-of-phase relationship with the eastern Pacific SST (on average) and (ii) the Pacific mean thermocline depth does not co-vary with the off-equatorial western Pacific thermocline depth (Niño6), rendering the unified model's recharge–discharge oscillator simplification, which uses the off-equatorial western Pacific thermocline depth (Niño6), questionable.

That the central Pacific wind stress may be omitted in the unified model (due to co-variability with the eastern Pacific SST) was also mentioned in Graham et al. (2015). However, they suggested that western Pacific wind stress (Niño5) could also be omitted from the unified model, but Fig. 6 suggests that on 1.5–4.5-year timescale this variable should be kept. This may be because Graham et al. (2015) used 1-year low-pass-filtered data, which could have obscured the signal on a 1.5–4.5-year timescale (see Sect. 6.2). Additionally, Graham et al. (2015) suggested using the thermocline depth in the western Pacific averaged over a region that lies on the Equator rather than off-Equator. From Fig. 5 we can see that this would likely yield similar results to the off-equatorial western Pacific thermocline depth (Niño6) (to ±1 phase).

6.2 Other timescales

ENSO is a phenomenon that occurs on timescales of 2–8 years, and previous work has often used a 1-year low-pass filter to obtain ENSO. Thus, we now test the relationships between the eastern Pacific SST (Niño3), central Pacific τ_x (Niño4), western Pacific τ_x (Niño5), western Pacific off-equatorial thermocline depth (Niño6), and Pacific mean thermocline depth on slightly shorter and slightly longer timescales. We do this to test how relationships between different variables change across timescales that are still somewhat within the ENSO range.

Figure 8 shows phase composites similar to those in Fig. 6b but for band-passed data over a 12–19-month range (shorter periods; Fig. 8a) and over a 42–135-month range (longer periods; Fig. 8b). These timescale bands are consistent with IMF 10 (shorter periods) and IMFs 13 and 14 (longer periods), which are the modes that together with IMF 11 and 12 explain the majority of the variance in the eastern Pacific SST (Niño3).

Figure 8 As in Fig. 6b but for (a) higher-frequency (band-pass filter over 12–19 months) data and (b) lower-frequency (band-pass filter over 42–135 months) data. Note that higher- and lower-frequency timescale bands were chosen based on timescale ranges of the 10th and 13th–14th IMFs, i.e. the IMFs with slightly smaller or larger (respectively) timescales than IMF 11 and IMF 12.

Interestingly, on these shorter and longer timescales the evolution is different than on quasi-oscillatory timescales (1.5–4.5-year periods of IMF 11 and 12). Namely, the western Pacific wind stress (Niño5) closely follows the western Pacific off-equatorial thermocline depth (Niño6) (Fig. 8a, b). This suggests a very different role of the western Pacific wind stress on different timescales, which may be relevant for ENSO diversity and asymmetry. However, other variables remain similar across timescales (to ±1 phase). Thus, recharge–discharge processes operate throughout the ENSO timescale range (2–8 years), since the Pacific mean thermocline depth and eastern Pacific SST (Niño3) remain out of phase also on these shorter and longer timescales. Also, the results from these shorter and longer timescales indicate that the western Pacific wind stress (Niño5) may be omitted from conceptual oscillator models (on these timescales) as it is anti-correlated with the eastern Pacific SST (Niño3). This is consistent with Graham et al. (2015), who suggested this revision of the unified oscillator, but using 1-year low-pass-filtered data, which may have obscured the different behaviour on quasi-periodic timescales (Sect. 6.1; Fig. 6).

The above analysis shows that it is important to filter the data to correct frequency bands as there may be different behaviour present on different timescales, even within the ENSO range of 2–8 years.

7 Conclusions

In this study we have used a recently developed nonlinear and nonstationary method for identifying intrinsic variability of multivariate systems, the multivariate empirical mode decomposition (MEMD; Rehman and Mandic, 2010). The method can objectively identify modes of variability on different timescales within a nonlinear and nonstationary dataset describing a complex system such as the climate system. The timescale identification is objective as it is done without any pre-selection of a timescale window in which we expect the quasi-periodic behaviour. It finds a signal that is synchronised across input time series (here PC time series of combined fields, i.e. SST, wind stress, thermocline depth over the tropical Pacific) for every timescale within the given system. Also, the multivariate modes of variability that emerge from the MEMD analysis have nonstationary (i.e. evolving in time) patterns of variability (e.g. Fig. 4) – a clear advantage over some other multivariate time series analysis tools whose patterns are stationary.

Additionally, a red noise significance test has been developed to robustly identify quasi-periodic modes of variability in the given data, which had not been used before in the framework of MEMD. This means that MEMD can now be used as an alternative method for objective identification of the timescale of quasi-periodic motions in the climate system. Since the red noise test can be applied on every grid point separately, MEMD together with the red noise test can also be used for identifying potential new regions of quasi-periodic variability (similar to Fig. 5).

We demonstrate that MEMD can identify physical quasi-periodic modes of variability within the climate system by analysing the tropical Pacific SST variability. We have identified a clear quasi-periodic behaviour on a timescale of about 2–3 years (16–53 months) in the tropical Pacific. This timescale falls within the typical timescale range of ENSO, i.e. 2–8 years, and the dynamics of this quasi-periodic variability is consistent with ENSO dynamics. While ENSO quasi-periodic variability is well known, an identification (via MEMD) of a frequency range linked to the two dominant quasi-periodic modes of variability (i.e. 16–53 months) has still led to a few interesting results.

By analysing composites (e.g. Figs. 5, 6) of the thermocline depth, wind stress, and SST, we have shown that the ∼ 2–3-year (eastern Pacific) ENSO variability is generally consistent with the recharge–discharge conceptual oscillator model of ENSO (e.g. Jin, 1997a; Burgers et al., 2005). This oscillator describes an interplay between the Pacific mean thermocline depth and eastern Pacific SST (Niño3), which are related to recharge and discharge of ocean heat content on a timescale of about 3 years. However, the unified model (Wang, 2001) may need to be revised (see also Graham et al., 2015) as most of the variables relevant for the unified oscillator model co-vary. In particular, the western Pacific off-equatorial thermocline depth (Niño6), central Pacific wind stress (Niño4), and eastern Pacific SST (Niño3) exhibit the same phase relation and thus describe the same dynamics. As also the western Pacific wind stress (Niño5; relevant for the unified oscillator) and Pacific mean thermocline depth co-vary, there seems to be a close relationship between the recharge–discharge and unified oscillator models, but it may be different than previously thought.

On shorter and longer timescales (12–19 months and 42–135 months; Fig. 8) the relationships between variables are different; in particular, the western Pacific wind stress (Niño5) on these timescales co-varies with the eastern Pacific SST (Niño3). This suggests that the role of the western Pacific atmosphere–ocean processes (via wind stress) in tropical Pacific variability (and ENSO) can be very different on different timescales, whereas other variables (and in other regions) largely keep their relationships across timescales. This means that the recharge–discharge oscillator model operates on all timescales considered (12–135 months), even if it has a characteristic timescale. However, the relations underlying the unified oscillator model exhibit different behaviour on different timescales, the implications of which should be explored further in the future.

MEMD analysis could be extended in several ways, for example (i) to assess ENSO dynamics in models as they typically struggle with the representation of the western Pacific processes (i.e. relevant for Niño5 wind stress; Guilyardi et al., 2020; Planton et al., 2021) and (ii) to study ENSO teleconnections on different timescales (e.g. Brönnimann, 2007; Fereday et al., 2008; Jiménez-Esteve and Domeisen, 2018, 2020; Hardiman et al., 2019). Also, future studies should involve an examination of sensitivity and causal links (not established here) between different variables and their links across scales (e.g. Runge et al., 2015; Jajcay et al., 2018; Jenney et al., 2019; Kretschmer et al., 2021), as well as dedicated model experiments.

Overall, this study has analysed the variability in the tropical Pacific (using MEMD with a red noise test), identified two quasi-periodic modes of variability (on ∼ 2–3-year timescale), related their physics to the recharge–discharge oscillator, and suggested a revision to the unified oscillator model (somewhat consistent with Graham et al., 2015). As the variability on this timescale is quasi-periodic, it may be predictable far in advance, which calls for further investigations of the tropical Pacific variability and related teleconnections, for their prediction, and for further model improvements (see also Chen et al., 2021; Lee et al., 2021).

Appendix A: MEMD for 3-D fields

To find the intrinsic variability of our 3-D field, i.e. ${\mathbf{A}}^{\prime }(L_t,L_y,L_x)$ mentioned in Sect. 3.3, we first reduce the dimensionality of our data by decomposing it using the singular value decomposition (SVD), which yields spatial patterns of our data (empirical orthogonal functions, EOFs) and corresponding time series (principal components, PCs). First, we multiply A^′ by $\sqrt{\cos \mathit{\varphi }}$ (area weighting; with ϕ latitude), divide by standard deviation (σ) at each grid point (and for each variable separately, if relevant), and reshape A^′ from a 3-D array $(L_t,L_y,L_x)$ to a 2-D array $(L_t,L_x\cdot L_y)$. Then A^′ can be expressed with a singular value decomposition as

$$\begin{array}{}\text{(A1)}& {\mathit{\sigma }}^{-\mathrm{1}}(L_x\cdot L_y){\mathbf{A}}^{\prime }(L_t,L_x\cdot L_y)=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T,\end{array}$$

where U and V represent the left and right singular vectors related to PCs and EOFs, $\mathbf{\Sigma }=\sqrt{(L_t-\mathrm{1})\mathbf{\Lambda }}$ is a diagonal matrix with square roots of the explained variance of each mode (denoted Λ, i.e. eigenvalues) along the diagonal, L_t is the length of time series, and superscript T denotes transpose. PCs are defined as $\mathbf{U}\sqrt{L_t-\mathrm{1}}$ and EOFs as $\mathbf{\Sigma }{\mathbf{V}}^T/\sqrt{L_t-\mathrm{1}}$, such that ${\mathit{\sigma }}^{-\mathrm{1}}(L_x\cdot L_y){\mathbf{A}}^{\prime }$ can be represented as a function of PCs and EOFs:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& {\mathit{\sigma }}^{-\mathrm{1}}(L_x\cdot L_y){\mathbf{A}}^{\prime }(L_t,L_x\cdot L_y)\\ & \simeq \sum _{m=\mathrm{1}}^{m=\mathrm{20}}\mathrm{EOF}(m,L_x\cdot L_y)\mathrm{PC}(m,L_t),\end{array}\end{array}$$

where m corresponds to the PC number and is ordered according to the eigenvalues (m=1 for the largest eigenvalue). We retain only the leading 20 PCs for the analysis (they generally describe the majority of the variance in A^′).

Now we can use the 20 PCs (PC(m,L_t)) as input to the MEMD algorithm (for details on the algorithm itself see Rehman and Mandic, 2010). This algorithm finds common timescales (i.e. intrinsic mode functions, IMFs) within the 20 PCs and splits each PC into several IMFs (the number of IMFs is not predetermined). Thus, each PC can be represented as a sum of IMFs

$$\begin{array}{}\text{(A3)}& \mathrm{PC}(m,L_t)=\sum _{s=\mathrm{1}}^{s=s_{\max }}\mathrm{IMF}(s,m,L_t),\end{array}$$

where s corresponds to IMF number and is ordered according to the timescale (s=1 for the shortest timescale and s_max for the longest timescale, which is usually a trend or a residual). Equation (A3) shows that each PC is a superposition of different IMFs (see also Table S2) and with it also a superposition of modes of variability in the selected field(s) with different timescales.

Each PC is associated with a spatial pattern (EOF$(m,L_x\cdot L_y)$), which allows a reconstruction of the time–space $(L_t,L_x\cdot L_y)$ structure/evolution for each IMF, yielding IMFs of the initial dataset A^′ (IMFA). We do this by multiplying PCs for each IMF with corresponding EOFs and summing over all 20 PCs/EOFs (similar to Eq. A2):

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\mathit{\sigma }}^{-\mathrm{1}}(L_x\cdot L_y)\mathrm{IMFA}(s,L_t,L_x\cdot L_y)\\ & \simeq \sum _{m=\mathrm{1}}^{m=\mathrm{20}}\mathrm{IMF}(s,m,L_t)\mathrm{EOF}(m,L_x\cdot L_y).\end{array}\end{array}$$

Here note that to get IMFA in the units of the input field we must multiply it by the field's standard deviation as the input data for the SVD algorithm were standardised (Eq. A1). Again, IMFAs are ordered by timescale, i.e. IMFA1 with the shortest timescale and IMFAs_max with the longest timescale (trend). From here we can reconstruct A^′ by summing over all IMFAs,

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& {\mathit{\sigma }}^{-\mathrm{1}}(L_x\cdot L_y){\mathbf{A}}^{\prime }(L_t,L_x\cdot L_y)\\ & \simeq \sum _{s=\mathrm{1}}^{s=s_{\max }}\mathrm{IMFA}(s,L_t,L_x\cdot L_y),\end{array}\end{array}$$

and ultimately one can also reshape A^′ from a 2-D array $(L_t,L_x\cdot L_y)$ to a 3-D array $(L_t,L_x,L_y)$. The importance of each IMFA for A^′ can be assessed by computing the explained variance of each IMFA or other significance methods. To find IMFA modes (and grid points) that correspond to potentially oscillatory behaviour we must perform a red noise test (see Appendix B).

Note that from here on (and in the main text) IMFAs are referred to as IMFs for simplicity.

Appendix B: Significance tests

Typically we can test if the modes (IMFs) are different from white or red noise, depending on the distribution of our data. In the climate system, variables often exhibit behaviour that resembles white or red noise. The IMFs that are significant (i.e. different from both red and white noise) likely represent quasi-oscillations, indicating higher potential for predictability of processes that correspond to the timescale of the significant IMF. Thus, this distinction is very important in climate science. Therefore, we discuss the white and red noise tests (for 1-D data, i.e. time series) below, whereas the robustness of IMFs from the MEMD analysis and the relevant significance tests are discussed in the main text (Sects. 4 and 5).

Note that the white and red noise tests are performed on 1-D time series; hence, EMD (univariate decomposition; see main text) is first used to test the performance of IMFs that arise from the EMD analysis. The multivariate data (via MEMD) in the main text are analysed with the simplest and most relevant test (i.e. theoretical red noise test described below).

B1 White noise test

The white noise significance test has been derived by Wu and Huang (2004), who showed that the energy density function of the sth IMF (E_s; i.e. average squared amplitude of the sth IMF) is

$$\begin{array}{}\text{(B1)}& E_s={\displaystyle \frac{\mathrm{1}}{L_t}}\sum _{j=\mathrm{1}}^{L_t}\left[{\mathrm{IMF}}_s\right(j){]}^{\mathrm{2}},\end{array}$$

where IMF_s(j) is the sth IMF at time step j ($j=\mathrm{1},\mathrm{\dots },L_t$), and L_t is the length of the time series. Wu and Huang (2004) further showed that the total energy density of the sth IMF can then be expressed as

$$\begin{array}{}\text{(B2)}& L_t{E}_s=\int S(\mathit{\nu }{)}_s\mathrm{d}\mathit{\nu },\end{array}$$

where ν is frequency, and S(ν)_s is the power spectrum of the sth IMF. From this they showed that for white noise,

$$\begin{array}{}\text{(B3)}& \ln E_s\approx -\ln T_s,\end{array}$$

where T_s is the average period of the sth IMF, and ln denotes the natural logarithm.

Note that the frequency (and thus also period) of each IMF is computed using Hilbert transform by first generating an analytical signal (e.g. Huang et al., 1998):

$$\begin{array}{}\text{(B4)}& Z\left(t\right)=X\left(t\right)+\mathrm{i}Y\left(t\right)=\left|Z\right(t\left)\right|e^{i\mathit{\theta }\left(t\right)},\end{array}$$

where t is time dimension, X(t) is our IMF time series, Y(t) is its Hilbert transform, Z(t) is the analytical signal, and $\mathit{\theta }\left(t\right)=\arctan \left(Y\right(t)/X(t\left)\right)$ is the instantaneous phase. The instantaneous frequency can then be computed by taking a time derivative of the phase, i.e. $\mathit{\nu }=\mathrm{d}\mathit{\theta }\left(t\right)/\mathrm{d}t/\mathrm{2}\mathit{\pi }$, and the average frequency of each IMF is computed by averaging the instantaneous frequency in time (note that period = $\mathrm{1}/\mathit{\nu }$, i.e. T_s).

The relationship between the logarithms of energy density and average period of the IMFs (Eq. B3) is then used in Fig. B1a (black solid line) to test whether an IMF (using EMD of the Niño3 index; blue dots in Fig. B1a; see also Sect. 3) is different from white noise or not. The mode is significant with respect to white noise if it exceeds the one-tailed 95th percentile threshold (denoted by the black dotted line). The percentile range serves as a significance test; i.e. if IMFs from our data are above, for example, the 95th percentile they are significant at the 95th percentile level. The percentile range can be expressed analytically as (Wu and Huang, 2004)

$$\begin{array}{}\text{(B5)}& \ln E_s=-\ln T_s\pm p\sqrt{{\displaystyle \frac{\mathrm{2}}{L_t}}}\exp (\ln T_s/\mathrm{2}),\end{array}$$

where p denotes a threshold (p=1.645 for the one-tailed 95th percentile of the Gaussian distribution). Note that typically the number of degrees of freedom (DoFs) for white noise data is expected to be equal to the total energy density of the sth IMF (i.e. DoF_s=L_tE_s; Wu and Huang, 2004).

Figure B1 Significance tests for EMD modes: (a) the white noise significance test and (b) the red noise significance test for EMD IMFs of the Niño3 index (blue dots). Blue dots are computed as in Fig. 2. (a) The black solid line represents the theoretical linear relationship between the logarithm of the average period (T_s) and the logarithm of energy density (E_s; i.e. average squared amplitude; Eq. B3); the black dotted line represents 5th–95th percentile (Eq. B5); the grey dots represent I=L_t realisations of IMFs of white noise time series (the length is the same as for the Niño3 index), whereas the grey solid and dotted lines represent their mean and the 5th–95th percentile, respectively. (b) The red solid line represents the theoretical red spectrum energy density (Eqs. B7–B9); the red dotted line represents the 95th percentile (via χ² test); the light pink dots represent I=L_t realisations of IMFs of the red noise time series (Eq. B6; the length L_t is the same as for the Niño3 index), whereas the pink solid and dotted lines represent their mean and the 95th percentile, respectively. Note that the x axis shows the logarithms of the period (log _e(T_s)) ordered from the shortest period (highest frequency) to the longest period (trend). For further descriptions of the figure see text.

Alternatively, we can test whether the input data are different from white noise by constructing multiple (I) realisations of synthetic white noise time series w_i (ith random normally distributed time series with standard deviation σ of 1). Then we can compute its IMFs via EMD (Sect. 3), and we can repeat the process I times. Employing Eq. B1 on these IMFs yields scattered grey dots in Fig. B1a (constructed in the same way as blue dots; see also Sect. 4, Fig. 2), where their mean and 5th–95th percentile are shown as grey solid and dotted lines, respectively.

A comparison with the IMFs from the input data (Niño3 index; blue dots in Fig. B1a) reveals that many IMFs lie outside the white noise range and that overall the data (blue dots) distribution does not resemble the white noise (grey dots) distribution (not noted in Wu and Huang, 2004). This suggests that a white noise test for such data is not a good test. Indeed, atmosphere–ocean coupled systems, such as ENSO, can often be represented as a red noise process (e.g. Hasselmann, 1976; Frankignoul and Hasselmann, 1977); thus, we now turn to a similar test but for data distributed as red noise.

B2 Red noise test

B2.1 Synthetic red noise data

To test if our data (e.g. Niño3 index) is purely red noise or it has inherent oscillations that we can identify through the (M)EMD analysis, we generate I realisations of synthetic red noise time series x (AR(1) process) as (e.g. Gilman et al., 1963)

$$\begin{array}{}\text{(B6)}& x_{i,j+\mathrm{1}}=r{x}_{i,j}+\sqrt{\mathrm{1}-r^{\mathrm{2}}}{w}_{i,j+\mathrm{1}},\end{array}$$

where r is the lag-1 auto-correlation from our data (e.g. Niño3 index), w is white noise (as in Appendix B1), i runs over I realisations of synthetic red noise data, and j (j>1 and j≤L_t; with L_t length of our data, e.g. the length of Niño3 record) is an index that runs over the time steps (one time step is one time unit, e.g. 1 month). For j=1 (the first time step) we set $x_{i,\mathrm{1}}=w_{i,\mathrm{1}}$.

Once we obtain the red noise time series x_i we can compute its IMFs via EMD (Sect. 3) and we can repeat the process I times (as for the white noise; Appendix B1). This yields the pink scattered dots in Fig. B1b. The mean over I cases for each IMF (frequency band) is shown by the pink solid line, and the (one-tailed) 95th percentile across the I cases is shown by the pink dotted line. Note that we plot logarithmic values in Fig. B1, as mentioned above.

Note that Franzke (2009) used a similar approach for indices such as the North Atlantic Oscillation and found a simple relationship between the power spectrum and frequency, consistent with Kolotkov et al. (2016). However, we follow Gilman et al. (1963) to define a relationship between the power spectrum and frequency. This incorporates the lag-1 auto-correlation of the time series into the theoretical red noise power spectrum (see below).

B2.2 Theoretical red noise test

Alternatively, one can compute a theoretical power spectrum of the red noise (see Gilman et al., 1963):

$$\begin{array}{}\text{(B7)}& S\left(\mathit{\nu }\right)={\displaystyle \frac{\mathrm{1}-r^{\mathrm{2}}}{\mathrm{1}-\mathrm{2}r\cos \mathrm{2}\mathit{\pi }\mathit{\nu }+r^{\mathrm{2}}}},\end{array}$$

where S is the power spectrum of red noise, $\mathit{\nu }=\mathrm{1}/t$ is the frequency, and r is again the lag-1 auto-correlation from our data. For each frequency estimate we must multiply S(ν) by the frequency range (Δν) (as in Eq. B2) to obtain a theoretical estimate of the (mean) energy of the red noise (E^red) (cf. Kolotkov et al., 2016):

$$\begin{array}{}\text{(B8)}& E^{\mathrm{red}}\left(\mathit{\nu }\right)=S\left(\mathit{\nu }\right)(\mathit{\nu }\mathit{\alpha }-\mathit{\nu }/\mathit{\beta }),\end{array}$$

where $\mathit{\beta }=\sqrt{{\mathit{\nu }}_s/{\mathit{\nu }}_{s+\mathrm{1}}}$ and $\mathit{\alpha }=\sqrt{{\mathit{\nu }}_{s-\mathrm{1}}/{\mathit{\nu }}_s}$ with s running over frequencies (from higher to lower frequency). Note that since EMD is a dyadic filter (each lower frequency is a half of the previous one; e.g. Flandrin et al., 2004; Rehman and Mandic, 2011) both α and β typically take a value of $\sqrt{\mathrm{2}}$ (consistent with Kolotkov et al., 2016). However, when mode mixing is present (e.g. here it is generally present at higher and lower frequencies of IMFs from MEMD analysis) this is not necessarily true, hence the use of α and β in Eq. (B8). Finally, E^red for every sth IMF (i.e. $E_s^{\mathrm{red}}$) must be scaled such that its total energy is the same as the total energy of our data (e.g. Madden and Julian, 1971; Bretherton et al., 1999):

$$\begin{array}{}\text{(B9)}& E_s^{\mathrm{red}}=E_s^{\mathrm{red}}{\displaystyle \frac{{\sum }_{s=\mathrm{2}}^{s_{\max }-\mathrm{1}}{E}_s}{{\sum }_{s=\mathrm{2}}^{s_{\max }-\mathrm{1}}{E}_s^{\mathrm{red}}}},\end{array}$$

where s_max is the number of IMFs (as above), and s represents the sth IMF of frequency ν_s; E_s was defined previously (Eq. B1). Note that we scale $E_s^{\mathrm{red}}$ from total energies of IMFs between s=2 and $s=s_{\max }-\mathrm{1}$ as the last IMF is typically a trend/residual and the first IMF does not necessarily follow the distribution correctly (but including the two usually does not significantly alter the results). $E_s^{\mathrm{red}}$ is shown in Fig. B1b as a red solid line.

This red noise test is typically used in climate science to determine the significance of power spectra peaks in our data (using S(ν) from Eq. B7), and it differs from the red noise test of Kolotkov et al. (2016) as it takes into consideration the lag-1 auto-correlation of the data. If the cosine function in the S(ν) (Eq. B7) is expanded into the Taylor series ($\cos \mathrm{2}\mathit{\pi }\mathit{\nu }\approx \mathrm{1}-(\mathrm{2}\mathit{\pi }\mathit{\nu }{)}^{\mathrm{2}}/\mathrm{2}+..$) one can realise that for large ν (high frequencies) S(ν) indeed reduces to the spectrum γν⁻² (with γ a constant) suggested by Kolotkov et al. (2016); Franzke (2009). However, for low frequencies (small ν) they do not agree well and ultimately S(ν) also becomes a constant (see Fig. B2 for comparison). Furthermore, as S(ν) depends on the lag-1 auto-correlation (r), we can see from Eq. (B7) that for r=0, S(ν)=1; i.e. it reduces to the power spectrum of the white noise. This means that this theoretical test can potentially be used for testing the significance of the data that correspond to either white or red noise.

Figure B2 Red noise power spectrum (S(ν)) for (black line) $S\left(\mathit{\nu }\right)=\mathit{\gamma }{\mathit{\nu }}^{-\mathrm{2}}$ and (grey line) S(ν,r) from Eq. (B7) for r=0.9. γ was estimated as a ratio between the integrated power spectra of the two spectra for frequencies higher than 0.02 per month ($\mathit{\gamma }={\sum }_{\mathit{\nu }}{\mathit{\nu }}^{-\mathrm{2}}/{\sum }_{\mathit{\nu }}S(\mathit{\nu },r)$) where the two power spectra generally agree well.

The significance of the IMFs from the input data is tested using a χ² test, where the sth IMF's ${\mathit{\chi }}_s^{\mathrm{2}}$ value for the (one-tailed) 95th percentile is computed from DoF_s=L_effE_s degrees of freedom (instead of L_tE_s as was the case for white noise, due to strong correlations between neighbouring data points; Bretherton et al., 1999; Wu and Huang, 2004; Kolotkov et al., 2016), where

$$\begin{array}{}\text{(B10)}& L_{\mathrm{eff}}={\displaystyle \frac{\mathrm{1}-r^{\mathrm{2}}}{\mathrm{1}+r^{\mathrm{2}}}}{L}_t.\end{array}$$

Then we multiply the expected red noise curve $E_s^{\mathrm{red}}$ by ${\mathit{\chi }}_s^{\mathrm{2}}/{\mathrm{DoF}}_s$ (e.g. Madden and Julian, 1971; Bretherton et al., 1999) to ultimately obtain a threshold for the 95th percentile (red dotted line in Fig. B1b). Note that for DoF_s<1 we set DoF_s=1 (to avoid numerical issues). The IMFs derived from the data (blue dots in Fig. B1b) that exceed the red noise threshold (i.e. lie above the red/pink dotted line in Fig. B1b) are considered significant at the 95th percentile (one-tailed).

Figure B1b shows that the two (synthetic and theoretical) red noise tests (for the Niño3 index via EMD) are somewhat comparable and that the majority of the input (e.g. Niño3 index) data (blue dots) lie within the red noise range (i.e. within the pink dots and below the pink/red dotted line). However, we can identify one IMF (period ∼ 31 months or ∼ 2.5 years) that is above the red noise threshold and well within the typical ENSO timescale (2–8 years), suggesting quasi-periodic behaviour (oscillations).

Similarly, we can use the theoretical red noise significance test on IMFs obtained via MEMD, which yields two significant modes of variability in the eastern Pacific SST (Niño3) index (see the main text). Note, however, that we do not necessarily expect exactly the same results from the EMD (Fig. B1b) and MEMD (Fig. 2) methods since the MEMD finds a synchronised signal within the tropical Pacific and across different variables, whereas EMD only analyses the 1-D Niño3 time series. This is also true for the number of IMFs obtained via the two different methods. MEMD yields significantly more IMFs as EMD (21 versus 10), which is likely a result of inputting several different time series with different timescales, especially in the high- and low-frequency range (i.e. periods shorter than about 8 months and longer than about 700 months).

Alternatively, one can also compute significant modes by computing the red noise test at each grid point and then average the results over all grid points, but we have not used this here. Instead, we use an additional test on map plots in Sect. 5 (Figs. 5 and S3), where we identify potentially oscillatory grid points and use grey shading on areas that are well represented with red noise alone (i.e. not significant).

Despite some differences between the MEMD and EMD modes of variability (primarily due to different input time series), the two methods agree on the quasi-periodic timescale of 2–3 years in the Niño3 region. This is also consistent with the significant periods inferred from the usual 1-D wavelet transform (not shown) and the power spectrum analysis of the Niño3 index (1-D) obtained via Fourier transform (Fig. B3). The latter confirms that the 2–3-year timescale is quasi-periodic as there is a significant peak in the power spectrum of the Niño3 index on those timescales; i.e. the black solid line in Fig. B3 (Niño3 index) is above the red dashed line (red noise one-tailed 95 % threshold).

Figure B3 Power spectrum of the Niño3 index. The power spectra are first computed for 500-month-long chunks (overlapped by 250 months) and then averaged over all cases (grey solid line). The black solid line represents a 10-point running mean of the black dotted line (to increase the number of degrees of freedom, which is $f_{\mathit{\omega }}{L}_t/\mathrm{0.5}{L}_{\mathrm{chunk}}=\mathrm{10}\times \mathrm{1680}/\mathrm{250}\approx \mathrm{67}$; see also Boljka et al., 2021). The red solid and dashed lines represent the theoretical red noise test and its (one-tailed) 95th percentile, respectively.