the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A skewed perspective of the Indian rainfall–El Niño–Southern Oscillation (ENSO) relationship
Justin Schulte
Frederick Policielli
Benjamin Zaitchik
Wavelet coherence is a method that is commonly used in hydrology to extract scaledependent, nonstationary relationships between time series. However, we show that the method cannot always determine why the timedomain correlation between two time series changes in time. We show that, even for stationary coherence, the timedomain correlation between two time series weakens if at least one of the time series has changing skewness. To overcome this drawback, a nonlinear coherence method is proposed to quantify the crosscorrelation between nonlinear modes embedded in the time series. It is shown that nonlinear coherence and autobicoherence spectra can provide additional insight into changing timedomain correlations. The new method is applied to the El Niño–Southern Oscillation (ENSO) and allIndia rainfall (AIR), which is intricately linked to hydrological processes across the Indian subcontinent. The nonlinear coherence analysis showed that the skewness of AIR is weakly correlated with that of two ENSO time series after the 1970s, indicating that increases in ENSO skewness after the 1970s at least partially contributed to the weakening ENSO–AIR relationship in recent decades. The implication of this result is that the intensity of skewed El Niño events is likely to overestimate India's drought severity, which was the case in the 1997 monsoon season, a time point when the nonlinear wavelet coherence between AIR and ENSO reached its lowest value in the 1871–2016 period. We determined that the association between the weakening ENSO–AIR relationship and ENSO nonlinearity could reflect the contribution of different nonlinear ENSO modes to ENSO diversity.
 Article
(10985 KB)  Fulltext XML

Supplement
(323 KB)  BibTeX
 EndNote
