1 A Skewed perspective of the Indian rainfall-ENSO Relationship 1

7 Wavelet coherence is a commonly used method in hydrology to extract scale-dependent, non-stationary relationships 8 between time series. However, we show that the method cannot always determine why the time-domain correlation 9 between two time series temporally changes. We show that even for stationary coherence, the time-domain correlation 10 between two time series weakens if at least one of the time series has changing nonlinear characteristics. To overcome 11 this drawback, a nonlinear coherence method is proposed for quantifying the cross-correlation between nonlinear 12 modes embedded in time series. It is shown that using nonlinear coherence spectra together with auto-bicoherence 13 spectra can provide additional insight into changing time-domain correlations. The new method is applied to the El 14 Niño /Southern Oscillation (ENSO) and All-India rainfall, which is closely linked to hydrological processes across 15 the India sub-continent. The nonlinear coherence analysis showed that the skewness of All-India rainfall is weakly 16 correlated with that of 4 ENSO time series after the 1970s, indicating that increases in ENSO skewness after the 1970s 17 at least partially contributed to the weakening All-India rainfall-ENSO relationship in recent decades. The implication 18 of this result is that the intensity of skewed El Niño events is likely to overestimate India drought severity, which was 19 the case in the 1997 monsoon season, a time point when the nonlinear wavelet coherence between All-India rainfall 20 and ENSO reached its lowest value in the 1871-2016 period. We determined that the association between the 21 weakening ENSO-All India rainfall relationship and ENSO nonlinearity could reflect the contribution of different 22 nonlinear ENSO modes to ENSO diversity. 23


24
The South Asian Monsoon, which is the dominant source of precipitation source for the Indian subcontinent, 25 has been a target for seasonal prediction for well over a century (Blanford, 1884). Despite this long heritage of 26 research, skillful prediction remains a challenge, driving extensive and ongoing research on statistically and 27 dynamically based prediction methods (e.g., REFS). It is difficult to overstate the importance of the South Asian

28
Monsoon to the well-being of citizens in India. Strong monsoon years have caused are associated with catastrophic 29 flooding (Kale, 2012;Sanyal and Lu, 2005) and large landslides (Dortch et al., 2009), while weak monsoons have led 30 to water shortages (Mishra et al., 2016) and crop losses (Parthasarathy et al., 1988;Prasanna, 2014) that, in historical 31 times, were known to resulted in significant food shortages in historical times (Fagan, 2009). Thus, while the majority 32 of monsoon forecast studies target prediction of rainfall totals, the hydrological and agricultural impacts of monsoon 33 variability provide the most pressing motivation for the work.

71
The nonlinear characteristics (e.g. skewness) of ENSO are also non-stationary and undergo interdecadal 72 changes (Wu and Hsieh, 2003). Numerous studies have reported an ENSO regime shift in the 1970s in which ENSO 73 began to evolve more nonlinearly than in previous decades (An, 2004;An and Jin 2004;An, 2009). It is a curious fact 74 that the ENSO regime shift of the 1970s coincided with the weakening ENSO-Indian rainfall relationship as 75 documented by Kumar et al. (1999). This observation begs the question as to whether nonlinear ENSO regime changes 76 are related to changes in the ENSO-Indian rainfall relationship.

77
Various mechanisms have been proposed for explaining the cause of ENSO skewness. Kang and Kug (2002) 78 suggested that the asymmetry between the magnitude of El Niño and La Niña events is related to the relative westward 79 displacement of zonal wind stress anomalies during La Niña events compared to El Niño events. Jin et al., (2003) and 80 An and Jin (2004) found that ENSO asymmetry is related to nonlinear dynamical heating (NDH), where the magnitude 81 of NDH is related to the propagation characteristics of ENSO. As shown by An and Jin (2004), NDH during strong 82 El Niño events like the 1982/1983 and 1997/1998 events tends to be stronger than that of weak El Niño events because 83 SST anomalies tend to propagate eastward. Since the late 1970s there has been a propensity for eastward propagation 84 characteristics of ENSO (Santoso et al., 2013), contrasting with the time period before the 1970s that consisted of the 85 relatively weak El Niño events of 1957Niño events of /1958Niño events of and 1972Niño events of /1973 (An and Jin, 2004;An, 2009). More recently, Su et al.

86
(2010) showed that vertical temperature advection may have an opposing effect on ENSO asymmetry and that the 87 asymmetry in the extreme eastern equatorial Pacific is related to meridional ocean temperature advection. is not a bounded quantity and so high bi-spectral power does not always mean strong phase dependence.

108
In this study, the deficiencies associated with the above-mentioned techniques are addressed using higher-109 order wavelet analysis, which allows for the quantification of frequency-dependent and non-stationary nonlinearities

114
Indian rainfall relationship in recent decades is related to the shift of ENSO from a linear regime to a nonlinear one.

115
The paper is organized as follows: In Section 2, data used are described. Section 3 includes the description of the 116 implemented methodologies. Results are presented in Section 4 and concluding remarks are provided in Section 5.

