Interactions in the observed ENSO

We examined synchronization and causal interactions in the NINO3.4 time series for the quasi-oscillatory modes with periods ranging from 5 to 96 months. We are looking for the observed interactions characterized by causality estimate exceeding the 95th percentile of the distribution of this quantity for the surrogate data samples (Fig. 1), where the surrogate time series were generated using a Monte Carlo method,24 yielding synthetic time series with the same spectrum, but void of any cross-scale interactions (see Methods for further details).

Fig. 1 Cross-scale phase synchronization (a), phase–phase causality (b), and phase–amplitude causality (c) in the observed NINO3.4 time series. The phase synchronization is a symmetrical relation, hence the plot is symmetric, while causality plots are shown with the period of the driver (master) time series on the x-axis and driven (slave) time series on the y-axis. Shown are (positive) significance-level deviations from the 95th percentile of the k-nearest neighbor estimates of (conditional) mutual information, tested using 500 Fourier transform surrogates (see Methods) Full size image

There are three pairs of modes that exhibit phase synchronization (a process by which two or more cyclic signals tend to oscillate with a repeating sequence of relative phase angles) (Fig. 1a). The quasi-biennial (QB) modes (periods of 1.8–2.1 year), which tend to peak in winter (not shown), are indeed synchronized with the annual cycle (AC) as well as with the so-called combination tones (CT; periods approximately 9 and 14 months). The combination tones result from the interaction of AC with low-frequency (LF) ENSO modes of periods 4–6 yr (see, for example, Stein et al.25 or Stuecker et al.26). Furthermore, the AC and CT modes with periods of 8–9 months are themselves phase locked. Finally, the LF modes with periods 5–6 yr and QB modes with periods 2–3 year also exhibit phase synchronization. These results reconfirm an important role of the annual cycle in ENSO dynamics, with strong ENSO events peaking in boreal winter,10,27 and point to the link between QB and LF modes, which may be responsible for extreme ENSO events14,15,16 (see below); thus, our synchronization analysis brings out known ENSO properties consistent with previous research.

On the other hand, the phase–phase causality diagnostics (Fig. 1b) provides an additional—this time new—information that complements the phase synchronization results and addresses the causes of these synchronizations by elucidating important directed connections between the LF, QB, and AC/CT ENSO modes. This analysis identifies a member of the pair of modes (time series) that is a skillful precursor of the other member in predictive context. For example, the phase of LF modes affects that of the AC/CT modes, which means that the “shape” of the annual cycle depends on whether the LF mode (periods of 4–6 year) is in its extreme warm or extreme cold phase.

Furthermore, the phase of the AC mode is a skillful precursor for the phase of QB modes with the periods of 1.8–2.1 year, while the phase of QB modes with the periods of 2–3-year dictates in part the phase of the CT modes (periods of 12–16 months).

The only pronounced phase–amplitude causality link (Fig. 1c) is the one between the phase of the LF ENSO mode (periods of 5–6 year) and the amplitude of QB modes (periods of 1.8–2.1 year).

Our analysis thus identifies the three fundamental time scales in ENSO dynamics—AC, QB, and LF—consistent with previous work,15,16,17 but offers further details on the interaction between these modes. Some of the interactions we identified rigorously here have been previously theorized to exist, but, to the best of our knowledge, were never detected in a data-driven way. Based on our results, it is natural to consider the AC and CT processes in combination to define the quasi-annual (QA) variability. The QB modes can be divided into two—‘‘faster’’ and ‘‘slower’’—sub-ranges, with the periods of 1.8–2.1 year and 2–3 year, respectively. Similarly, the LF processes can be divided into the ones associated with 4–5-year and 5–6-year periods.

