*4.2. CWT*

The wavelet transform can be seen as a bandpass filter of uniform shape and varying location and width [20]. The continuous wavelet transform (CWT) *Wx*(*a*, *τ*) of a time series *x*(*t*) is given as follows:

$$\mathcal{W}\_{\mathbf{x}}(a,\tau) = \int\_{-\infty}^{+\infty} \mathbf{x}(t) \boldsymbol{\uppsi}^\*(t; a, \tau) d\_{\mathbf{f}} \tag{4}$$

$$
\psi(t;a,\tau) = \frac{1}{\sqrt{a}} \psi(\frac{t-\tau}{a}) \tag{5}
$$

$$
\psi(t) = \sqrt[4]{\frac{1}{\pi}} \cos(kt) e^{-\frac{t^2}{2}} \tag{6}
$$

where *Wx*(*a*, *τ*) represents a group of wavelet functions, *W*(*a*,*b*), based on a mother wavelet *ψ* which can be scaled and translated, modifying the scale parameter a and the translation parameter *τ,* respectively. *ψ*∗(*t*; *a*, *τ*) corresponds to the complex conjugate of *ψ*(*t*; *a*, *τ*). *ψ*(*t*) is the Morlet wavelet function, *k* is the non-dimensional frequency, here taken to be 6 to satisfy the admissibility condition, and *t* is time.

The statistical significance of wavelet power can be assessed relative to the null hypothesis that the signal is generated by a stationary process with a given background power spectrum (*Pη*). Many geophysical time series have distinctive red noise characteristics that can be modeled very well by a first-order autoregressive (AR1) process. The Fourier power spectrum of an AR1 process with lag-1 autocorrelation *α* [48] is given by

$$P\_{\eta} = \frac{1 - a^2}{\left| 1 - \alpha e^{-2i\pi\eta} \right|^2} \tag{7}$$

The probability that the wavelet power of a process with a given power spectrum (*Pη*) [19] is greater than *p* is

$$D\left(\frac{\left|W\_{\eta}^{X}(s)\right|^{2}}{\sigma\_{X}^{2}} < p\right) = \frac{1}{2}P\_{\eta}\chi\_{v}^{2}(p) \tag{8}$$

where *η* is the Fourier frequency index. *ν* is equal to 1 for real and 2 for complex wavelets.

*4.3. XWT*

The XWT of two time series *X* and *Y* is defined as

$$\mathcal{W}^{XY} = \mathcal{W}^X \mathcal{W}^{Y^\*} \tag{9}$$

where \* denotes complex conjugation. We further define the cross wavelet power as |*WXY*|. The complex argument arg(*WXY*) can be interpreted as the local relative phase between *X* and *Y* in time–frequency space. The theoretical distribution of the cross wavelet power of two time series with background power spectra *P<sup>X</sup> <sup>k</sup>* and *<sup>P</sup><sup>Y</sup> <sup>k</sup>* is given as

$$D\left(\frac{\mathcal{W}\_{\boldsymbol{n}}^{X}(\boldsymbol{s})\mathcal{W}\_{\boldsymbol{n}}^{Y}(\boldsymbol{s})}{\sigma\_{X}\sigma\_{Y}} < p\right) = \frac{Z\_{\boldsymbol{v}}(p)}{\boldsymbol{v}}\sqrt{P\_{k}^{X}P\_{k}^{Y}}\tag{10}$$

where *Zv*(*p*) is the confidence level associated with the probability *p* defined by the square root of the product of two X2 distributions [19].

We use the circular mean of the phase over regions with higher than 5% statistical significance that are outside the cone of influence (COI) to quantify the phase relationship. This is a useful and general method for calculating the mean phase. The circular mean of a set of angles (*ai*, *i* =1... *n*) is defined as

$$a\_m = \arg(X, Y) \text{ with } X = \sum\_{i=1}^{n} \cos(a\_i) \text{ and } Y = \sum\_{i=1}^{n} \sin(a\_i) \tag{11}$$

It is difficult to calculate the confidence interval of the mean angle reliably since the phase angles are not independent. The number of angles used in the calculation can be set arbitrarily high simply by increasing the scale resolution. However, it is interesting to know the scatter of angles around the mean. For this, we define the circular standard deviation as

$$S = \sqrt{-2\ln(R/n)}\tag{12}$$

where *<sup>R</sup>* <sup>=</sup> <sup>√</sup>*X*<sup>2</sup> <sup>+</sup> *<sup>Y</sup>*2. The circular standard deviation is analogous to the linear standard deviation in that it varies from zero to infinity.

## *4.4. WTC*

Cross wavelet power reveals areas with high common power. Another useful measure is how coherent the cross wavelet transform is in time–frequency space. We define the wavelet coherence of two time series as

$$R\_n^2(s) = \frac{\left|\mathcal{S}(s^{-1}\mathcal{W}\_n^{XY}(s))\right|^2}{\mathcal{S}(s^{-1}\left|\mathcal{W}\_n^X(s)\right|^2)\cdot\mathcal{S}(s^{-1}\left|\mathcal{W}\_n^Y(s)\right|^2)}\tag{13}$$

where *S* is a smoothing operator. Notice that this definition closely resembles that of a traditional correlation coefficient, and it is useful to think of the wavelet coherence as a localized correlation coefficient in time–frequency space. We write the smoothing operator *S* as

$$S(\mathcal{W}) = S\_{scale}(S\_{time}(\mathcal{W}\_n(s)))\tag{14}$$

where *S* scale denotes smoothing along the wavelet scale axis and *S* time denotes smoothing in time. It is natural to design the smoothing operator so that it has a similar footprint as the wavelet used. For the Morlet wavelet, a suitable smoothing operator is given by Torrence and Compo [20]:

$$S\_{time}\left(\mathcal{W}\right)|\_s = \left. \left( \mathcal{W}\_n(s) \cdot c\_1^{\frac{-\rho^2}{2s^2}} \right) \right|\_s \tag{15}$$

$$\left.S\_{time}\left(W\right)\right|\_{s} = \left.\left(\mathcal{W}\_{n}\left(\mathbf{s}\right)\cdot\mathbf{c}\_{2}\left.\Pi\left(\mathbf{0}.6\mathbf{s}\right)\right)\right|\_{n} \tag{16}$$

where *c*<sup>1</sup> and *c*<sup>2</sup> are normalization constants and Π is the rectangle function. The factor of 0.6 is the empirically determined scale decorrelation length for the Morlet wavelet [20]. In practice, both convolutions are conducted discretely and therefore the normalization coefficients are determined numerically.

In this study, the cross wavelet energy, wavelet correlation agglomeration and phase spectra were calculated for monthly temperature, evaporation and rainfall and monthly runoff series to analyze the multi-temporal correlation amongst temperature, evaporation, rainfall and runoff. The hydrometeorological variables were used as input and output signals to characterize the responses of runoff changes to climatic factors in the Yinjiang River watershed. Climatic factors (rainfall, temperature and evaporation) were taken as input signals, and runoff was taken as an output signal. The correlation between climatic factors and runoff signals in the frequency and time domains at different energies was analyzed by using XWT and WTC. We had focused our analysis on the P–Q relationship to assess the impact of rainfall on runoff changes. Besides that, the E–Q and T–Q relationship were used to assess the impacts of evaporation and temperature on runoff changes, respectively.

#### **5. Result Analysis**