136
The monthly SST data from 1871-2016 were based on the Hadley Centre Global Sea Ice and Sea Surface 137 Temperature (HadISST1; Rayner et al., 2003) The data at each grid point were converted to monthly anomalies in the 138 same way as they were computed for the ENSO and AIR All-India time series.  To better diagnose changes in time series statistics associated with AIR and ENSO, we adopted a wavelet 142 analysis. For a time series, X, comprising data points 1 , 2 , … , , the continuous wavelet transform is given by (1)

144
where s is wavelet scale, 0 is an analyzing wavelet, is a time step (1 month in this study), and n is time. The  Linear wavelet coherence (Table 1) was used to quantify the linear relationship between two time series as a 151 function of frequency and time. Linear wavelet coherence between two time series X and Y is given by

156
In the context of the Indian monsoon, strong coherence between rainfall and a climate pattern (e.g. ENSO) 157 at a scale s indicates shared temporal characteristics between a climate pattern and rainfall. Because theory supports a  proportional to the input (King, 1998). As ENSO is nonlinear, we adopted higher-order wavelet methods to address 167 the deficiencies of traditional wavelet methods.

168
The type of nonlinearities considered in this study were quadratic nonlinearities in which the scales 1 , 2 , 169 and 3 satisfied the sum rule

173
These types of nonlinearities arise, for example, when a sinusoid is squared, in which case a harmonic is produced.

174
More generally, quadratic nonlinearities induce time series skewness, which was computed in this study using

176
where S in the standard deviation of the time series and ̅ is the mean of the time series. Positive (negative) skewness 177 meant that the right (left) tail of the time series distribution is longer than the left (right). In other words, positive 5 (negative) skewness meant that there was a tendency for positive (negative) time series events (i.e. anomalies) to be 179 more intense than negative (positive) ones.

180
In this paper, quadratic nonlinearities giving rise to time series skewness were quantified using local and 181 global wavelet-based auto-bicoherence methods (Schulte, 2016b). Global auto-bicoherence (Table 1) was computed 182 using the equation

183
( 1 , 2 ) = 186 is the global bi-spectrum and the hat denotes the complex conjugate. Identical to wavelet coherence, auto-bicoherence 187 is bounded by 0 and 1, a value of 1 indicating the strongest possible phase coupling among the phases ( 3 ), ( 2 ),

188
and ( 1 ) such that sum rules Eq. (3-4) is satisfied. A peak in the auto-coherence spectrum at ( 1 , 2 ) indicated that 189 there was quadratic phase coupling among oscillatory modes with scales 1 , 2 , and 3 so the oscillatory modes was 190 contributing to time series skewness. It is important to note that the auto-bicoherence method cannot detect other types 191 of nonlinearities such as cubic nonlinearities whose detection would require trisepctra (Collis et al., 1998).

192
To determine if the strength of the quadratic phase coupling was a function of time changed temporally, the 193 local auto-bicoherence spectrum (Schulte, 2016b) given by was computed, where ( 1 , 1 ) is the local bi-spectrum given as

206
To be consistent with the wavelet power and coherence analyses, results for the higher-order wavelet analysis 207 were casted in terms of Fourier period rather than wavelet scale. The Fourier period corresponding to was denoted 208 by , where the Fourier period is obtained by multiplying by 1.03 for the Morlet wavelet (Torrence and Compo, 209 1998). Thus, the local diagonal slice of the auto-bicoherence spectra were plotted using the Fourier period 1 210 corresponding to 1 as the vertical axis and time as the horizonal axis.

