Runoff Dynamics and Associated Multi-Scale Responses to Climate Changes in the Middle Reach of the Yarlung Zangbo River Basin , China

Long-term hydro-climatic datasets and sophisticated change detection methods are essential for estimating hydro-climatic trends at regional and global scales. Here, we use the ensemble empirical mode decomposition method (EEMD) to investigate runoff oscillations at different time scales and its response to climatic fluctuations in the middle of the Yarlung Zangbo River Basin (MYZRB) over the period 1961–2009. In the study region, results revealed that the runoff presented an overall nonlinear and nonstationary decreasing-increasing alternative trend with weak quasi-three-year and unobvious quasi-five-year cycles at the inter-annual scale, while, significance was discovered with quasi-12-year and quasi-46-year cycles at the inter-decadal scale. Variance contribution rates of the hydrological components suggested that the inter-annual oscillations played an essential role in the runoff variations in the MYZRB. According to the reconstructed inter-decadal runoff series, the runoff may keep declining in future. For the response of runoff to climate change, overall, the runoff had a positive correlation with precipitation and a negative correlation with extreme temperature. But the runoff did not show obvious correlation with mean temperature. Furthermore, from a temporal scale point of view, the inter-annual runoff showed significant response to the inter-annual precipitation. The inter-decadal runoff strongly responded to the inter-annual extreme temperature. These findings will help us understand the hydro-climatic intrinsic mechanism in the MYZRB and develop better water resources management to account for climate change impact.