The South Asian monsoon, which is the dominant precipitation source for the Indian subcontinent, has been a target for seasonal prediction for well over a century (Blanford, 1884). Despite this long heritage of research, skillful prediction remains a challenge, driving extensive and ongoing research on statistically and dynamically based prediction methods (e.g., REFS). It is difficult to overstate the importance of the South Asian monsoon to the wellbeing of citizens in India. Strong monsoon years have caused catastrophic flooding (Kale, 2012; Sanyal and Lu, 2005) and large landslides (Dortch et al., 2009), while weak monsoons have led to water shortages (Mishra et al., 2016) and crop losses (Parthasarathy et al., 1988; Prasanna, 2014) that resulted in significant food shortages in the past (Fagan, 2009). Thus, while the majority of monsoon forecast studies target the prediction of rainfall totals, the hydrological and agricultural impacts of monsoon variability provide the most pressing motivation for our work.
Much of the research on South Asian monsoon prediction has focused on the relationship between the El Niño–Southern Oscillation (ENSO; Walker and Bliss, 1932) and monsoon strength. During El Niño years, droughts are favored, while rainfall surpluses are favored during La Niña years (Shukla and Paolino, 1983; Kripalani and Kulkarni, 1997). However, there is no onetoone relationship between ENSO and Indian rainfall. As a result, summer rainfall predictions based on ENSO have proven challenging. For example, the 1997–1998 El Niño event was extremely strong, yet climatological Indian monsoon conditions were observed (Shen and Kimoto, 1999; Slingo and Annamalai, 2000). It is therefore important to understand why certain El Niño events are not accompanied by monsoon failures.
There are a few reasons for the challenges faced when predicting Indian rainfall using ENSO. The first reason is that the relationship between ENSO and India's rainfall is nonstationary. As shown by Torrence and Webster (1999), the relationship between ENSO and India's rainfall cycles between periods of high and low coherence. Kumar et al. (1999) found that the relationship between Indian rainfall and ENSO weakened in the 1970s and hypothesized that a southward shift in Walker circulation anomalies associated with ENSO events and increased Eurasian spring and winter surface temperatures were responsible for the weakening relationship. Other works suggest that the changing ENSO–Indian rainfall relationship is the result of the modulating influence of tropical Atlantic sea surface temperatures (SSTs) and the Atlantic Multidecadal Oscillation (AMO; Lu et al., 2006; Kucharski et al., 2007, 2009; Chen et al., 2010). In contrast, Kumar et al. (2006) and Fan et al. (2017) argued that the occurrence of different ENSO flavors (Johnson, 2013), such as the eastern Pacific and central Pacific types, could explain the ENSO–Indian rainfall relationship changes. Other investigators adopted another perspective to explain changes in the ENSO–Indian rainfall relationship and concluded that temporal undulations in the ENSO–Indian rainfall relationship are related to statistical undersampling and stochastic fluctuations (Gershunov et al., 2001; van Oldenborgh and Burgers, 2005; Delsole and Shukla, Cash et al., 2017). In a recent analysis, Yun and Timmermann (2018) showed that changes in the ENSO–Indian rainfall relationship are consistent with a stochastically perturbed ENSO signal and argued that changes in the ENSO–Indian monsoon relationship may not be related to external climateforcing mechanisms.
The second reason for the ENSOrelated prediction challenges is that ENSO itself is a nonstationary phenomenon. Using wavelet analysis, Kestin et al. (1998) found that the interannual variability of ENSO from 1930 to 1960 was dominated by a 4 to 7year periodicity, whereas from 1960 to 1990 the interannual variability was also dominated by a 2 to 5year periodicity. A wavelet power spectral analysis conducted by Torrence and Webster (1999) and Schulte (2016a) showed that ENSO signal energy in the 2 to 7year band undulates throughout the historical period and increases after the 1960s (Schulte, 2016a). These changes in spectral characteristics are relevant to Indian monsoon prediction because differing spectral characteristics among predictors (e.g., ENSO) and predictands (e.g., Indian rainfall) can negatively impact the predictive skill of statistical models (Jiang et al., 2020). Using a waveletbased variance transformation method, Jiang et al. (2020) demonstrated that accounting for differences in spectral characteristics can improve prognostic skill. That study suggests that Indian monsoon prediction could be improved by using waveletbased methods instead of timedomain correlation and regression methods.
The nonlinear characteristics (e.g., skewness) of ENSO are also nonstationary and undergo interdecadal changes (Wu and Hsieh, 2003). Numerous studies have reported an ENSO regime shift in the 1970s in which ENSO began to evolve more nonlinearly than in previous decades (An, 2004, 2009; An and Jin, 2004). It is a curious fact that the ENSO regime shift in the 1970s coincided with the weakening ENSO–Indian rainfall relationship, as documented by Kumar et al. (1999). This observation creates a question that begs to be answered, namely whether nonlinear ENSO regime changes are related to changes in the ENSO–Indian rainfall relationship.
Various mechanisms have been proposed to explain the cause of ENSO skewness. Kang and Kug (2002) suggested that the asymmetry between the magnitude of El Niño and La Niña events is related to the relative westward displacement of zonal wind stress anomalies during La Niña events compared to El Niño events. Jin et al. (2003) and An and Jin (2004) found that ENSO asymmetry is related to nonlinear dynamical heating (NDH), where the magnitude of NDH is related to the propagation characteristics of ENSO. As shown by An and Jin (2004), NDH during strong El Niño events like the 1982–1983 and 1997–1998 events tends to be stronger than during weak El Niño events because SST anomalies tend to propagate eastward. Since the late 1970s there has been a propensity for eastward propagation characteristics of ENSO (Santoso et al., 2013), contrasting with the time period before the 1970s that consisted of the relatively weak El Niño events of 1957–1958 and 1972–1973 (An and Jin, 2004; An, 2009). More recently, Su et al. (2010) showed that vertical temperature advection may have an opposing effect on ENSO asymmetry, and that the asymmetry in the extreme eastern equatorial Pacific is related to meridional ocean temperature advection.
Previous investigators have used different metrics to quantify ENSO asymmetry. To measure the nonlinear character of ENSO, An and Jin (2004) used timedomain metrics, such as skewness and maximum potential intensity (MPI), to quantify the skewness of SST anomalies and the skewness of individual ENSO events, respectively. An (2004) applied a principal component analysis (PCA) to a 21year moving window of tropical Pacific SST skewness and found that the first PCA mode is characterized by positive skewness across the eastern equatorial Pacific and negative skewness across the central equatorial Pacific. This pattern means that interdecadal changes in the nonlinearity of ENSO are associated with positively skewed SST anomalies across the eastern equatorial Pacific, implying that El Niño events are stronger than La Niña events. While the methods implemented in the aforementioned studies provided important insights, they cannot reveal the frequency modes of ENSO that are contributing to the skewness.
Recognizing the limitations of timedomain approaches, Timmermann (2003) conducted a bispectral analysis of the Niño 3 anomaly time series, where a peak (f_{1}; f_{2}) in the bispectrum means there is statistical phase dependence among oscillators with frequencies of f_{1} and f_{2} and f_{1}+f_{2}. That bispectral analysis revealed statistically significant bispectral power at several frequency pairs, including (0.038, 0.038), (0.028, 0.028), (0.0225, 0.0225), (0.0045, 0.032), and (0.0045, 0.045) [months^{−1}]. The peaks (0.0045, 0.032) and (0.0045, 0.045) [months^{−1}] were identified with the nonlinear interactions among the 18 and 2year variability. Although the analysis provided new insights, the Fourierbased analysis could not reveal how the nonlinear nature of ENSO changed with time, which is an important property to capture given how the nonlinear characteristics of ENSO are nonstationary (Santoso et al., 2013). Much like the crosswavelet power (Maraun and Kurths, 2004) and timedomain covariance, bispectral power is not a bounded quantity, and so high bispectral power does not always mean strong phase dependence.
In this study, the deficiencies associated with the abovementioned techniques are addressed using higherorder wavelet analysis, which allows for the quantification of frequencydependent and nonstationary nonlinearities in time series (Van Milligen et al., 1995; Elsayed, 2006; Schulte, 2016b). More specifically, the objectives of the paper are the following: (1) to quantify the nonlinearity of ENSO using higherorder wavelet analysis together with recently developed statistical tests, (2) to determine if different nonlinear modes of ENSO are associated with distinct SST patterns, and (3) to develop nonlinear wavelet coherence methods to test the hypothesis that the breakdown of the ENSO–Indian rainfall relationship in recent decades is related to the shift of ENSO from a linear regime to a nonlinear one. The paper is organized as follows: in Sect. 2, the data used are described. Section 3 includes the description of the implemented methodologies. Results are presented in Sect. 4, and concluding remarks are provided in Sect. 5.
The variability in India rainfall from 1871 to 2016 was analyzed using the allIndia rainfall (AIR; Parthasarathy et al., 1994) time series, which was created by averaging representative rain gauges at various locations across India. The full monsoon season (June–September) and the late monsoon (August–September) season were used to identify possible withinseason variations in the ENSO–allIndia rainfall (ENSO–AIR) relationships. To remove the influence of the annual cycle, AIR time series were converted into anomaly time series by subtracting the 1871–2016 longterm mean for each month from the individual monthly values. The AIR anomaly (hereafter AIR) time series were subsequently standardized by dividing them by their 1871–2016 standard deviations. Because wavelet analysis focuses on specific frequency components that are not impacted by longterm timedomain trends, no detrending of the data was performed.
The monthly data for the Niño 1 and 2, Niño 3, Niño 3.4, and Niño 4 indices (available at: https://www.psl.noaa.gov/gcos_wgsp/Timeseries/, last access: 15 March 2019) from 1871 to 2016 were used to understand how the nonlinear characteristics of SSTs varied from one ENSO region to another. The Niño 1+2 index is the average SST in the region, with latitudinal boundaries of 0^{∘} and 10^{∘} S and longitudinal boundaries of 90 and 80^{∘} W. The Niño 3 index is the average SST in the region, with latitudinal boundaries of 5^{∘} N and 5^{∘} S and longitudinal boundaries of 150 and 90^{∘} W. Variations in SSTs further west were described using the Niño 3.4 and Niño 4 indices, where the Niño 3.4 index is defined as the average SST in the region bounded by 5^{∘} N and 5^{∘} S and 170 and 120^{∘} W. The Niño 4 index is defined as average SSTs in the region bounded by 5^{∘} N and 5^{∘} S and 160^{∘} E and 150^{∘} W. The seasonal cycle was removed from these time series in the same way as it was removed from the AIR time series.
The monthly SST data from 1871 to 2016 were based on the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST1; Rayner et al., 2003). The data at each grid point were converted to monthly anomalies in the same way as they were computed for the ENSO and AIR time series.
3.1 Wavelet analysis
To better diagnose the changes in time series statistics associated with AIR and ENSO, we adopted a wavelet analysis. For a time series, X, comprising data points ${x}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{x}_{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}{x}_{N}$, the continuous wavelet transform is given by the following:
where s is the wavelet scale, ψ_{0} is an analyzing wavelet, δt is a time step (1 month in this study), and n is time. The sample wavelet power spectrum ${\left{W}_{n}\left(s\right)\right}^{\mathrm{2}}$ measured the energy content of a signal at time n and scale s. The commonly used Morlet wavelet, with angular frequency ω=6, was used throughout this paper because it balances time and frequency localization, and because it is commonly used in hydrological and climate studies (Schaefli et al., 2007; Zhang et al., 2007; Holman et al., 2011; Carey et al., 2013). The readers are referred to Torrence and Compo (1998) and Grinsted et al. (2004) for more details about wavelet analysis.
Linear wavelet coherence (Table 1) was used to quantify the linear relationship between two time series as a function of frequency and time. The linear wavelet coherence between two time series X and Y is given by the following:
where S is a smoothing operator (Grinsted et al., 2004) and ${W}_{n}^{\mathrm{XY}}\left(s\right)$ is the crosswavelet power spectrum. Two time series are perfectly coherent (${R}_{n}^{\mathrm{2}}\left(s\right)=\mathrm{1}$) at s if ${\mathit{\varphi}}_{n}^{X}\left(s\right){\mathit{\varphi}}_{n}^{Y}\left(s\right)=c$ over a sufficiently long time interval, where c is a constant, ${\mathit{\varphi}}_{n}^{X}\left(s\right)$ is the phase associated with X, and ${\mathit{\varphi}}_{n}^{Y}\left(s\right)$ is the wavelet phase associated with Y.
In the context of the Indian monsoon, strong coherence between rainfall and a climate pattern (e.g., ENSO) at a scale s indicates shared temporal characteristics between a climate pattern and rainfall. Because the theory supports a causal link between ENSO and monsoon variability through changes in the Walker circulation (Ropelewski and Halpert, 1987), strong coherence means that when ENSO is in a warm (cool) phase at the scale s, negative (positive) rainfall anomalies are preferred. Thus, a periodic climate forcing could create periodicities in an otherwise noisy rainfall time series.
3.2 Higherorder wavelet analysis
Although the wavelet power spectrum is useful for quantifying the signal energy at a scale s and time n, it cannot determine if there is a nonlinear relationship among different frequency components. In fact, the power spectrum can only fully describe time series in frequency space in the case of linear systems in which the output is proportional to the input (King, 1996). As ENSO is nonlinear, we adopted higherorder wavelet methods to address the deficiencies in traditional wavelet methods. The type of nonlinearities considered in this study were quadratic nonlinearities in which the scales s_{1}, s_{2}, and s_{3} satisfied the sum rule as follows:
and the wavelet phases satisfied the following:
These types of nonlinearities arise, for example, when a sinusoid is squared, in which case a harmonic is produced. More generally, quadratic nonlinearities induce time series skewness, which was computed in this study using the following:
where SD in the standard deviation of the time series, and $\phantom{\rule{0.125em}{0ex}}\overline{x}$ is the mean of the time series. Positive (negative) skewness meant that the right (left) tail of the time series distribution is longer than the left (right). In other words, positive (negative) skewness meant that there was a tendency for positive (negative) time series events (i.e., anomalies) to be more intense than negative (positive) ones.
In this paper, quadratic nonlinearities giving rise to time series skewness were quantified using local and global waveletbased autobicoherence methods (Schulte, 2016b). Global autobicoherence (Table 1) was computed as follows:
where
is the global bispectrum, and the hat denotes the complex conjugate. Identical to wavelet coherence, autobicoherence is bounded by 0 and 1, with a value of 1 indicating the strongest possible phase coupling among the phases ϕ_{n}(s_{3}), ϕ_{n}(s_{2}), and ϕ_{n}(s_{1}), such that sum rules of Eqs. (3–4) are satisfied. A peak in the autocoherence spectrum at (s_{1}, s_{2}) indicated that there was quadraticphase coupling among oscillatory modes that was contributing to time series skewness. It is important to note that the autobicoherence method cannot detect other types of nonlinearities such as cubic nonlinearities, whose detection would require trispectra (Collis et al., 1998).
To determine if the strength of the quadraticphase coupling was a function of time, the local autobicoherence spectrum (Schulte, 2016b) given by the following:
was computed, where ${B}_{n}^{X}\left({s}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{s}_{\mathrm{1}}\right)$ is the local bispectrum, given as follows:
and ${s}_{\mathrm{3}}={s}_{\mathrm{1}}/\mathrm{2}$. In this special case, the local autobicoherence spectrum revealed the time evolution of the autobicoherence estimates located along the diagonal slices of the global autobicoherence spectra. Local biphase is as follows:
and was used to measure the skewness and asymmetries of waveforms. A biphase of 0^{∘} meant that the relationship among the scale components produced positive skewness with respect to a horizonal axis so that positive deviations from the mean are larger than negative deviations from the mean (King, 1996; Maccarone, 2013; Schulte, 2016b). On the other hand, a biphase of 180^{∘} indicated negative skewness with respect to the mean. Biphases near −90^{∘} or 90^{∘} indicated that a time series rose (fell) more quickly than it fell (rose).
To be consistent with the wavelet power and coherence analyses, results for the higherorder wavelet analysis were casted in terms of the Fourier period rather than wavelet scale. The Fourier period corresponding to s_{i} was denoted by p_{i}, where the Fourier period is obtained by multiplying s_{i} by 1.03 for the Morlet wavelet (Torrence and Compo, 1998). Thus, the local diagonal slice of the autobicoherence spectra were plotted using the Fourier period p_{1} corresponding to s_{1} as the vertical axis and time as the horizonal axis.
3.3 Statistical hypothesis testing
The statistical significance of all wavelet spectra was evaluated using the cumulative areawise test (Schulte, 2016a, 2019a) to account for the simultaneous testing of multiple hypotheses (Maraun and Kurths, 2004; Maraun et al., 2007). This test evaluated the statistical significance of points in the wavelet domain based on the area of contiguous regions of pointwise significance (i.e., patches) to which they belong, so a larger area implies a greater statistical significance. Given that the patch area can change as the pointwise significance changes, the cumulative areawise test was used to evaluate significance based on the patch area averaged across a set of pointwise significance levels (Schulte, 2019a). The test was applied at the 5 % cumulative areawise significance level, using pointwise significance levels ranging from 0.02 to 0.18, because this choice of pointwise significance levels was shown to result in the cumulative areawise test outperforming the pointwise test in terms of true positive detection for moderate to high signaltonoise ratios, even though the cumulative areawise test is more stringent. The test was performed using the Advanced Biwavelet Wavelet R software package (available at: http://justinschulte.com/wavelets/advbiwavelet.html, last access: 20 July 2020). Technical details of the testing procedure can be found in Schulte (2019a) and Appendix A.
To assess the statistical significance of the global autobicoherence estimates, a modified version of the cumulative areawise test was applied. In the modified version of the cumulative areawise test, the normalized area of patches was computed by dividing the patch area by the product $\widehat{{s}_{\mathrm{1}}}\widehat{{s}_{\mathrm{2}}}$, where $\widehat{{s}_{\mathrm{1}}}$ is the mean first coordinate of the patch, and $\widehat{{s}_{\mathrm{2}}}$ is the mean second coordinate. The reason for this modified normalized area is that dividing area by, say, ${\widehat{{s}_{\mathrm{1}}}}^{\mathrm{2}}$ retained the correlation between the normalized area and s_{2}. The test was applied using the same pointwise significance levels that were used to assess the statistical significance of wavelet power and coherence.
3.4 Higherorder coherence
Although wavelet coherence spectra can provide information regarding how the relationship between two climate variables changes at a scale s, it cannot completely explain why the timedomain correlation between the climate variables temporally fluctuates. The reason is that linear wavelet coherence only examines how well the variance of one time series corresponds to the variance of another at a scale s (Table 1) because linear coherence is determined by the wavelet power spectra of the time series. However, for two time series to be perfectly correlated in the time domain, higher skewness of one time series must also correspond to higher skewness of the other time series.
Recognizing the skewness associated with two time series was correlated in frequency space using the following:
which was referred to as thirdorder coherence (nonlinear coherence, hereafter; the readers are referred to Appendix B for a more general definition of higherorder coherence). In Eq. (11), s_{smooth} is one of the three scales, and ${B}_{n}^{\mathrm{XY}}({s}_{\mathrm{1}},{s}_{\mathrm{2}})$ is the thirdorder crosswavelet power spectrum, which is the product of the bispectrum of X and the conjugate of the bispectrum of Y, the higherorder analog of the crosswavelet power spectrum. The word crossbispectrum was not used to avoid confusion with crossbicoherence analysis (Van Milligen,1995). Like wavelet coherence, the nonlinear coherence is bounded by 0 and 1, with a value of 1 indicating that the bispectra of X and Y at (s_{1}, s_{2}) are perfectly and linearly correlated. The statistical significance of nonlinear coherence was assessed using the cumulative areawise test in the same way as it was used to assess the statistical significance of linear wavelet coherence.
Higherorder wavelet analysis can also be interpreted in terms of linear and nonlinear modes. A linear mode ${\mathit{\gamma}}_{{s}_{i}}^{X}$ is the signal component of X at the scale s_{i} obtained by setting all wavelet coefficients to zero, except those at s_{i}, and taking the inverse wavelet transform of the result. Because linear modes are only composed of a single frequency component, the local crosscorrelation (coherence) between ${\mathit{\gamma}}_{{s}_{i}}^{X}$ and ${\mathit{\gamma}}_{{s}_{i}}^{Y}$ is only impacted by the variances in X and Y at s_{i}. On the other hand, nonlinear coherence measures the local crosscorrelation between the skewness of ${\mathit{\gamma}}_{{s}_{\mathrm{1}}}^{X}+{\mathit{\gamma}}_{{s}_{\mathrm{2}}}^{X}+{\mathit{\gamma}}_{{s}_{\mathrm{3}}}^{X}$ and ${\mathit{\gamma}}_{{s}_{\mathrm{1}}}^{Y}+{\mathit{\gamma}}_{{s}_{\mathrm{2}}}^{Y}+{\mathit{\gamma}}_{{s}_{\mathrm{3}}}^{Y}$ or between ${\mathit{\gamma}}_{{s}_{\mathrm{1}}}^{X}+{\mathit{\gamma}}_{{s}_{\mathrm{1}}/\mathrm{2}}^{X}$ and ${\mathit{\gamma}}_{{s}_{\mathrm{1}}}^{Y}+{\mathit{\gamma}}_{{s}_{\mathrm{1}}/\mathrm{2}}^{Y}$ in the case that s_{1}=s_{2}.
To better understand nonlinear coherence, we supposed that, in the following:
was applied for constants c_{1}, c_{2}, and c_{3}. Adding Eqs. (12) and (13) and subtracting Eq. (14) from the result produced the equality, as follows:
for some constant $\mathrm{K}={\mathrm{c}}_{\mathrm{1}}+{\mathrm{c}}_{\mathrm{2}}{\mathrm{c}}_{\mathrm{3}}.$ Thus, if X was found to be perfectly nonlinearly coherent with Y,then Xand Ywere perfectly coherent at the three scales participating in the quadraticphase coupling. Even if the coherence was perfect at two scales, the relative biphase ${\mathit{\psi}}_{n}^{\mathrm{bi}}\left({s}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{s}_{\mathrm{2}}\right)$ fluctuated randomly if the relative phase difference at the remaining scale fluctuated randomly, so the nonlinear coherence was low. Thus, if nonlinear coherence was high, then there was some nonrandom relationship between X and Y at all three scales, even if high linear coherence was not identified at one or more scales. This theoretical idea indicates that nonlinear coherence can uncover relationships that linear coherence cannot (see Fig. S1 in the Supplement).
The relative biphase difference ${\mathit{\psi}}_{n}^{\mathrm{bi}}\left({s}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{s}_{\mathrm{2}}\right)$ is the higherorder analog of the relative phase difference between two time series. It measures how much the cycle geometry of one time series lags that of another. A lagged biphase of 180^{∘} means that the skewness or asymmetry of the forcing time series is opposite to that of the response. For example, if the forcing has positive skewness, then the response will have negative skewness. If the relative biphase is 0^{∘}, then a negative (positive) skewness of the forcing produces a negative (positive) skewness of the response, contributing to the positive timedomain correlation between the time series. Scales and time points for which nonlinear coherence is high are where the relative biphase is stable.
In this paper, we focused on nonlinear coherence computed along the diagonal slices (p_{1}=p_{2}) of the time series bispectra so that ${s}_{\mathrm{smooth}}={s}_{\mathrm{1}}={s}_{\mathrm{2}}$ in this study. The nonlinear coherence spectra were then plotted using p_{1} as the vertical axis and time as the horizonal axis. High nonlinear coherence at p_{1} and n meant that the skewness or asymmetry between ${\mathit{\gamma}}_{{p}_{\mathrm{1}}}^{X}+{\mathit{\gamma}}_{{p}_{\mathrm{1}}/\mathrm{2}}^{X}$ and ${\mathit{\gamma}}_{{p}_{\mathrm{1}}}^{Y}+{\mathit{\gamma}}_{{p}_{\mathrm{1}}/\mathrm{2}}^{Y}$ were locally crosscorrelated.
To demonstrate the concept of nonlinear coherence, we considered a simple example in which the nonlinear climate forcing time series was given by the following:
and the response to the forcing was given as follows:
In Eqs. (16) and (17), γ(t) is a timevarying nonlinear coefficient, W_{F}(t) is Gaussian white noise associated with the forcing, W_{R}(t) is Gaussian white noise associated with the response, φ= 0 is phase, and ${p}_{\mathrm{1}}=\mathrm{2}{p}_{\mathrm{3}}=\mathrm{32}$. The nonlinear coefficient was assumed to be a linear function of time, as follows:
The effect of the coefficient was to linearly increase the variance of F(t) at p_{3}=16 and increase the strength of the quadraticphase coupling between the modes with periods ${p}_{\mathrm{3}}={p}_{\mathrm{1}}/\mathrm{2}=\mathrm{16}$ and p_{1}=32.
As shown in Fig. 1a, F(t) (black curve) and R(t) (thick green curve) evolve coherently from t=0 to t=200. After t=200, F(t) begins to noticeably exceed R(t) at certain time points (e.g., t=430), while the relationship between them at other points is reversed (e.g., t=450) in the sense that a positive forcing produces a negative response. As a result, the correlation between F(t) and R(t) weakens (Fig. 1b). An inspection of the wavelet coherence spectrum (Fig. 2a) reveals that the coherence at p_{1}=32 is strong and stable, so changes in the relationship strength at that timescale are not the cause of the weakening timedomain correlation. The coherence at all other periods is also stationary by construction, so it is not the changing relationship strength at any timescale that is causing the timedomain correlation weakening. However, the variance of F(t) at p_{3}=16 increases with time (not shown), and the coherence between F(t) and R(t) is also weak at that timescale, implying that larger fluctuations in F(t) at p_{3}=16 are not accompanied by larger fluctuations in R(t). Thus, a variance increase of F(t) is one reason for the weakening timedomain correlation, though the linear coherence and wavelet power methods cannot explain why the skewness of F(t) increases without a corresponding increase in the skewness of R(t) (Fig. 1c).
An inspection of the local autobicoherence spectrum of F(t) (Fig. 2b) reveals that the autobicoherence at p_{1}=32 is increasing with time, indicating that the phase coupling between modes with periods p_{3}=16 and p_{1}=32 is strengthening with time. The biphase of 0^{∘}, as indicated by arrows pointing to the right, confirms that the quadraticphase coupling is contributing to the positive skewness, seen in Fig. 1a, to an increasing degree. Furthermore, the nonlinear coherence between R(t) and F(t) is weak and mostly statistically insignificant at p_{3}=32 (Fig. 2c), implying the skewness of ${\mathit{\gamma}}_{\mathrm{16}}^{\mathrm{F}}+{\mathit{\gamma}}_{\mathrm{32}}^{\mathrm{F}}$ is uncorrelated with the skewness of ${\mathit{\gamma}}_{\mathrm{16}}^{\mathrm{R}}+{\mathit{\gamma}}_{\mathrm{32}}^{\mathrm{R}}$, where ${\mathit{\gamma}}_{\mathrm{16}}^{\mathrm{F}}+{\mathit{\gamma}}_{\mathrm{32}}^{\mathrm{F}}$ is the sum of the cosines in Eq. (16) and the components of W_{F}(t) at p_{3}=16 and p_{1}=32. The nonlinear mode ${\mathit{\gamma}}_{\mathrm{16}}^{\mathrm{R}}+{\mathit{\gamma}}_{\mathrm{32}}^{\mathrm{R}}$ is the sum of the cosine in Eq. (17) and the components of W_{R}(t) at p_{3}=16 and p_{1}=32. Thus, the skewness of R(t) in the time domain is practically uncorrelated with the skewness of F(t) because the skewness of F(t) is solely related to the phase coupling between the modes with periods p_{3}=16 and p_{1}=32. Thus, the increase in the skewness of F(t) also contributes to the weakening timedomain correlation.
The lack of nonlinear coherence at timescales for which F(t) is nonlinear has implications for empirical prediction. At time points where F(t) is positively skewed, R(t) is overestimated because R(t) is not inheriting the skewness of F(t). That is, if one created a linear regression model based on the relationship between F(t) and R(t) from t=0 to t=200, one would find that a forcing value of, say, 1 would produce a response close to 1. If the same model was used to predict R(t) at, say, t=430, one would predict that the forcing with a value of around 2 should result in a response near 2. However, because the relatively large value F(430) results from skewness, and R(t) is uncorrelated with its skewness, the response is only as strong as the part of F(t) not resulting from the quadraticphase coupling. The more nonlinear F(t) becomes, the more F(t) will overestimate R(t) when F(t) is positively skewed. Similarly, the positive forcing produces a negative response at t=450 because of skewness and not simply because of a change in variance. Nonlinear coherence allows for the quantification and identification of these timedomain aberrations.
The weakening relationship shown in Fig. 1b could lead a researcher studying a hydrological process to believe that another directforcing mechanism must be influencing the hydrological process. This belief could lead to the application of partial wavelet coherence (Ng and Chan, 2012) and partial correlation analyses to identify another influential forcing mechanism. However, in this case, there are no other directforcing mechanisms; the weakening timedomain relationship is solely related to how F(t) transitioned from a linear regime to a nonlinear regime. This theoretical result suggests that hydrological studies using wavelet coherence should also consider the nonlinearity of the times series.
4.1 The ENSO and Indian monsoon time series and their timedomain relationship
The contrasting Niño 1+2 and Niño 4 indices are shown in Fig. 3. For the Niño 1+2 time series (Fig. 3a), a few recent and notably intense warm events are located around 1982–1983, 1997–1998, and 2015–2016, coinciding with the strongest El Niño events in recent decades (McPhaden, 1999; Hu and Fedorov, 2017; Santoso et al, 2017). A few notably intense Niño 1+2 events are also seen in the 1800s, indicating that intense ENSO events are not unique to recent decades. An inspection of Fig. 3a also reveals that the recent intense warm Niño 1+2 events are also skewed in the sense that they are stronger than the surrounding cool Niño 1+2 events. Unlike the 1982–1983 and 1997–1998 Niño 1+2 events, the 1982–1983 and 1997–1998 warm events for the Niño 4 time series are unremarkable (Fig. 3b). Furthermore, cool Niño 4 events are preferentially stronger than warm events after the 1960s, suggesting an intensification of the negative skewness.
The 20year sliding skewness time series of the Niño 1+2 index (Fig. 4) reveals enhanced skewness during the 1800s, nearzero skewness around the 1930s and early 1940s, and especially enhanced skewness after the 1970s associated with an upward trend in skewness beginning around the 1940s. In contrast to the Niño 1+2 index, the skewness of the Niño 4 index becomes more negative after the 1960s, and the magnitude of the skewness is generally smaller than that of the Niño 1+2 time series. This finding suggests that the transition of the Niño 1+2 time series to a nonlinear regime is more pronounced than the transition associated with the Niño 4 time series.
Interestingly, a 20year sliding skewness analysis of AIR (Fig. 4) reveals that the skewness of June–September AIR remains close to zero until the 1990s, despite the upward trend in Niño 1+2 skewness beginning in the 1940s (Fig. 4a). Similarly, the skewness of August–September AIR does not increase to the extent that Niño 1+2 skewness does (Fig. 4b). On the other hand, the skewness of August–September AIR could be negatively correlated with the skewness of the Niño 4 index after the 1960s, consistent with how the August–September AIR and the Niño 4 index are negatively correlated. The skewness of June–September AIR becomes more negative in the 1990s and 2000s, but it is unclear if that negative skewness is related to ENSO, noise, or another climate pattern because the skewness of the Niño 1+2 and Niño 4 indices do not change as abruptly. Negative June–September AIR skewness is accompanied by enhanced positive skewness of the Niño 1+2 index prior to the 1940s, which is consistent with how June–September AIR is negatively correlated with the Niño 1+2 index during that time period (Fig. 5a).
The differences in skewness shown in Fig. 4 suggests that the correlation between the ENSO time series and AIR degrades after the 1970s, which is confirmed by the 20year sliding correlation between the June–September AIR and ENSO time series (Fig. 5a). The relationship with the Niño 1+2 generally weakens from the 1800s to the 2000s. In contrast, the June–September AIR–Niño 4 index relationship appears to have no longterm trend, resulting in the Niño 4 index becoming more strongly correlated with AIR than the Niño 1+2 index after the 1970s.
The stronger AIR–Niño 4 index relationship compared to the AIR–Niño 1+2 index relationship after the 1970s is more evident in the August–September analysis (Fig. 5b). An abrupt weakening of the August–September AIR–Niño 1+2 index relationship occurs around the 1970s, with the relationship reversing around the 1990s. A comparison of Figs. 4b and 5b reveals that the weakening and reversal of the relationship occurs during the time period when the Niño 1+2 index is especially skewed. The difference in the magnitudes of Niño 1+2 and Niño 4 skewness after the 1970s could explain why the August–September AIR–Niño 1+2 index relationship weakens more abruptly than the AIR–Niño 4 index relationship. Thus, a further investigation is needed to better understand the temporal changes in ENSO statistics and their impact on the ENSO–AIR relationship.
4.2 Wavelet power analysis and coherence
The wavelet power spectra associated with the Niño 1+2 and Niño 4 time series (Fig. 6) reveal enhanced variance in the 16 to 64 month band after 1965 for all the time series. The appearance of holes in the contoured regions suggests that there are oscillatory modes with nearby frequencies (Schulte, et al., 2015), though the wavelet power spectra cannot determine if there is quadraticphase coupling between the oscillatory modes.
The linear wavelet coherence spectrum, shown in Fig. 7, indicates that the AIR relationship with the Niño 1+2 and Niño 4 indices in the 16 to 64 month band breaks down after 1995, which is consistent with the findings from the sliding correlation analysis shown in Fig. 5. The relationship between AIR and these ENSO indices also weakens around 1925, but this weakening does not appear in the sliding correlation analysis. Note that the lack of linear coherence after 1995 coincides with the enhanced ENSO variance (Fig. 6), implying that higher ENSO variance need not be associated with higher AIR variance at those timescales, so changes in ENSO variance could be contributing to the weakening ENSO–AIR timedomain correlation. However, ENSO skewness is also enhanced during this time period (Fig. 4), so the weakening relationships may not simply be related to ENSO variance. Thus, further analysis is needed to extract information unrevealed by the linear wavelet power and coherence methods.
4.3 Local autobicoherence of ENSO
Figure 8 shows that the local autobicoherence spectra of all ENSO time series contain statistically significant local autobicoherence, but the spectrum of the Niño 4 index is only associated with a few statistically significant regions, such as the one around 2015 at a period of 32 months. For the Niño 1+2 index, there is an intensification in the autobicoherence spectrum after the 1970s, which is consistent with the ENSO regime shift (Santoso et al., 2013). A comparison of Figs. 4 and 8 reveals that enhanced skewness coincides with stronger autobicoherence on the 32 to 64 month timescale, suggesting that the skewness partially arises from the stronger quadraticphase coupling among oscillatory modes with periods ranging from 32 to 64 months.
To confirm that the nonlinearphase coupling identified in Fig. 8 is associated with skewed waveforms, we inspected the corresponding local biphase spectra (not shown). It was found that the biphase in the 42 to 64 month band is generally 0^{∘}, so the nonlinearphase coupling in that band contributes to the positive skewness of the 1982–1983, 1997–1998, and 2015–2016 events.
4.4 Nonlinear coherence between allIndia rainfall and ENSO
The results shown in Fig. 9 indicate that the nonlinear wavelet coherence between AIR and the time series for the Niño 1+2 and Niño 4 indices is statistically significant in the 32 to 64 month band, mainly prior to the 1980s. The nonlinear coherence in this band appears to peak around the 1972–1973 El Niño event, indicating that an increase in the positive skewness of ENSO should tend to coincide with the enhanced negative skewness of AIR around this time. However, much of the statistically nonlinear coherence is located during the time period when ENSO is more linear than it has been in recent decades (Fig. 8), so the effects of the nonlinearities are small regardless of the nonlinear wavelet coherence. In contrast, the autobicoherence of the Niño 1+2 time series in the 32 to 64 month band is statistically significant and high after the 1970s (Fig. 8), so the lack of nonlinear coherence after the 1980s, shown in Fig. 9, is expected to impact the timedomain correlation more strongly, much like the theoretical situation shown in Figs. 1 and 2. Our results are consistent with this theoretical idea because the AIR–Niño 1+2 relationship weakens more than the AIR–Niño 4 relationship after the 1970s (Fig. 5), which is expected because the Niño 1+2 index is more nonlinear than the Niño 4 index during this time period. However, unlike the theoretical example shown in Fig. 2, the linear coherence between the ENSO time series and AIR also weakens around the 1990s (Fig. 7), so the weakening relationship could be the result of a combination of factors that includes ENSO nonlinearity.
The 20year sliding mean of the ENSO autobicoherence, coherence, and nonlinear coherence averaged in the 32 to 64 month band further highlights the impact of ENSO nonlinearity. As shown in Fig. 10a, the sliding mean nonlinear coherence between the Niño 1+2 index and AIR fluctuates less than the linear coherence and reaches a clear global maximum around the 1970s before rapidly declining to a global minimum around the late 1990s when the Niño 1+2 index is very nonlinear. As shown in Fig. 10a, the Niño 1+2 autobicoherence peaks around the same time that the August–September AIR–Niño 1+2 index correlation is positive. In fact, the correlation between the sliding September–August AIR–Niño 1+2 correlation time series and the sliding Niño 1+2 autobicoherence time series is 0.81, much higher than the correlation with the linear coherence ($r=\mathrm{0.11}$) and nonlinear coherence ($r=\mathrm{0.34}$) time series. These results support the idea that the Niño 1+2 regime shift impacted the weakening timedomain correlation. On the other hand, the correlation between the Niño 4 autobicoherence and the August–September AIR–Niño 4 correlation time series is weak, so changes in the nonlinearity of the Niño 4 index are unlikely to contribute substantially to changes in the AIR–Niño 4 relationship. Nevertheless, this result agrees with theory that suggests that nonlinearity is only an important contributor when the time series is highly nonlinear, which is not the case for the Niño 4 index because of the low autobicoherence (Figs. 8b and 10b). Because the nonlinear coherence between AIR and indices for the Niño 1+2 and Niño 4 is weak (Figs. 9 and 10), the more pronounced change in the August–September AIR–Niño 1+2 correlation reflects the more intense increase in Niño 1+2 nonlinearity compared to that of the Niño 4 index in recent decades.
4.5 A possible explanation for the ENSO nonlinearity impacts
To better understand the association between ENSO nonlinearity and the AIR–ENSO relationship, the global autobicoherence spectra associated with the ENSO time series were first computed (Fig. 11). Then, the autobicoherence of SSTs associated with a few select peaks (p_{1} and p_{2}) in Fig. 11 were computed at each grid point in the domain bounded by 20^{∘} N and 20^{∘} S and by 146^{∘} E and 80^{∘} W. The peaks were selected based on the autobicoherence spectra of the Niño 3.4 and Niño 1+2 indices. To select the peaks, local maxima in the autobicoherence within the statistically significance regions shown in Fig. 11 were identified, where points associated with local maxima were chosen because they were associated with the clearest patterns.
The spatial structure of global autobicoherence corresponding to the peaks in the Niño 3.4 autobicoherence spectrum are shown in Fig. 12. The autobicoherence associated with the pair (31, 31) is greatest across the central equatorial Pacific, with the overall spatial pattern being reminiscent of a central Pacific El Niño (Lee and McPhaden, 2010). This result suggests that the phase coupling between the 31 month mode and the 15.5 month mode could be related to the occurrence of central Pacific El Niño events (Sect. 5). In contrast, the autobicoherence pattern associated with the pair (56, 56) is more uniform, with autobicoherence slightly greater across the extreme eastern equatorial Pacific than the central equatorial Pacific. This pattern is reminiscent of an eastern Pacific El Niño. Like the pattern corresponding to the pair (31, 31), the autobicoherence for the pair (105, 57) tends to be greater across the central equatorial Pacific. Our findings suggest that different nonlinear modes contribute to different ENSO flavors. Although An and Jin (2004) and Burgers and Stephenson (1999) showed that the skewness is greatest across the eastern equatorial Pacific, we determined that such a timedomain approach is unable to capture frequencydependent patterns in nonlinearity.
The spatial autobicoherence plots associated with the peaks in the Niño 1+2 autobicoherence spectrum are shown in Fig. 13. The autobicoherence associated with the pairs (148, 53) and (148, 105) is strong across the eastern equatorial Pacific but weak across the central equatorial Pacific, suggesting that the quadraticphase coupling between the 148 and 105 month modes and between the 148 and 53 month modes is associated with the skewness of eastern equatorial Pacific SSTs. The pattern associated with the pair (62, 44) is reminiscent of an eastern Pacific El Niño, and the autobicoherence associated with the pair (88, 88) is relatively weak across the entire equatorial Pacific. A comparison of Figs. 12 and 13 shows that there is a tendency for autobicoherence to be greater across the eastern equatorial Pacific than the central equatorial Pacific, which is consistent with how SSTs across the eastern equatorial Pacific are most skewed (Burgers and Stephenson, 1999; An and Jin, 2004).
The nonlinear nature of ENSO was examined using higherorder wavelet methods. The autobicoherence spectra of the Niño 1+2 and Niño 4 indices time series revealed that ENSO skewness arose from the quadraticphase coupling of modes with various periods. For the Niño 1+2 index, the quadraticphase coupling after the 1970s was especially strong, which is consistent with how ENSO underwent a regime shift around the 1970s (Santoso et al., 2013) that was marked by an increase in ENSO skewness. Although the Niño 3.4 time series was not considered in detail in this study, an autobicoherence analysis of the time series (not shown) revealed phase coupling between modes with periods of 31 and 15.5 months in addition to coupling between modes with periods of 61 and 30.5 months. The phase coupling between the 31 and 15.5 modes was found to be especially strong after 1995, whereas the quadraticphase coupling between the 61 and 30.5 month modes was found to intensify after the 1970s. These additional results suggest that nonlinear modes vary in intensity and time of occurrence.
The evolution of SSTs across the Niño 4, Niño 3.4, Niño 3, and Niño 1+2 regions was found to be nonlinear, but the degree to which the time series are nonlinear is different (Fig. 11). Overall, the Niño 1+2 time series were found to be the most nonlinear, while the Niño 4 index was found to be the most linear. The spatial patterns associated with the nonlinearities depend on the frequency components contributing to the nonlinearities. For example, the quadraticphase coupling between the modes with periods of 31 and 15.5 months was found to be the strongest in the central equatorial Pacific and weakest across the eastern equatorial Pacific. This finding suggests that the more frequent occurrence of central Pacific El Niño events in recent decades (Lee and McPhaden, 2010) could be linked to the strengthening of this quadraticphase coupling, which could explain the relationship between ENSO nonlinearity and changes in the ENSO–AIR relationship because central Pacific El Niño events have been shown to be more effective at creating droughtinducing subsidence over India (Kumar et al., 2006).
The results from the present and previous studies (Fan et al., 2017) support the idea that changes in the ENSO–AIR relationship are related to ENSO flavors because ENSO nonlinearity appears to be related to ENSO flavors (Figs. 12 and 13), which is in opposition to the findings of other work showing that the changes are related to sampling variability or to noise. According to Yun and Timmermann (2018), the changes in the timedomain correlation between AIR and ENSO are consistent with the assumption that AIR is the sum of the ENSO signal and Gaussian white noise (i.e., AIR = ENSO and white noise). However, for this hypothesis to hold, the difference in AIR–ENSO must be Gaussian white noise. As shown in this study, the nonlinear wavelet coherence between ENSO metrics and AIR is weak, and ENSO–AIR contains periodicities (Fig. S2), which means that AIR is not simply a stochastically perturbed ENSO signal as noise does not contain periodicities. The retention of nonGaussian noise features was certainly the case for R(t)–F(t) in the example in Sect. 3.5 because the difference retains the cosine function with a period of 16.
The fact that nonlinear coherence between rainfall and ENSO is determined by linear coherence between ENSO and rainfall at two or three frequencies means that the changing timedomain correlation could be more fully understood by determining why linear coherence changes at the frequencies that contribute to ENSO skewness. Such an analysis could provide a more mechanistic perspective than the theoretical perspective adopted in this study. A preliminary analysis showed that there was enhanced linear coherence between the North Atlantic Oscillation index and AIR after 1995 in the 16 to 64 month band associated with ENSO nonlinearity. This result suggests that conditions across the North Atlantic (Kakade and Dugam, 2000; Bhatla et al., 2016) could influence the nonlinear coherence between ENSO and AIR and, thus, the corresponding timedomain correlation.
The tools used and developed in this study may have important applications for understanding how forecasting systems replicate Indian rainfall and its associated teleconnections. These methods, for example, could determine if forecasting systems can reproduce nonlinear characteristics of climate time series. As such, an R software package has been developed to implement these methods (available at: http://justinschulte.com/wavelets/advbiwavelet.html; last access: 22 November 2019). These methods could provide new directions for improving current forecasting systems and, ultimately, predictions of Indian rainfall.
The first step (Step 1) in assessing the cumulativearea significance of a point was the calculation of the N=12 sets as follows:
where each set is the subset of the wavelet domain consisting of points for which the wavelet quantities are pointwise significant at the α_{i} level (i.e., the pointwise test pvalue, ρ_{pw}, is less than α_{i}). In this paper, α_{1}=0.02, α_{12}= 0.18, and ${\mathit{\alpha}}_{i+\mathrm{1}}{\mathit{\alpha}}_{i}=\mathrm{0.02}$. In the second step (Step 2), a geometric pathway about x was computed, where a geometric pathway is a nested sequence, as follows:
such that, in the following:
are pathcomponents of ${P}_{\mathrm{pw}}^{i}$ containing x. The equivalence relation ∼ on ${P}_{\mathrm{pw}}^{i}$ makes two points in ${P}_{\mathrm{pw}}^{i}$ equivalent if they can be connected by a continuous path in ${P}_{\mathrm{pw}}^{i}$. The third step (Step 3) involved the calculation of the normalized area corresponding to ${P}_{i}^{x}$. The normalized area is defined as the patch area divided by the square of the meanscale coordinate of the patch, where ${A}_{i}^{x}$ was assumed to be 0 if ${P}_{i}^{x}=\mathit{\varphi}$ or ${P}_{i}^{x}=\left\{x\right\}$. The critical area ${A}_{i}^{\mathrm{crit}}$ was obtained by computing the (1−α_{c})th percentile of the null distribution of normalized areas corresponding to the significance level α_{i}, where α_{c} is the significance level of the cumulative areawise test. The null distributions were constructed by generating 1000 patches at the α_{i} significance level under the null hypothesis of red noise. The final step (Step 4) was to compute the following:
where ${\mathit{\lambda}}_{i}^{x}=$ 2 if ${A}_{i}^{x}/{A}_{i}^{\mathrm{crit}}\mathit{>}\mathrm{1}$ and ${\mathit{\lambda}}_{i}^{x}=$ 0 if ${A}_{i}^{x}/{A}_{i}^{\mathrm{crit}}\mathit{<}=\mathrm{1}$. The wavelet quantity at the point x was deemed statistically significant at the α_{c} cumulative areawise level if r^{x}>1.
For p>1, the (p+1)thorder polyspectrum of a time series X is given by the following:
where
The thirdorder polyspectrum is the bispectrum, and the fourthorder polyspectrum is the trispectrum (Collis et al., 1998) which identifies the frequency components contributing to kurtosis. The (p+1)thorder coherence between two time series is given as follows:
where ${B}_{n}^{\mathrm{XY}}({s}_{\mathrm{1}}{s}_{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}{s}_{p})$ is the (p+1)thorder cross spectrum given by the following:
When p=2, Eq. (B3) measures the local crosscorrelation between skewness, and when p=3, the equation measures the local crosscorrelation between kurtosis.
Data for the Indian rainfall can be accessed through https://www.tropmet.res.in/DataArchival51Page (Indian Institute of Tropical Meteorology, 2019). The monthly ENSO indices are available at https://www.psl.noaa.gov/gcos_wgsp/Timeseries (NOAA/OAR/ESRL PSD, 2019). An R software package used to implement the new methods can be found on the author's web page (available at: http://justinschulte.com/wavelets/advbiwavelet.html; Schulte, 2019b).
The supplement related to this article is available online at: https://doi.org/10.5194/hess2454732020supplement.
JS conceived the study design and performed the experiments. BZ and FP provided feedback about the experiments and helped with writing the paper.
The authors declare that they have no conflict of interest.
We thankful for the helpful suggestions provided by James DossGollin and the anonymous reviewer.
This research has been supported by the National Aeronautics and Space Administration (grant no. NNX26AN38G.).
This paper was edited by Pierre Gentine and reviewed by James DossGollin and one anonymous referee.
An, S.I.: Interdecadal changes in the El NiñoLa Niña symmetry, Geophys Res. Lett., 31, L23210, https://doi.org/10.1029/2004GL021699, 2004.
An, S.I.: A review of interdecadal changes in the nonlinearity of the El Nino–Southern Oscillation, Theor. Appl. Climatol., 97, 29–40, 2009.
An, S.I. and Jin, F.F: Nonlinearity and asymmetry of ENSO., J. Climate, 17, 2399–2412, 2004.
Ashok, K., Guan, Z., and Yamagata, T.: Impact on the Indian Ocean dipole on the relationship between the Indian monsoon rainfall and ENSO, Geophys Res. Lett., 28, 4499–4502, 2001.
Ashok, K., Guan, Z, Saji, N. H., and Yamagata, T.: Individual and combined influences of ENSO and the Indian Ocean dipole on the Indian summer monsoon, J. Climate, 17, 3141–3155, 2004.
Bhatla, R., Singh, A. K., Mandal, B., Ghosh, S., Pandey, S. N., and Abhijit, S.: Influence of North Atlantic Oscillation on Indian Summer Monsoon Rainfall in Relation to QuasiBinneal Oscillation, Pure Appl. Geophys., 173, 2959–2970, 2016.
Blanford, H. F.: On the connexion of the Himalaya snowfall with dry wind and seasons of drought in India, Proc. R Soc. Lond., 37, 3–22, 1884.
Burgers, G. and Stephenson, D. B: The “Normality” of ENSO, Geophys Res. Lett., 26, 1027–1030, 1999.
Carey, S. K., Tetzlaff, D., Buttle, J., Laudon, H., McDonnell, J,. McGuire, K., Seibert, J., Soulsby, C., and Shanley, J.: Use of color maps and wavelet coherence to discern seasonal and inter annual climate influences on streamflow variability in northern catchments, Water Resour. Res., 49, 6194–6207, 2013.
Cash, B. A., Barimalala, R., Kinter, J. L., Altshuler, E. L., Fennessy, M. J., Manganello, J. V., Molteni, F., Towers, P., and Vitart, F.: Sampling variability and the changing ENSO–monsoon relationship, Clim. Dynam., 48, 4071–4079, 2017.
Chen, W., Dong, B., and Lu, R.: Impact of the Atlantic Ocean on the multidecadal fluctuation of El NiñoSouthern OscillationSouth Asian monsoon relationship in a Coupled General Circulation Model, J. Geophys. Res., 115, D17109, https://doi.org/10.1029/2009JD013596, 2010.
Collis, W. B., White, P. R., and Hammond, J. K.: Higherorder Spectra: The Bispectrum and Trispectrum, Mech. Syst. Signal Pr., 12, 375–394, 1998.
DelSole, T. and Shukla, J.: Climate models produce skillful predictions of Indian summer monsoon rainfall, Geophys. Res. Lett., 39, L09703, https://doi.org/10.1029/2012GL051, 2012.
Dortch, J. M., Owen, L. A., Haneberg, W. C., Caffee, M. W., Dietsch, C., and Kamp, U.: Nature and timing of large landslides in the Himalaya and Transhimalaya of northern India, Quaternary Sci. Rev., 28, 1037–1054, 2009.
DelSole, T. and Shukla, J.: Climate models produce skillful predictions of Indian summer monsoon rainfall, Geophys Res. Lett, 39, L09703, https://doi.org/10.1029/2012GL051279, 2012.
Elsayed, M. A. K.: Wavelet Bicoherence Analysis of Wind–wave Interaction, Ocean Eng., 33, 458–470, 2006.
Fagan, B.: Floods, famines, and emperors: El Niño and the fate of civilizations, Basic Books, New York, 2009.
Fan, F., Dong, X., Fang, X., Xue, F., Zheng, F., and Zhu, J.: Revisiting the relationship between the south Asian summer monsoon drought and El Niño warming pattern, Atmos. Sci. Lett., 18, 175–182, 2017.
Gershunov, A., Schneider, N., and Barnett, T.: Lowfrequency modulation of the ENSOIndian monsoon rainfall relationship: Signal or noise?, J. Climate, 14, 2486–2492, 2001.
Grinsted, A., Moore, J. C., and Jevrejeva, S.: Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlin. Processes Geophys., 11, 561–566, https://doi.org/10.5194/npg115612004, 2004.
Holman, I. P., RivasCasado, M., Bloomfield, J. P., and Gurdak, J. J.: Identifying nonstationary groundwater level response to North Atlantic ocean–atmosphere teleconnection patterns using wavelet coherence, Hydrogeol. J., 19, 1269–1278, https://doi.org/10.1007/s1004001107559, 2011.
Hu, S. and Fedorov, A. V.: The extreme El Nino of 2015–2016 and the end of global warming hiatus, Geophys. Res. Lett., 4415 3816–24, 2017
Indian Institute of Tropical Meteorology: available at: https://www.tropmet.res.in/DataArchival51Page, last access: 10 October 2019.
Jiang, Z., Sharma, A., and Johnson, F.: Refining Predictor Spectral Representation Using Wavelet Theory for Improved Natural System Modeling, Water Resour. Res., 56, e2019WR026962, https://doi.org/10.1029/2019wr026962, 2020.
Jin, F.F., An, SI., Timmermann, A., and Zhao, J.: Strong El Nino events and nonlinear dynamical heating, Geophys. Res. Lett., 30, 1120, https://doi.org/10.1029/2002GL016356, 2003.
Johnson, N. C.: How Many ENSO Flavors Can We Distinguish?, J. Climate, 26, 4816–4827, 2013.
Kakade, S. B. and Dugam, S. S.: The simultaneous effect of NAO and SO on the monsoon activity over India, Geophys. Res. Lett., 27, 3501–3504, 2000.
Kale, V.: On the link between extreme floods and excess monsoon epochs in South Asia, Clim. Dynam., 39, 1107–1122, 2012.
Kang, I.S. and Kug, J.S.: El Nino and La Niña sea surface temperature anomalies: Asymmetry characteristics associated with their wind stress anomalies, J. Geophys. Res., 107, 4372, https://doi.org/10.1029/2001JD000393, 2002.
Kestin, T. A., Karoly, D. J., Yano, J.I., and Rayner, N. A.: Time–frequency variability of ENSO and stochastic simulations, J. Climate, 11, 2258–2272, 1998.
King, T.: Quantifying Nonlinearity and Geometry in Time Series of Climate, Quaternary Sci. Rev., 15, 247–266, 1996.
Kripalani, R. H. and Kulkarni, A.: Climatic impact of El Nino/La Nina on the Indian monsoon: A new perspective, Weather, 52, 39–46, 1997.
Kucharski, F., Bracco, A., Yoo, J. H., and Molteni, F.: Lowfrequency variability of the Indian monsoon–ENSO relationship and the tropical Atlantic: The “weakening” of the 1980s and 1990s, J. Climate, 20, 4255–4266, 2007.
Kucharski, F., Bracco, A., Yoo, J. H., Tompkins, A. M., Feudale, L., Ruti, P., and Dell'Aquila, A.: A gillMatsunotype mechanism explains the tropical Atlantic influence on African and Indian monsoon rainfall, Q. J. Roy. Meteor. Soc., 135, 569–579, 2009.
Kripalani, R. H. and Kulkarni, A.: Climatic impact of El Nino/La Nina on the Indian monsoon: A new perspective, Weather, 52, 39–46, 1997.
Kumar, K. K., Rajagopalan, B., and Cane, M. A.: On the weakening relationship between the Indian monsoon and ENSO, Science, 284, 2156–2159, 1999.
Kumar, K. K., Rajagopalan, B., Hoerling, M., Bates, G., and Cane, M. A.: Unraveling the Mystery of Indian Monsoon Failure During El Niño, Science, 314, 115–119, 2006.
Lee, T. and McPhaden, M. J.: Increasing intensity of El Niño in the centralequatorial Pacific, Geophys. Res. Lett., 37, L14603, https://doi.org/10.1029/2010GL044007, 2010.
Lu, R., Dong, B., and Ding, H.: Impact of the Atlantic Multidecadal Oscillation on the Asian summer monsoon, Geophys. Res. Lett., 33, L24701, https://doi.org/10.1029/2006GL027655, 2006.
McPhaden, M. J.: Genesis and evolution of the 1997–98 El Nino, Science, 283, 950–954, 1999.
Maccarone, T. J.: The Biphase Explained: Understanding the Asymmetries in Coupled Fourier Components of Astronomical Timeseries, Mon. Not. R. Astron. Soc., 435, 3547, https://doi.org/10.1093/mnras/stu1824, 2013.
Maraun, D. and Kurths, J.: Cross wavelet analysis: significance testing and pitfalls, Nonlin. Processes Geophys., 11, 505–514, https://doi.org/10.5194/npg115052004, 2004.
Maraun, D., Kurths, J., and Holschneider, M.: Nonstationary Gaussian Processes in the Wavelet Domain: Synthesis, Estimation, and Significance Testing, Phys. Rev. E, 75, 016707, https://doi.org/10.1103/PhysRevE.75.016707, 2007.
Mishra, V., Aadhar, S., Asoka, A., Pai, S., and Kumar, R.: On the frequency of the 2015 monsoon season drought in the IndoGangetic Plain, Geophys. Res. Lett., 43, 12–102, 2016.
Ng, E. K. W. and Chan, J. C. L.: Geophysical applications of partial wavelet coherence and multiple waveletcoherence, J. Atmos. Ocean. Tech., 29, 1845–1853, 2012.
NOAA/OAR/ESRL PSD: https://www.psl.noaa.gov/gcos_wgsp/Timeseries/, last access: 15 March 2019.
Parthasarathy, B., Munot, A. A., and Kothawale, D. R.: Regression model for estimation of Indian foodgrain production from summer monsoon rainfall, Agr. Forest Meteorol., 42, 167–182, 1988.
Parthasarathy, B., Munot, A. A., and Kothawale, D. R.: AllIndia monthly and seasonal rainfall series: 1871–1993, Theor. Appl. Climatol., 49, 217–224, 1994.
Prasanna, V.: Impact of monsoon rainfall on the total foodgrain yield over India, J. Earth Syst. Sci., 123, 1129–1145, 2014.
Rayner, N. A., Parker, D. E., Horton, E. B., Folland, C. K., Alexander, L. V., Rowell, D. P., Kent, E. C., and Kaplan, A.: Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century, J. Geophys. Res., 108, 4407, https://doi.org/10.1029/2002JD002670,2003.
Ropelewski, C. F and Halpert, M. S.: Global and regional scale precipitation patterns associated with the El Niño/Southern Oscillation, Mon. Weather Rev., 115, 1606–1626, 1987.
Santoso, A., McGregor, S., Jin, F.F., Cai, W., England, M. H., An, S.I., McPhaden, M. J., and Guilyardi, E.: Latetwentiethcentury emergence of the El Niño propagation asymmetry and future projections, Nature, 504, 126–130, 2013.
Santoso, A., McPhaden, M. J., and Cai, W.: The defining characteristics of ENSO extremes and the strong 2015/2016 El Niño, Rev. Geophys, 55, 1079–1129, 2017.
Sanyal, J. and Lu, X. X.: Remote sensing and GISbased flood vulnerability assessment of human settlements: a case study of Gangetic West Bengal, India, Hydrol. Process., 19, 3699–3716, 2005.
Schaefli, B., Maraun, D., and Holschneider, M.: What drives high flow events in the Swiss Alps? Recent developments in wavelet spectral analysis and their application to hydrology, Adv. Water Resour., 30, 2511–2525, 2007.
Schulte, J. A.: Cumulative areawise testing in wavelet analysis and its application to geophysical time series, Nonlin. Processes Geophys., 23, 45–57, https://doi.org/10.5194/npg23452016, 2016a.
Schulte, J. A.: Wavelet analysis for nonstationary, nonlinear time series, Nonlin. Processes Geophys., 23, 257–267, https://doi.org/10.5194/npg232572016, 2016b.
Schulte, J. A.: Statistical hypothesis testing in wavelet analysis: theoretical developments and applications to Indian rainfall, Nonlin. Processes Geophys., 26, 91–108, https://doi.org/10.5194/npg26912019, 2019a.
Schulte, J. A.: R software package, available at: http://justinschulte.com/wavelets/advbiwavelet.html (last access: 22 November 2019), 2019b.
Schulte, J. A., Duffy, C., and Najjar, R. G.: Geometric and topological approaches to significance testing in wavelet analysis, Nonlin. Processes Geophys., 22, 139–156, https://doi.org/10.5194/npg221392015, 2015.
Shen, X. and Kimoto, M.: Influence of El Niño on the 1997 Indian summer monsoon, J. Meteor. Soc. Japan, 77, 1023–1037, 1999.
Shukla, J. and Paolino, D. A.: The Southern Oscillation and longrange forecasting of the summer monsoon rainfall over India, Mon. Weather Rev., 111, 1830–1837, 1983.
Slingo, J. M. and Annamalai, H.: 1997: The El Niño of the century and the response of the Indian summer monsoon, Mon. Weather Rev., 128, 1778–1797, 2000.
Su, J., Zhang, R., Li, T., Rong, X., Kug, J.S., and Hong, C.C.: Causes of the El Niño and La Niña amplitude asymmetry in the equatorial eastern Pacific, J. Climate, 23, 605–617, https://doi.org/10.1175/2009JCLI2894.1, 2010.
Timmermann, A.: Decadal ENSO amplitude modulations: A nonlinear mechanism, Global Planet. Change, 37, 135–156, 2003.
Torrence, C. and Compo, G. P.: A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61–78, 1998.
Torrence, C. and Webster, P. J.: Interdecadal changes in the ENSOmonsoon system, J. Climate, 12, 2679–2690, 1999.
Van Milligen, B. P., Saìnchez, E., Estrada, T., Hidalgo, C., Brañas,B., Carreras, B., and Garcìa, L.: Wavelet Bicoherence: A New Turbulence Analysis Tool, Phys. Plasmas, 2, 3017–3032, 1995.
van Oldenborgh, G. J. and Burgers, G.: Searching for decadal variations in ENSO precipitation teleconnections, Geophys Res Lett., 32, L15701, https://doi.org/10.1029/2005GL023110, 2005.
Walker, G. T. and Bliss, E. W.: World weather V, Mem., R. Meteorol. Soc., 4, 53–84, 1932.
Wu, A. and Hsieh, W. W.: Nonlinear interdecadal changes of the El NiñoSouthern Oscillation, Clim. Dynam., 21, 719–730, 2003.
Yun, K. S and Timmermann, A.: Decadal monsoonENSO relationships reexamined, Geophys Res Lett., 45, 2014–2021, 2018.
Zhang, Q., Xu, C., Jiang, T., and Wu, Y.: Possible influence of ENSO on annual maximum streamflow of the Yangtze River, China, J. Hydrol., 333, 265–274, https://doi.org/10.1016/j.jhydrol.2006.08.010, 2007.