212
The statistical significance of all wavelet spectra was evaluated using the cumulative area-wise test (

225
To assess the statistical significance of the global auto-bicoherence estimates, a modified version of the 226 cumulative area-wise test was applied. In the modified version of the cumulative area-wise test, the normalized area 227 of patches was computed by dividing patch area by the product 1 2 , where 1 is the mean first-coordinate of the patch 228 and 2 is the mean second coordinate. The reason for this modified normalized area is that dividing area by say, 1 2 , 229 retained the correlation between normalized area and 2 . The test was applied using the same point-wise significance 230 levels that were used to assess the statistical significance of wavelet power and coherence.  coherence measures the local cross-correlation between the skewness of 1 + 2 + 3 and 1 + 2 + 3 or between 257 1 + 1 /2 and 1 + 1 /2 in the case that 1 = 2 .

258
To better understand nonlinear coherence, we supposed that

261
( 3 ) − ( 3 ) = 3 (14) 262 for constants 1 , 2 , and 3 . Adding Eqs. (12) and (13) and subtracting Eq. (14) from the result produced the equality 265 for some constant K = 1 + 2 − 3 . Thus, if X was found to be perfectly nonlinear coherent with Y, then X and Y must 266 be were perfectly coherent at the three scales participating in the quadratic phase coupling. Even if the coherence was 267 perfect at two scales, the relative bi-phase ( 1 , 2 ) will fluctuated randomly if the relative phase difference at the 268 remaining scale fluctuated randomly so that the nonlinear coherence was will be low. Thus, if nonlinear coherence 269 was high, then there was some non-random relationship between X and Y at all three scales even if high linear 270 coherence was not identified at one or more scales. This theoretical idea indicates that nonlinear coherence can uncover 271 relationships that linear coherence cannot (see Figure S1 in supplementary material).

272
The relative bi-phase difference ( 1 , 2 ) is the higher-order analog of the relative phase difference

279
In this paper, we focused on nonlinear coherence computed along the diagonal slices ( 1 = 2 ) of the time 280 series bi-spectra. The nonlinear coherence spectra was then plotted using 1 as the vertical axis and time as the 281 horizonal axis. High nonlinear coherence at 1 and n meant that the skewness or asymmetry between 1 + 1 /2 and 282 1 + 1 /2 were locally cross-correlated.

283
To demonstrate the concept of nonlinear coherence, we considered a simple example in which the nonlinear 284 climate forcing time series was given by 285 286 and the response to the forcing was given as 287   (Figure 1c).

306
An inspection of the local auto-bicoherence spectrum of F(t) (Figure 2b) reveals that the auto-bicoherence at 307 1 = 32 is increasing with time, indicating that the phase coupling between modes with periods 3 = 16 and 1 = 32 308 is strengthening with time. The bi-phase of 0•, as indicated by arrows pointing to the right, confirms that the quadratic 309 phase coupling is contributing to the positive skewness seen in Figure 1a to an increasing degree. Furthermore, the 310 nonlinear coherence between R(t) and F(t) is weak and mostly statistically insignificant at 3 = 32 (Figure 2c

324
The more nonlinear F(t) becomes, the more F(t) will overestimate R(t) when F(t) is positively skewed. Similarly, the 325 positive forcing produces a negative response at t = 450 because of skewness and not simply a change in variance.

326
Nonlinear coherence allows for the quantification and identification of these time-domain aberrations.

327
The weakening relationship shown in Figure 1b could lead a researcher studying a hydrological process to 328 believe that another direct forcing mechanism must be influencing the hydrological process. This belief could lead to

362
The differences in skewness shown in Figure 4 suggests that the correlation between the ENSO time series

412
To confirm that the nonlinear phase coupling identified in Figure 8 is associated with skewed waveforms, we 413 inspected the corresponding local bi-phase spectra (not shown). It was found that the bi-phase in the 42-to 64-month 414 period band is generally 0° so that the nonlinear phase coupling in that period band contributes to the positive skewness

417
The results shown in Figure 9 indicate that the nonlinear wavelet coherence between AIR and the time series

429
However, unlike the theoretical example shown in Figure 2, the linear coherence between the ENSO time series and 430 AIR also weakens around the 1990s (Figure 7) so that the weakening relationship could be the result of a combination 431 of factors that includes ENSO nonlinearity.

432
The 20-year sliding mean of the ENSO auto-bicoherence, coherence, and nonlinear coherence averaged in 433 the 32 to 64-month period band further highlights the impact of ENSO nonlinearity. As shown in Figure 10a, the 434 sliding mean nonlinear coherence between the Niño 1+2 index and AIR fluctuates less than linear coherence and 435 reaches a clear global maximum around the 1970s before rapidly declining to a global minimum around the late 1990s 436 when the Niño 1+2 index is very nonlinear. As shown in Figure 10a, the Niño 1+2 auto-bicoherence peaks around the nonlinearity is only an important contributor when the timeseries is highly nonlinear, which is not the case for the 445 Niño 4 index because of the low auto-bicoherence (Figures 8b and 10b). Because the nonlinear coherence between 446 AIR and indices for the Niño 1+2 and Niño 4 is weak (Figures 9 and 10), the more pronounced change in the August-

447
September AIR-Nino 1+2 correlation reflects the more intense increase in Niño 1+2 nonlinearity compared to that of 448 the Niño 4 index in recent decades.

466
Although An and Jin (2004) and Burgers and Stephenson (1999) showed that skewness is greatest across the eastern 467 equatorial Pacific, we determined that such a time-domain approach is unable to capture frequency-dependent patterns 468 in nonlinearity.

469
The spatial auto-bicoherence plots associated with the peaks in the Niño 1+2 auto-bicoherence spectrum are 470 shown in Figure 13.

12
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, 490 but the degree to which the time series are nonlinear are different ( Figure 11)

506
-AIR contains periodicities ( Figure S2), which means that AIR is not simply a stochastically perturbed ENSO signal, 507 as noise does not contain periodicities. The retention of non-Gaussian noise features was certainly the case for R(t) -508 F(t) in the example in Section 3.5 because the difference retains the cosine function with a period of 16.

509
The fact that nonlinear coherence between rainfall and ENSO is determined by linear coherence between 510 ENSO and rainfall at two or three frequencies means that the changing time-domain correlation could be more fully 511 understood by determining why linear coherence changes at the frequencies that contribute to ENSO skewness. Such 512 an analysis could provide a more mechanistic perspective than the theoretical perspective adopted in this study

554
The third-order poly spectrum is the bi-spectrum, and the fourth-order poly spectrum is the tri-spectrum (Collis et al.,