**3. Materials**

The monthly rainfall (P) data of 8 rainfall observation stations and the monthly runoff (Q) data (long data from January 1984 to December 2015) at the hydrological station were collected from the Guizhou Provincial Hydrology and Water Resources Bureau (http://www.gzswj.gov.cn/hydrology\_gz\_new/index.phtml) (accessed on 15 September 2016). The monthly evaporation (E) data (January 1984 to December 2015), which were also derived from the Guizhou Provincial Hydrology and Water Resources Bureau, were measured by the evaporation dish of the hydrological observation station and represented the water evaporation of the water surface or soil. The monthly temperature (T) data

of the weather station corresponding to the runoff time series were obtained from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/) (accessed on 5 August 2021). The average annual rainfall data for the watershed were interpolated by the Kriging method with the annual rainfall data of all 8 rainfall stations (both inside and outside the watershed). The DEM data with a spatial resolution of 30 m were obtained from the International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences, which could be downloaded from the Geospatial Data Cloud (http://www.gscloud.cn) (accessed on 12 September 2021). The lithology data were derived from the Karst Scientific Data Center (http://www.karstdata.cn/) (accessed on 12 September 2021), Institute of Geochemistry, Chinese Academy of Sciences.

#### **4. Methodology**

A series of wavelet analysis methods were used to identify the multi-scale impact of climate factors on runoff change in the Yinjiang River watershed. These wavelet analysis methods included MRA [47], CWT [19–21,48,49], XWT [19–21], cross wavelet phase angle (CWPA) [50] and WTC [19]. MRA was carried out using a free MATLAB software package provided by the WaveLab Development Team and available at http://statweb.stanford. edu/~wavelab/, (accessed on 12 September 2021). Other methods, including CWT, XWT, WTC, and CWPA, were carried out using a free MATLAB software package (Mathworks, Natick, MA, USA) kindly provided by Grinsted et al. [19] at http://noc.ac.uk/usingscience/crosswavelet-wavelet-coherence, (accessed on 12 September 2021). The package includes code originally written by Torrence and Compo [20] of the University of Alaska, available at http://paos.colorado.edu/research/wavelets/, (accessed on 12 September 2021). The flowchart for identifying the multi-scale influences of climate factors on runoff changes is shown in Figure 2.

**Figure 2.** Flowchart for identifying the multi-scale influences of climate factors on runoff changes.

#### *4.1. MRA*

Choosing the particular values of a0 = 2 and s0 = 1 corresponds to the dyadic case used in MRA. The aim is to reduce/increase the resolution by a factor of 2 between two scales. Therefore, the approximation of a signal *x*(*t*) at a resolution *j*, denoted by *A<sup>j</sup> <sup>x</sup>*, and the detail of the same function at a resolution *j*, denoted by *D<sup>j</sup> <sup>x</sup>*, are defined by:

$$A\_{\mathcal{X}}^{\dot{j}}(t) = \sum\_{k=-\infty}^{+\infty} C\_{j,k} \psi\_{j,k}(t) \tag{1}$$

$$D\_x^j(t) = \sum\_{k=-\infty}^{+\infty} D\_{j,k} \Phi\_{j,k}(t) \tag{2}$$

where Φ*j*,*k*(*t*) is a scaled and translated basis function called the scaling function [47], which is determined with *ψj*,*k*(*t*) when a wavelet is selected. *Cj*,*<sup>k</sup>* is the scaling coefficient given the discrete sampled values of *x*(*t*) at resolution *j* and location *k*. It is calculated from Φ*j*,*k*(*t*) in a similar way for the wavelet coefficient *Dj*,*<sup>k</sup>* from *ψj*,*k*(*t*) for detailed mathematical expressions [49].

The signal *x*(*t*) can be reconstructed from the approximation and detail components as:

$$\mathbf{x}(t) = A\_x^j(t) + \sum\_{j=1}^{J} D\_x^j \tag{3}$$

where *J* is the highest resolution level considered. Since the MRA ensures variance is well captured in a limited number of resolution levels, the analysis of energy distribution in the sampling time series across scales gives a good idea of the energy distribution across frequencies.