We observe a pronounced connection between the (phase of) the slower LF mode and both the phase and amplitude of the faster QB mode. In particular, the slower LF mode affects the phase of the QA mode, and, therefore—indirectly—the phase of the faster QB mode, which tends to be affected by and phase-synchronized with the QA mode; the slow LF mode also directly affects the amplitude of the faster QB mode. The connections between the phase of slow LF mode and the phase of QA mode important in the causal sequence above are both direct and indirect. In the latter indirect case, the connection works through the phase synchronization between the slow LF mode and the slow QB mode and subsequent causal effect of the latter on the phase of the QA mode. The faster LF modes add to the picture by also affecting the phase of the QA mode, and, therefore, indirectly, the phase of the faster QB mode.

Consequences of causal connections

One of the main findings of our study is the apparent importance, in ENSO dynamics, of the LF phase → QA phase → QB phase causal linkages, as well as LF phase → QB amplitude causal linkages. To illustrate these causal connections further, we utilized an approach of conditional composites, in which we first identified three distinct phases (by dividing the span between maximum and minimum values into three bins) of the low-frequency (LF) ENSO component: LF−, LF0 and LF+; and subsequently computed the composites (mean values) of any variable of interest over the data points associated with these three phases.

Figure 2 visualizes the response of the AC frequency to changes in the phase of the LF ENSO variability, thus illustrating the LF phase → QA phase causal linkage. Here, we composited the instantaneous frequencies of the annual cycle computed as the slope of the continuous 12-month-long snippets of the AC-phase time series via the robust regression. These results are presented in the form of histograms and suggest that the positive phase of the LF cycle speeds up the annual cycle (thus shortens its period and increases its frequency), while the negative period of LF cycle causes the AC period to become longer than a year. Not surprisingly, the annual cycles associated with neutral LF0 conditions, as well as climatological annual cycles, have the average period of exactly 1 year.

Fig. 2 Histograms of instantaneous frequencies of the NINO3.4 annual cycle. Top: unconditioned (using all data); bottom: histograms conditioned on the phase of the low-frequency (LF) ENSO mode. The black horizontal line marks exactly 1 year period, and also given are respective means ± one standard deviation of instantaneous frequencies for all four panels Full size image

In fact, not only the effective frequency, but also the entire shape of the annual cycle changes depending on the phase of the LF ENSO mode. Figure 3 shows conditional composites of the annual cycle associated with the LF−, LF0, and LF+ phases of the low-frequency ENSO cycle; these composites were computed by averaging the raw NINO3.4 data associated with a given phase of LF variability for each month. The neutral LF conditions correspond with the seasonal cycle of NINO3.4 temperatures that is close to climatological seasonal cycle. The latter cycle is not purely harmonic, and is characterized by relatively fast warming between January and May and a slower cooling afterwards. The annual cycle conditioned on LF− phase has the same general character, but is on average colder than the climatological AC. The annual cycle associated with the LF+ phase of interannual ENSO signal is, of course, warmer on average due to El Niño-type conditions, but also has a very different, more harmonic shape, with September-through-December warming absent from the LF−, LF0, and climatological annual cycles. Also shown in Fig. 3 are the NINO3.4 time series during two extreme El Niño events (namely, 1982/83 and 1997/98 events). Note how a very strong biennial ENSO signal during those years completely masks any visible AC variability that may be present in the NINO3.4 time series at that time. An episodic character of such pronounced biennial extreme events makes it difficult to associate the QB ENSO variability with alterations between weak and strong annual cycles, as was suggested previously.12 We argued that such extreme ENSO events arise due to synchronization between a suite of different QB modes, which are individually characterized by a relatively small variance.

Fig. 3 Annual cycle composites of NINO3.4 index conditioned on the phase of the low-frequency (LF) ENSO mode. Also shown are snippets of the NINO3.4 time series during two particular extreme ENSO events: 1982/83 and 1997/98 Full size image

We hypothesize here that these ‘‘internal’’ QB synchronizations arise due to causal interactions represented by the QA phase → QB phase causal linkage identified by our conditional mutual information analysis. This linkage is further illustrated, albeit indirectly, in Fig. S2 in the Supplementary Material, which shows that the composite annual cycles associated with QB–, QB0, and QB+ phases of QB variability closely match those associated with LF−, LF0, and LF+ phases, respectively, consistent with LF phase → QA phase → QB phase-directed connections.