Introduction
The hydrologic cycle is one of the major components of the planetary system regulating human, animal, and plant life [1].It is desirable for an overall balance among the hydrologic cycle components (i.e., precipitation, evapotranspiration, runoff, infiltration, and glaciers) at the basin-scale to sustain the ecosystem functions [2].However, these hydrologic cycle components dynamically interact with each other to a certain extent, which apparently correlates with the effects of climate change [2,3].Following the third and fourth assessment reports of the Intergovernmental Panel on Climate Change (IPCC), the fifth assessment report pointed out that the global warming is in evidence [4].Climate warming can greatly change the hydrologic cycle and redistribute water resources over time and space [5].The Tibetan Plateau is considered as a driver and amplifier of global climatic changes and is also the headwater region of many rivers.In the Tibetan Plateau, the change to runoff processes caused by climate change is directly related to the formation and evolution of local water resources and posed profound impacts on the water resources of its middle and low reaches.In the Tibet Autonomous Region, the Yarlung Zangbo River Basin (YZRB) is the largest river basin, and in recent years, the issues of climate change, glacier retreat, and hydrologic alternation in this region have drawn much attention [5][6][7].In the context of global warming, temperature and precipitation have shown an increasing trend since 1961 [6,8], while glacier accumulation and potential evapotranspiration underwent a decrease trend [6,7].Among them, the changes to precipitation and total water resources were even more significant due to agricultural and industrial development in YZRB [9], especially in the MYZRB.The MYZRB represents the agricultural, industrial, political, and cultural center of Tibet, with several large cities (Lhasa, Shigatse, Shannan, and Nyingchi) located there, and over half of the Tibetan population lives there.Therefore, studying the hydrological cycle of MYZRB and its response to climate change in the context of global warming has important practical significance and scientific value.
Hydro-climatic change detection plays an important role in accurately estimating global and regional hydro-climatic change trends and understanding their complex mechanism [10].Most hydro-climatic variables (e.g., runoff, precipitation, and temperature) with long-term variations illustrate nonlinear and nonstationary complex processes accompanied by different scales of periodic oscillations [11].Although studies for investigating the hydrologic and climatic changes over the YZRB are available in the literature, few can be found for exploring the nonlinear runoff changes and associated multi-scale responses to climate fluctuations, especially in the MYZRB.Some conventional change detection methods such as moving average (MA), linear regression (LR), Mann-Kendall trend (MK), empirical orthogonal functions (EOF), and artificial neural network (ANN) were used to analyze the change trends of hydro-climatic factors [6,[8][9][10].However, their real fluctuation still remains unclear [5].A new signal detection technology, EEMD, developed on the basis of empirical mode decomposition (EMD) [12] with stronger self-adaptability and local variation characteristics to overcome the scale-mixing problem through assist of white noise and multiple mean of an ensemble of trials [13,14], has been gradually applied in different regions to detect the fluctuations and oscillations of hydro-climatic changes for exploring some significant characteristics [15][16][17].In the present study, the EEMD was adopted to reveal the nonlinear change characteristics of annual runoff and to investigate the runoff responses to climatic changes at inter-annual and inter-decadal scales in the MYZRB.
This study focused on runoff dynamics and associated multi-scale responses to climate change in the MYZRB with statistical analyses in three aspects: (1) characterizing the oscillation and dynamics of runoff at different time scales, (2) investigating the contributions of oscillations at different time scales to runoff dynamics and their significance level, and (3) projecting the response of runoff to mean and extreme of climate changes.The rest of this paper is organized as follows: Section 2 describes the study area, data utilized, and methodology, respectively.Section 3 presents the results.Section 4 provides discussion and is followed by the conclusions given in Section 5.

Study Area
The YZRB is the highest river basin in the world with its average elevation higher than 4600 m a.s.l.[9].A certain amount of bare rocks and glaciers occupy of 11.8% and 1.2% of the total area of the basin developed in this watershed.MYZRB, as shown in Figure 1, is located between 27 • 49 -31 • 17 N and 85 • 23 -93 • 21 E and is gauged mainly by the Lazi hydrometric station at the inlet and the Yangcun hydrometric station at the outlet of MYZRB [18].The drainage of the MYZRB accounts for 47.9% of the YZRB in area within China, approximately about 118,000 km 2 .The MYZRB, situated in the monsoon climate zone of the temperate semi-humid and semi-arid plateau, characterizes strong radiation, long sunshine time, low average temperature, and low accumulated temperature.The moisture from Bengal Bay transported by Indian monsoons and Western Pacific subtropical high from the east to the west crossing the Yarlung Zangbo Suture Zone constitutes the primary source for precipitation of the area.Affected by Indian monsoons with thier abundant moisture and high temperature, the MYZRB receives the most precipitation from June to September and generates its maximum discharge in the same duration.
There are only 16 national meteorological stations around and within the study area.Thus, it is difficult to make an accurate analysis for climate change only with these stationary observations.Given this problem, we adopt a more detailed gridded climate datasets (Figure 1) produced from intensive meteorological stations to investigate the temperature and precipitation changes in the MYZRB.
from the east to the west crossing the Yarlung Zangbo Suture Zone constitutes the primary source for precipitation of the area.Affected by Indian monsoons with thier abundant moisture and high temperature, the MYZRB receives the most precipitation from June to September and generates its maximum discharge in the same duration.
There are only 16 national meteorological stations around and within the study area.Thus, it is difficult to make an accurate analysis for climate change only with these stationary observations.Given this problem, we adopt a more detailed gridded climate datasets (Figure 1) produced from intensive meteorological stations to investigate the temperature and precipitation changes in the MYZRB.

Dataset Utilized
The China daily Temperature/Precipitation Analysis Product (CTAP and CPAP, respectively) with a 0.5° spatial resolution derived from the National Meteorological Information Center (NMIC) and China Meteorological Administration (CMA) were used in this study.Based on the observations from 2472 national meteorological stations, the CTAP and CPAP were produced by the Thin Plate Spline (TPS) and climatology-based Optimal Interpolation (OI) technique, respectively [8,19].All the observations used in CTAP and CPAP have undergone strict quality control (QC), including the extrema detection, internal consistency check, and spatial consistency check, before the interpolation operation [20,21].It is worth mentioning that snowfall was measured accurately by a heated device attached to rain gauges through the manual heating method [21].To further improve the accuracy of CPAP, the Parameter-Elevation Regression on Independent Slopes Model (PRISM) was adopted to correct orographic errors [22].Previous studies on the validation of accuracy for CPAP using stationary observations showed that relative bias was around 3.21% and the absolute error of 91% observations was less than 1.0 mm/d [23][24][25].The CPAP has been successfully utilized for validating high spatio-temporal resolution satellite-based precipitation products in recent years [22,23,26,27].For CTAP, including maximum (Tmax), minimum (Tmin), and mean temperature (Tmean) datasets, the average deviations are mainly concentrated in the ±0.2 °C and the RMSEs fluctuate around 0.25 °C, 0.35 °C, and 0.25 °C, respectively [28].The CTAP and CPAP were selected in this study to obtain temperature and precipitation indices in the MYZRB over the period 1961-2009.With these climate

Dataset Utilized
The China daily Temperature/Precipitation Analysis Product (CTAP and CPAP, respectively) with a 0.5 • spatial resolution derived from the National Meteorological Information Center (NMIC) and China Meteorological Administration (CMA) were used in this study.Based on the observations from 2472 national meteorological stations, the CTAP and CPAP were produced by the Thin Plate Spline (TPS) and climatology-based Optimal Interpolation (OI) technique, respectively [8,19].All the observations used in CTAP and CPAP have undergone strict quality control (QC), including the extrema detection, internal consistency check, and spatial consistency check, before the interpolation operation [20,21].It is worth mentioning that snowfall was measured accurately by a heated device attached to rain gauges through the manual heating method [21].To further improve the accuracy of CPAP, the Parameter-Elevation Regression on Independent Slopes Model (PRISM) was adopted to correct orographic errors [22].Previous studies on the validation of accuracy for CPAP using stationary observations showed that relative bias was around 3.21% and the absolute error of 91% observations was less than 1.0 mm/d [23][24][25].The CPAP has been successfully utilized for validating high spatio-temporal resolution satellite-based precipitation products in recent years [22,23,26,27].For CTAP, including maximum (Tmax), minimum (Tmin), and mean temperature (Tmean) datasets, the average deviations are mainly concentrated in the ±0.2 • C and the RMSEs fluctuate around 0.25 • C, 0.35 • C, and 0.25 • C, respectively [28].The CTAP and CPAP were selected in this study to obtain temperature and precipitation indices in the MYZRB over the period 1961-2009.With these climate data, we obtained basin-scale annual mean precipitation (Pav) and average of annual Tmean (Tavg).Moreover, two basin-scale extreme temperature indices, including annual minimum Tmin (TNn) and annual maximum Tmax (TXx), were computed using the RClimDex1.1 (available on http://etccdi.pacificclimate.org/software.shtml).
The observed annual runoff at the national hydrological stations Lazi and Yangcun (Figure 1) over the period 1961 to 2009 were used in this study.These data have been subjected to strict quality control standard.The statistical information, including mean, standard deviation (Std.Deviation), Skewness, and Kurtosis, is illustrated in Table 1.With these datasets, we obtained the runoff yield of MYZRB by calculating the runoff difference between the Lazi and the Yangcun hydrological stations for analyzing the variability of the yielded runoff and its responses to climate changes over the MYZRB.The runoff in this study represents the runoff yield in the MYZRB, namely the runoff from upper reach was deducted.

Methodology
The EEMD, designed to provide a more accurate expression of a time series in a time-frequency-energy domain, is an adaptive and temporally local filter method and does not rely on any a priori determined basis functions.Several studies indicated that it is more suitable for analyzing nonlinear and non-stationary time series than many traditional decomposition methods (e.g., Fourier Transform, wavelet decomposition, and empirical orthogonal functions) [5,12,14,15].
In this study, the method was employed to display the oscillation cycles and trends of climate and runoff changes over the period 1961-2009.It decomposes a time series X into a finite number of intrinsic mode functions (IMFs) and a residual (RES) as Equation (1).
where n represents the number of the IMFs.The IMFs reveal the oscillation characteristics from high frequency (HF, less than 10-year period) to low frequency (LF, not less than 10-year period) at different time scales [5], and the RES reflects the trend of the original time series.For a more comprehensive review of the EEMD, please refer to Huang et al. [12] and Wu et al. [14].
A posteriori test method, as reported in [13,29,30], was employed to test the significance of the IMF components.Before that, we should obtain the distribution of the energy respecting the period in the form of spectral function [5].The energy density of the kth IMF denoted as E k can be calculated as Equation ( 2).
where N represents the lengths of IMF.Next, before examining the periods of IMFs, the properties of an IMF should be listed as follows: an IMF is any function having symmetric envelopes defined by the local maxima and minima separately and also having the same number of zero-crossings and extrema.Based on this definition, the mean period of the function can be computed through counting the number of peaks (local maxima) of the function [29].The mean period of the kth IMF denoted as T k can be described as Equation (3).
where NP k denotes the number of peaks for each IMF.Since the white noise sequence is tested by the Monte Carlo method [13,14], a series of T k and E k are generated.Thus, we can obtain a simple equation that relates the averaged energy density ( Ēk ) and the averaged period (T k ) as the Equation (4).
where α is the significance level (e.g., α = 0.05).Theoretically, for the IMF of the white noise series, the relation between the energy density and the averaged period can be expressed by a straight line with the slope of −1.However, there is little deviation in the actual applications.Therefore, the confidence limits for the energy spectrum distribution of white noise is listed as Equation ( 5) [29].
Through decomposition, if the energy of an IMF component is located above the confidence curve, indicating the periodic oscillation is statistically significant at corresponding level of α, or it is considered to contain less actual physical meaning in the specified IMF component [5,13,14].
In addition, it is noteworthy that the mirror-symmetric extension was used to solve the overshooting and undershooting phenomenon of the impact of the boundary on the decomposition process in EEMD [5,14].

Runoff Dynamic Analysis
As shown in Figure 2, an overall decreasing trend with an obvious turn point in the mid-1980s over the period 1961-2009 can be found in the annual runoff variations of MYZRB.The runoff in the MYZRB over the period 1965-1997 in average was relatively less than that in the other years, which implied that the MYZRB was relatively drier over that period.During the 2000s, the runoff in the MYZRB was more abundant, with an anomaly value up to about 189.83 × 10 8 m 3 , which suggested that an extreme hydrological event with a broken out frequency has occurred.As a whole, the runoff in the MYZRB presented a decreasing-increasing trend alternatively over the period 1961-2009, with the annual mean value of about 299.93 × 10 8 m 3 and the ratio between the maximum and the minimum annual runoff reached to 3.69:1.All these results suggest that the annual runoff in the MYZRB varied in nonlinear and nonstationary processes over the period 1961-2009, and this was supported by the five-year moving mean of runoff anomaly variations as exhibited in Figure 2.
After decomposing the original annual runoff series of the MYZRB over the period 1961-2009 with EEMD, four IMFs (IMF1-4) and one trend component RES were obtained, as shown in Figure 3. Certain level of inherent oscillation in the original runoff series and the significance level of the actual physical meaning thus were recognized from HF to LF at different time scales, and the trend of the original runoff series over time can be estimated from the RES [30].The periodic significances of IMFs were illustrated in Figure 4, which exhibited that the oscillation period of the decomposed components IMF1 and IMF2 cannot reach to 90% significance level, while IMF3 and IFM4 reach to above 90% confidence levels, implying that IMF3 and IMF4 were statistically of relatively more actual physical meaning.As shown in Figure 4, a weak quasi-three-year cycle appeared in IMF1 and an unobvious quasi-five-year cycle appeared in IMF2 at the inter-annual scale, while a significant quasi-12-year cycle appeared in IMF3 and a quasi-46-year year appeared in IMF4 at the inter-decadal scale.When combining Figures 3 and 4, we can find that the decomposed components such as IMF1-4 contained the nonlinear feedback of the hydrologic circulations and the periodic oscillations of hydrologic circulations under external forcing.The variance contribution rate of each IMFs for the original runoff series is listed in Table 2.It is worthwhile mentioning that although the cycles revealed in IMF1 and IMF2 were not significant, their variance contribution rates were higher compared with IMF3 and IMF4, which suggested that they cannot be ignored in maintaining the total energy of the original runoff signals.Combining Figure 3 with Table 2, one can find that the runoff oscillation characterizing an increasing-decreasing alternative trend in runoff change amplitude was clearly appeared in IMF1, and its variance contribution rate reached to the largest of 43.2%.Similar to the original runoff series, the runoff oscillation amplitudes over the period 1960-1964 and in 2000s were relatively higher than other periods.The variance contribution rate of IMF2 was about 16.7%, primarily reflecting that the runoff amplitude became larger after the turn point appeared in the mid-1980s.For inter-decadal signals of IMF3 and IMF4, their variance contribution rates were 18.6% and 8.0%, respectively, in which the former mainly reflected that the runoff oscillation amplitude was larger and experienced increasingdecreasing trend over the period 1997-2009 over the MYZRB.Although the lowest variance   The variance contribution rate of each IMFs for the original runoff series is listed in Table 2.It is worthwhile mentioning that although the cycles revealed in IMF1 and IMF2 were not significant, their variance contribution rates were higher compared with IMF3 and IMF4, which suggested that they cannot be ignored in maintaining the total energy of the original runoff signals.Combining Figure 3 with Table 2, one can find that the runoff oscillation characterizing an increasing-decreasing alternative trend in runoff change amplitude was clearly appeared in IMF1, and its variance contribution rate reached to the largest of 43.2%.Similar to the original runoff series, the runoff oscillation amplitudes over the period 1960-1964 and in 2000s were relatively higher than other periods.The variance contribution rate of IMF2 was about 16.7%, primarily reflecting that the runoff amplitude became larger after the turn point appeared in the mid-1980s.For inter-decadal signals of IMF3 and IMF4, their variance contribution rates were 18.6% and 8.0%, respectively, in which the former mainly reflected that the runoff oscillation amplitude was larger and experienced increasingdecreasing trend over the period 1997-2009 over the MYZRB.Although the lowest variance The variance contribution rate of each IMFs for the original runoff series is listed in Table 2.It is worthwhile mentioning that although the cycles revealed in IMF1 and IMF2 were not significant, their variance contribution rates were higher compared with IMF3 and IMF4, which suggested that they cannot be ignored in maintaining the total energy of the original runoff signals.Combining Figure 3 with Table 2, one can find that the runoff oscillation characterizing an increasing-decreasing alternative trend in runoff change amplitude was clearly appeared in IMF1, and its variance contribution rate reached to the largest of 43.2%.Similar to the original runoff series, the runoff oscillation amplitudes over the period 1960-1964 and in 2000s were relatively higher than other periods.The variance contribution rate of IMF2 was about 16.7%, primarily reflecting that the runoff amplitude became larger after the turn point appeared in the mid-1980s.For inter-decadal signals of IMF3 and IMF4, their variance contribution rates were 18.6% and 8.0%, respectively, in which the former mainly reflected that the runoff oscillation amplitude was larger and experienced increasing-decreasing trend over the period 1997-2009 over the MYZRB.Although the lowest variance contribution rate was found in IMF4 for original runoff series, the most significant periodic oscillation presented similar to cosine curve with about quasi-46-year cycle was detected.The RES component variations suggested a nonlinear decreasing-increasing trend over the period 1961-2009, with the turn point in the mid-1980s.As listed in Table 2, the variance contribution rates of each IMFs and RES revealed that the interannual oscillations were stronger than inter-decadal oscillations in runoff dynamics over the period 1961-2009 over the MYZRB.According to the theory of EEMD [13,14], we reconstructed the interannual and inter-decadal runoff series by IMFs and RES, in which the inter-annual runoff series was obtained by summing the IMF1, IMF2, and RES, and the inter-decadal runoff series was derived by summing the IMF3, IMF4, and RES.As seen in Figure 5, it was clear that the inter-annual runoff series was basically consistent with the original runoff anomaly series, which characterized the oscillation of original runoff anomaly series to a large extent, indicating the inter-annual oscillation playing an important role in the overall runoff dynamics over the MYZRB.However, the change trend was not very consistent between the inter-annual runoff series and original runoff anomaly series over the period 1997-2009, where a possible reason was that small-scale oscillations such as IMF1-2 contained lots of outside system noise [5].The inter-decadal runoff series effectively revealed the alternative appearances of wet and dry period over the period 1961-2009 over the MYZRB.In other words, the runoff series at the decadal scale featured an alternating positive and negative phase over the study period.According to the inter-decadal runoff series, we can infer that the runoff over MYZRB may present a reduction stage in future.In terms of RES, it can fully reflect the characteristic of a decreasing-increasing change trend over the period 1961-2009 over the MYZRB.As listed in Table 2, the variance contribution rates of each IMFs and RES revealed that the inter-annual oscillations were stronger than inter-decadal oscillations in runoff dynamics over the period 1961-2009 over the MYZRB.According to the theory of EEMD [13,14], we reconstructed the inter-annual and inter-decadal runoff series by IMFs and RES, in which the inter-annual runoff series was obtained by summing the IMF1, IMF2, and RES, and the inter-decadal runoff series was derived by summing the IMF3, IMF4, and RES.As seen in Figure 5, it was clear that the inter-annual runoff series was basically consistent with the original runoff anomaly series, which characterized the oscillation of original runoff anomaly series to a large extent, indicating the inter-annual oscillation playing an important role in the overall runoff dynamics over the MYZRB.However, the change trend was not very consistent between the inter-annual runoff series and original runoff anomaly series over the period 1997-2009, where a possible reason was that small-scale oscillations such as IMF1-2 contained lots of outside system noise [5].The inter-decadal runoff series effectively revealed the alternative appearances of wet and dry period over the period 1961-2009 over the MYZRB.In other words, the runoff series at the decadal scale featured an alternating positive and negative phase over the study period.According to the inter-decadal runoff series, we can infer that the runoff over MYZRB may present a reduction stage in future.In terms of RES, it can fully reflect the characteristic of a decreasing-increasing change trend over the period 1961-2009 over the MYZRB.

The Multi-Scale Response of Runoff Dynamics to Climate Changes
In the context of global warming, temperature and precipitation changes in a mountainous basin have a significant influence on runoff variation [11,31].We obtained four IMF (IMF1-4) and one trend component for the Pav anomaly series, Tavg anomaly series, TNn anomaly series, and TXx anomaly series over the same study period using EEMD.As shown in Figure 6a, the significance test results of Pav anomaly series for IMF1-4 reveals that precipitation has strong quasi-3-year period appeared in IMF1 and unobvious quasi-seven-year period (IMF2) at inter-annual scale, and unobvious quasi-12-year (IMF3) and quasi-44-year (IMF4) periods at inter-decadal scale.When comparing the IMFs between runoff and precipitation, we can find that the runoff variation shows a similar inter-annual scale but a different inter-decadal scale over the MYZRB.Similarly, as seen in Figure 6b, the Tavg has weak quasi-three-year period and unobvious five-year period at the inter-annual scale, and unobvious quasi-12-year and quasi-39-year periods at the inter-decadal periods.The TNn has weak quasi-3-year period and unobvious quasi-five-year period at the inter-annual scale, and unobvious quasi-12-year and quasi-32-year periods at the inter-decadal scale (Figure 6c).From Figure 6d, we can find that the TXx has weak quasi-three-year period and unobvious quasi-six-year period, and unobvious quasi-12-year and quasi-30-year periods.Therefore, it is noted that Tavg, TNn and TXx have a similar inter-annual scale to original runoff series but a different inter-decadal scale.
In combination with Figures 5 and 7a, we can find that the overall change trend (decreasingincreasing) of runoff is consistent with precipitation over the MYZRB.However, they have different turning points.The turning points in runoff and precipitation series are in the mid-1980s and the early 1980s, respectively.Moreover, after turning points the precipitation has a larger increasing trend than runoff, which is probably caused by continuously increasing temperature shown in Figure 7b-d; the turning points of Tavg, TNn, and TXx appear in the late 1980s.Higher temperature tends to raise the actual evapotranspiration (AET), resulting in mitigating the increasing amplitude of runoff [32].These results reveal that the runoff variation of the MYZRB over the period 1961-2009 has a certain response to the regional climate change, but there are some differences among their turning points, which indicate that this response has a certain time lag.

The Multi-Scale Response of Runoff Dynamics to Climate Changes
In the context of global warming, temperature and precipitation changes in a mountainous basin have a significant influence on runoff variation [11,31].We obtained four IMF (IMF1-4) and one trend component for the Pav anomaly series, Tavg anomaly series, TNn anomaly series, and TXx anomaly series over the same study period using EEMD.As shown in Figure 6a, the significance test results of Pav anomaly series for IMF1-4 reveals that precipitation has strong quasi-3-year period appeared in IMF1 and unobvious quasi-seven-year period (IMF2) at inter-annual scale, and unobvious quasi-12-year (IMF3) and quasi-44-year (IMF4) periods at inter-decadal scale.When comparing the IMFs between runoff and precipitation, we can find that the runoff variation shows a similar inter-annual scale but a different inter-decadal scale over the MYZRB.Similarly, as seen in Figure 6b, the Tavg has weak quasi-three-year period and unobvious five-year period at the inter-annual scale, and unobvious quasi-12-year and quasi-39-year periods at the inter-decadal periods.The TNn has weak quasi-3-year period and unobvious quasi-five-year period at the inter-annual scale, and unobvious quasi-12-year and quasi-32-year periods at the inter-decadal scale (Figure 6c).From Figure 6d, we can find that the TXx has weak quasi-three-year period and unobvious quasi-six-year period, and unobvious quasi-12-year and quasi-30-year periods.Therefore, it is noted that Tavg, TNn and TXx have a similar inter-annual scale to original runoff series but a different inter-decadal scale.
In combination with Figures 5 and 7a, we can find that the overall change trend (decreasing-increasing) of runoff is consistent with precipitation over the MYZRB.However, they have different turning points.The turning points in runoff and precipitation series are in the mid-1980s and the early 1980s, respectively.Moreover, after turning points the precipitation has a larger increasing trend than runoff, which is probably caused by continuously increasing temperature shown in Figure 7b-d; the turning points of Tavg, TNn, and TXx appear in the late 1980s.Higher temperature tends to raise the actual evapotranspiration (AET), resulting in mitigating the increasing amplitude of runoff [32].These results reveal that the runoff variation of the MYZRB over the period 1961-2009 has a certain response to the regional climate change, but there are some differences among their turning points, which indicate that this response has a certain time lag.In order to quantitatively describe the response of runoff to precipitation and temperature, we obtained the inter-annual and inter-decadal Pav, Tavg, TNn and TXx anomaly series using the same method as runoff, with IMFs and RES.Through a multi-scale correlation analysis, the correlation coefficients (CC) and associated significance level among annual runoff, Pav, Tavg, TNn and TXx over the period 1961-2009 over the MYZRB are listed in Table 3.It can be found that runoff has a positive correlation with precipitation at different time scales, but the most significant one is interannual precipitation versus inter-annual runoff, indicating that the t runoff variation is sensitive to response to precipitation in the inter-annual time scale.For the response of runoff to temperature,  In order to quantitatively describe the response of runoff to precipitation and temperature, we obtained the inter-annual and inter-decadal Pav, Tavg, TNn and TXx anomaly series using the same method as runoff, with IMFs and RES.Through a multi-scale correlation analysis, the correlation coefficients (CC) and associated significance level among annual runoff, Pav, Tavg, TNn and TXx over the period 1961-2009 over the MYZRB are listed in Table 3.It can be found that runoff has a positive correlation with precipitation at different time scales, but the most significant one is interannual precipitation versus inter-annual runoff, indicating that the t runoff variation is sensitive to response to precipitation in the inter-annual time scale.For the response of runoff to temperature, In order to quantitatively describe the response of runoff to precipitation and temperature, we obtained the inter-annual and inter-decadal Pav, Tavg, TNn and TXx anomaly series using the same method as runoff, with IMFs and RES.Through a multi-scale correlation analysis, the correlation coefficients (CC) and associated significance level among annual runoff, Pav, Tavg, TNn and TXx over the period 1961-2009 over the MYZRB are listed in Table 3.It can be found that runoff has a positive correlation with precipitation at different time scales, but the most significant one is inter-annual precipitation versus inter-annual runoff, indicating that the t runoff variation is sensitive to response to precipitation in the inter-annual time scale.For the response of runoff to temperature, there are some differences among Tavg, TNn, and TXx.Tavg has not an obvious correlation with runoff over the period 1961-2009 at any time scale.It suggests that Tavg may not be suitable for analyzing the response of runoff variation to temperature.The most significant correlation between runoff and TNn is at inter-annual and inter-decadal scale, respectively, which reveals that the the response of runoff to TNn has a long time lag.In other words, the TNn affects runoff variation slowly over the study period.For the TXx, it has statistically significant correlation with runoff at a different time except inter-annual versus inter-decadal scale, indicating that the TXx plays an important role in runoff variation over the study period in this region.Overall, to investigate the response of runoff to temperature, it is suitable for the inter-annual TNn and/or TXx (extreme temperature) versus the inter-decadal runoff over the MYZRB.As shown in Table 3, it can be noted that the precipitation has higher CCs with runoff than temperature at both inter-annual versus inter-annual scale and inter-decadal versus inter-decadal scale, however, the influence of extreme temperature (TNn and TXx) on runoff is stronger than precipitation at inter-decadal versus inter-annual scale.Overall, the precipitation and extreme temperature are both the major climate factors affecting runoff variation at different time scales over the MYZRB, moreover, the mean state of temperature (Tavg) may not be suitable for investigating the response of runoff to climate change over the MYZRB.

Discussion
Detailed analyses are performed regarding the runoff dynamic and associated multi-scale response to climate fluctuation in the MYZRB over the period 1961-2009 in this study.
The results showed that temperature had an obvious increasing trend, especially the extreme minimum temperature TNn (Figure 7), which may be caused by the large scale atmospheric circulation and human activity [31,33,34].The continuously increasing temperature resulted in a great rise in the AET [32] and thus lead to the increasing amplitude of runoff being less than precipitation (Figures 5  and 7a).A more important reason that increasing temperature reduced the runoff yield may be due to little snow and glacier covered in the study region.To confirm this assumption, we obtained the glacial area and proportion in the MYZRB were 1115 km 2 and 0.9%, respectively based on the Second Chinese Glacier Inventory (SCGI) datasets [35].Also, we computed the mean snow cover duration (SCD) in the MYZRB over the period 2001-2014, with the improved MODIS daily snow-cover dataset through cloud-removal algorithm [36].As shown in Figure 8, we can find that most regions of the MYZRB have nearly not snow cover except small part of northeastern MYZRB, and the basin-scale average value of SCD is approximately 34 days.These results confirm our assumption of little glacier and snow covered the study area.Thus, we can infer that, in this area, the stream runoff is mainly derived from precipitation, which is consistent with the literature [37].The multi-scale response of runoff to climate changes exhibited the most significant correlation between runoff and precipitation at both inter-annual versus inter-annual and inter-decadal versus inter-decadal scale, indicating that precipitation played a more important role in the overall runoff dynamic [38].For temperature, the response to runoff showed a large difference, in which the runoff has no obvious correlation with the mean state of temperature Tavg, indicating that Tavg is not suitable for investigating the response of runoff to temperature in the context of global warming [39].The extreme minimum temperature TNn has a significant correlation with runoff only at the interdecadal versus inter-annual, while the extreme maximum temperature Txx has significant correlation with runoff at a different time scale, except for the inter-annual versus inter-decadal, which indicates that TNn has a slow influence on runoff, but TXx has a direct and fast influence on runoff [31].
In this study, we do not apply any detection performance for turning point [40,41] because the trend component represents medium and long-term variations of hydro-climatic variables.We obtained the approximate turning points (e.g., early 1980s and mid-1980s) by checking RES, which should be sufficient for current study.
Although this study investigated the runoff dynamics and associated responses to precipitation and temperature (mean and extreme states), the hydro-meteorological data after 2009 is difficult to obtain due to data confidentiality, which may have an influence on the conclusion of this study.In addition, the quantitative responses of runoff to climate change and their more detailed mechanisms require further study and discussion.

Conclusions
In this study, based on the hydrological and meteorological data over the period 1961-2009 in the middle reach of the Yarlung Zangbo River Basin, we characterized the oscillation and dynamics of runoff at different time scales with EEMD, then investigated the contributions of oscillations at a different time scale to runoff dynamic and their significance level, and finally estimated the response of runoff to mean and extreme states of climate fluctuation.In addition, we illustrated the potential causes of the multi-scale responses.The main findings and physical processes involved can be summarized as follows.
The runoff of MYZRB showed an obvious nonlinear and nonstationary decreasing-increasing fluctuation, accompanying with a turning point in the mid-1980s.Its variation illustrated a weak quasi-three-year period (IMF1) and an unobvious quasi-five-year period (IMF2) at the inter-annual scale, but strong significant quasi-12-year (IMF3) and quasi-46-year (IMF4) periods at the interdecadal scale.Although the HF (IMF1-2) has weak periodic oscillation, they have larger variance contribution rates (43.2% and 16.7%, respectively) than LF (IMF3-4, 18.6% and 8.0%, respectively), The multi-scale response of runoff to climate changes exhibited the most significant correlation between runoff and precipitation at both inter-annual versus inter-annual and inter-decadal versus inter-decadal scale, indicating that precipitation played a more important role in the overall runoff dynamic [38].For temperature, the response to runoff showed a large difference, in which the runoff has no obvious correlation with the mean state of temperature Tavg, indicating that Tavg is not suitable for investigating the response of runoff to temperature in the context of global warming [39].The extreme minimum temperature TNn has a significant correlation with runoff only at the inter-decadal versus inter-annual, while the extreme maximum temperature Txx has significant correlation with runoff at a different time scale, except for the inter-annual versus inter-decadal, which indicates that TNn has a slow influence on runoff, but TXx has a direct and fast influence on runoff [31].
In this study, we do not apply any detection performance for turning point [40,41] because the trend component represents medium and long-term variations of hydro-climatic variables.We obtained the approximate turning points (e.g., early 1980s and mid-1980s) by checking RES, which should be sufficient for current study.
Although this study investigated the runoff dynamics and associated responses to precipitation and temperature (mean and extreme states), the hydro-meteorological data after 2009 is difficult to obtain due to data confidentiality, which may have an influence on the conclusion of this study.In addition, the quantitative responses of runoff to climate change and their more detailed mechanisms require further study and discussion.

Conclusions
In this study, based on the hydrological and meteorological data over the period 1961-2009 in the middle reach of the Yarlung Zangbo River Basin, we characterized the oscillation and dynamics of runoff at different time scales with EEMD, then investigated the contributions of oscillations at a different time scale to runoff dynamic and their significance level, and finally estimated the response of runoff to mean and extreme states of climate fluctuation.In addition, we illustrated the potential causes of the multi-scale responses.The main findings and physical processes involved can be summarized as follows.
The runoff of MYZRB showed an obvious nonlinear and nonstationary decreasing-increasing fluctuation, accompanying with a turning point in the mid-1980s.Its variation illustrated a weak quasi-three-year period (IMF1) and an unobvious quasi-five-year period (IMF2) at the inter-annual scale, but strong significant quasi-12-year (IMF3) and quasi-46-year (IMF4) periods at the inter-decadal scale.Although the HF (IMF1-2) has weak periodic oscillation, they have larger variance contribution rates (43.2% and 16.7%, respectively) than LF (IMF3-4, 18.6% and 8.0%, respectively), which implies that the inter-annual oscillations had been playing a more important role in the overall runoff change than the inter-decadal oscillations over the MYZRB.Also, the reconstructed inter-annual runoff series basically consistent with the fluctuation of original runoff anomaly series confirmed this conclusion.The reconstructed inter-decadal runoff series effectively revealed that the runoff series at this scale featured an alternating positive and negative phase over the period 1961-2009 in the MYZRB.According to the inter-decadal runoff series, we can infer that the runoff over MYZRB may present a reduction trend in the future.
In comparison with the variation of precipitation over the study period, the runoff variation showed a basically consistent trend (decreasing-increasing).However, they had different turning points in early 1980s and mid-1980s, respectively.Besides that, there were similar inter-annual oscillations and different inter-decadal oscillations between runoff and precipitation.As for temperature, the Tavg, TNn, and TXx presented almost continuously increasing trends and also have similar inter-annual oscillations and different inter-decadal oscillations.Higher temperature tends to rise the AET, resulting in mitigating the increasing amplitude of runoff after turning point.
Precipitation has a more significant influence on runoff dynamic than temperature at both inter-annual versus inter-annual scale and inter-decadal versus inter-decadal scale, however, at the inter-decadal versus inter-annual scale, extreme temperature was dominated, which indicated that precipitation played a more important role in runoff variation in the MYZRB over the period 1961-2009 than temperature.Moreover, the most significant scale for response of runoff to precipitation was inter-annual versus inter-annual.Additionally, the response of runoff to extreme minimum temperature (TNn) had a long time lag, namely the TNn affected runoff dynamic slowly over the period 1961-2009 in the MYZRB.
These results indicate that the EEMD is a highly effective approach to extracting IMFs and RES of the nonlinear and nonstationary hydro-climatic variables, such as runoff, precipitation, and temperature.Multi-scale periodic oscillations in IMFs and nonlinear RES confirm the conclusion that long-term runoff, precipitation, and temperature series feature nonlinear and nonstationary changing process.With IMFs and RES of hydro-climatic variables, this study could investigate runoff dynamics and associated multi-scale responses to climate changes in the MYZRB.All these findings help to better understand long-term runoff, precipitation, and temperature (mean and extreme states) variations and their internal interaction and contribute to developing better management measures to account for climate change.

Figure 1 .
Figure 1.Geographical location of the middle reach of the Yarlung Zangbo River Basin and hydrometric stations settled up in the area.The 0.5° grids denote the gridded climatic products.

Figure 1 .
Figure 1.Geographical location of the middle reach of the Yarlung Zangbo River Basin and hydrometric stations settled up in the area.The 0.5 • grids denote the gridded climatic products.

Figure 2 .
Figure 2. Variations of runoff anomaly and its five-year moving mean for the period 1961-2009 in the middle of the Yarlung Zangbo River Basin (MYZRB).

Figure 3 .
Figure 3.The IMF1-4 and trend component of runoff anomaly normalized by normal distribution over the period 1961-2009 in the MYZRB.

Figure 2 .
Figure 2. Variations of runoff anomaly and its five-year moving mean for the period 1961-2009 in the middle of the Yarlung Zangbo River Basin (MYZRB).

14 Figure 2 .
Figure 2. Variations of runoff anomaly and its five-year moving mean for the period 1961-2009 in the middle of the Yarlung Zangbo River Basin (MYZRB).

Figure 3 .
Figure 3.The IMF1-4 and trend component of runoff anomaly normalized by normal distribution over the period 1961-2009 in the MYZRB.

Figure 3 .
Figure 3.The IMF1-4 and trend component of runoff anomaly normalized by normal distribution over the period 1961-2009 in the MYZRB.

Water 2018 ,
10, x FOR PEER REVIEW 7 of 14 contribution rate was found in IMF4 for original runoff series, the most significant periodic oscillation presented similar to cosine curve with about quasi-46-year cycle was detected.The RES component variations suggested a nonlinear decreasing-increasing trend over the period 1961-2009, with the turn point in the mid-1980s.

Figure 4 .
Figure 4.A significance test for the IMF1-4 of runoff anomaly over the period 1961-2009 over the MYZRB.lnTk represents the natural logarithm of the inherent time scale (period) for each IMF, and lnEk represents the natural logarithm of the energy spectral density for each intrinsic mode function (IMF).

Figure 4 .
Figure 4.A significance test for the IMF1-4 of runoff anomaly over the period 1961-2009 over the MYZRB.lnT k represents the natural logarithm of the inherent time scale (period) for each IMF, and lnE k represents the natural logarithm of the energy spectral density for each intrinsic mode function (IMF).

Water 2018 , 14 Figure 5 .
Figure 5.A comparison among the reconstructed inter-annual and inter-decadal runoff series and the original runoff anomaly series, as well as trend component.

Figure 5 .
Figure 5.A comparison among the reconstructed inter-annual and inter-decadal runoff series and the original runoff anomaly series, as well as trend component.

Figure 8 .
Figure 8. Mean snow cover duration (SCD) derived from annual SCDs from 2001 to 2014 in the MYZRB.

Figure 8 .
Figure 8. Mean snow cover duration (SCD) derived from annual SCDs from 2001 to 2014 in the MYZRB.

Table 1 .
Statistical description for the observed annual runoff data at the national hydrological stations Lazi and Yangcun over the period 1961 to 2009.

Table 2 .
Variance contribution rates of IMF1-4 and trend component for runoff anomaly.

Table 2 .
Variance contribution rates of IMF1-4 and trend component for runoff anomaly.

Table 3 .
Correlation coefficients and significance level between runoff and climate factors over the period 1961-2009 over the MYZRB.: IA and ID denote inter-annual and inter-decadal scale, respectively; vs. is versus."***" represents significant at the level of 0.01, "**" is significant at the level of 0.05, and "*" is significant at the level of 0.1. Notes