Moreover, Fig. S3 in the Supplementary Material identifies clear growth in the amplitude of the BC and QB variability as the low-frequency phase changes from La Niña (LF−) to El Niño (LF+) conditions, consistent with the LF phase → QB amplitude causal interaction. On the other hand, the changes in the amplitude of AC conditioned on the LF phases are non-monotonic and somewhat less pronounced compared to those in BC and QB amplitudes.

The interactions identified above are instrumental in setting up extreme ENSO events (Fig. 4). In particular, during all of the strong El Niño events of years 72/73, 82/83, and 97/98, the QA, QB, and LF modes were characterized by synchronous pronounced maxima. Note, however, that strictly speaking, the synchronization with LF mode is not really necessary for an extreme ENSO event to materialize, since the ‘‘peak’’ of this mode spans a good part of the year (for example, the peak of LF mode did not really occur in winter 1982/83, but the LF wintertime ‘‘background’’ during that time was still abnormally warm), while the amplitude of the QA mode is, in general, small, so that this mode does not contribute much directly to the magnitude of a given event. Instead, what appears to be essential for an extreme ENSO to occur is the synchronization of multiple QB modes with each other. We believe that this ‘‘internal’’ QB synchronization is what has been picked up by our conditional mutual information analysis in the form of LF → QA → QB phase connections and also LF phase → QB amplitude connections (since synchronization of phases of different QB modes should automatically result in a large-amplitude event.) By contrast, during a moderate El Niño of 87/88, the LF, QB, and QA modes exhibited phase shifts, with lower-frequency modes leading the higher-frequency modes (in particular a suite of QB modes) instead of being ‘‘stacked’’ on top of one another, thus limiting the magnitude of this event. Notably, strong La Niña events do not seem to be associated with the minimum of the LF mode, but instead occur during near-neutral LF conditions when the minima of the QA modes and the minima of the whole range of QB modes synchronize. Thus, in both El Niño and La Niña cases, the behavior of the QB modes has a vital control on the magnitude of the ENSO events.

Fig. 4 Top: wavelet reconstructions of the observed NINO3.4 index (1970–1999): reconstruction of the annual (black), quasi-biennial (for a range of periods from 18 to 30 months, with 2 month step; shades of blue to green) and 5-year ENSO cycles (red). All of these reconstructions were computed via continuous complex wavelet transform (CCWT) as A p (t) cos ϕ p (t), for the corresponding central wavelet periods p. Bottom: the full observed normalized NINO3.4 index. The years of strong El Niño and La Niña events are marked with the red and blue shading, respectively Full size image

Interactions in CMIP5 models

The Coupled Model Intercomparison Project Phase 5 (CMIP5)28 is a framework for coordinated climate change experiments providing global circulation model (GCM) outputs from various modeling groups. We analyzed time series of the sea-surface temperature obtained from the individual runs of CMIP5 models and compared the simulated NINO3.4 characteristics with the observed characteristics (Fig. 5). To start with, we measured the similarity between the observed and simulated wavelet spectra29 using root-mean-square distance and Pearson correlation coefficient, with zero distance and unit correlation coefficient indicating the perfect match. The ensemble-mean values of these two measures for individual CMIP5 models are shown in the first two columns of Fig. 5. The models exhibit great variations in their ability to match the observed NINO3.4 spectra, with correlation coefficients ranging from 0.2 to 0.8 and rms distances from 30 to 120; furthermore, there are also substantial sampling variations in the NINO3.4 spectra from multiple runs of a single model (not shown). This means that the ENSO tends to exhibit different epochs of sampling variability in models, in which its strength, spectrum and other properties may vary significantly.30

Fig. 5 Measures of similarity between various characteristics of the observed and CMIP5 simulated NINO3.4 time series. Shown are ensemble averages of these characteristics over multiple runs of individual models (see the model acronyms on the left). The first two columns compare the observed and simulated wavelet spectra, using the root-mean-square (rms) distance and Pearson correlation, both computed in the frequency space. Next three pairs of columns display measures of similarity between the observed and simulated interaction maps analogous to the observed maps of Fig. 1, namely the phase synchronization, phase–phase causality and phase–amplitude causality maps. Each pair presents two distinct similarity measures: the Pearson correlation (corr) in the phase space of the corresponding map, as well as the Adjusted Rand Index (ARI; see text) Full size image

Similarly to the wavelet spectra, the synchronization and causality maps (see examples in Fig. 6 analogous to the observed maps of Fig. 1), vary considerably from model to model (Fig. 5, right columns), as well as between individual runs of a single model (see Supplementary Excel Table). We compared the observed and simulated maps using, once again, the standard Pearson correlation coefficient, as well as the so-called Adjusted Rand Index (ARI), which is especially well suited to measure the similarity of clustered data.31 Both measures were computed for the pairs of interaction maps (observation vs. individual simulation) filled with ones or zeros depending on whether the significant interaction between the processes of different time scales was identified or not (so the colored areas of maps in Figs. 1 and 6 would be filled with ones, and white areas—with zeros); hereafter, we will call the maps so constructed the thresholded binary maps.

Fig. 6 The same as in Fig. 1, but for the individual simulations of CMIP5 models that best match the observed structures: phase synchronization in CanESM model (a); phase–phase causality in CCSM4 model (b); and phase–amplitude causality in CSIRO-mk360 model (c) Full size image

The above similarity measures averaged over the ensembles of individual model simulations (Fig. 5, right columns) are in fact fairly low. For phase synchronization, the highest similarity was detected in the CanESM2 model, at the 0.24 level. The phase synchronization map for the best run of the CanESM2 model (Fig. 6a) indicates synchronization between the processes with the same time scales as in the observed data (Fig. 1a), that is, between the LF, QB, and AC/CT modes. This, however, is more of an exception than a rule, as the time scales of significant phase synchronization in most of the runs do not match the observed time scales.

The similarity levels between phase–phase causality maps from observations and model simulations are about the same as for phase synchronization maps, with maximum ensemble-mean correlation of 0.23 obtained for the CCSM4 model. The causality map for the best run of this model (Fig. 6b) is correlated with the observed map (Fig. 1b) at the 0.41 level, and captures correctly the observed LF–QA, LF–QB, and QA–QB connections. Note that the best matches to the observed maps in terms of phase synchronization and phase–phase causality come from individual simulations of different models, meaning that neither model run was able to capture the entirety of the observed interactions. Finally, no model was able to capture the observed phase–amplitude causal connection between the LF and QB modes (Fig. 1c). The highest ensemble-mean correlation between the observed and simulated maps is only 0.03 (MRI-CGCM3), with the highest correlation of 0.1 for one of the CSIRO-mk360 simulations (the corresponding causality map is shown in Fig. 6c). The phase-synchronization and phase–phase causality maps for the latter run are, however, inferior to those from other models in terms of their similarity to the observed maps (not shown).

To summarize, the CMIP5 models exhibit great sampling variations in the simulated ENSO characteristics. Some of the simulations do exhibit certain aspects of interactions between the processes of different time scales, which match the observed interactions. However, no single simulation is able to reproduce the entire sequence of causal connections inferred from the observed data.

Interactions in conceptual and statistical models

The same interactions were also sought in a conceptual, low-dimensional dynamical model (parametric recharge oscillator: PRO25), as well as in the empirical stochastic (EMR) model of Pacific sea-surface temperatures due to Kondrashov et al.,32 which were both used to produce synthetic NINO3.4 time series. The results for both models are presented in the Supplementary Material. In a nutshell, the PRO dynamical model fails to reproduce any of the observed causal interactions (Supplementary Figs. S4 and S5), whereas the EMR statistical model does offer promising results, especially with regards to reproducing the observed LF → QA → QB phase causality (Supplementary Fig. S6). Some of the initial analyses suggest that the QB modes in the EMR model are inherently present in the algebraic structure of this model’s deterministic propagator, but are strongly influenced by the feedbacks that involve state-dependent (multiplicative) noise (Supplementary Fig. S7).