Abstract

Watershed discharge (WD) in the alpine regions, such as the upper reach of Yarlung Zangbo River Basin (YZRB), China, could have changed severely in response to climate changes. Yet, how hydrometeorological variables varied at different time scales and how WD varied in response to hydrometeorological variables in the alpine regions remained questions to be answered. The ensemble empirical mode decomposition (EEMD) method was employed in this study to investigate the nonlinear climate change trends (averaged and extreme states) and the associated multiscale impacts on WD variations over the upper reach of the YZRB during 1961–2009. All investigated hydroclimatic variables, i.e., precipitation, temperature, and WD, were found to be varied nonlinearly with clear multiscale oscillations characterizing great differences in the oscillation periods, corresponding significance levels, and variance contribution rates, among which precipitation posed a weak impact on WD variations, while temperature played a significant role in WD fluctuations. Furthermore, among all temperature extremes, the dominant index affecting WD variations was TXm (annual mean of the daily maximum temperature) but not TXx (annual maximum of the daily maximum temperature) at both interannual and interdecadal scales, which might be caused by that TXx increased evapotranspiration and reduced WD. A significant correlation between temperature (both averaged and partial extreme states) and annual WD at both interannual and interdecadal scales indicated that a synchronous change existed between them. The present study provided first insight into how hydrometeorological variables varied at different time scales and how WD fluctuated in response to hydrometeorological variables over the upper reach of the YZRB, China.

1. Introduction

There is a broad consensus that the rapid global warming is an indisputable fact and it will increase the climate variability, frequency, and intensity of extreme events and even the redistribution of water resources [13]. Therefore, hydroclimatic change detection has been of great concern and has been playing a significant role in accurately investigating hydroclimatic change trends and exploring their complex mechanism at different time and space scales [4]. Most previous studies adopted conventional trend analysis methods such as linear regression (LR) or spline function fitting (SFF) [5, 6], wavelet analysis (WA) [7], and empirical orthogonal functions (EOFs) [8]. However, the change processes of hydroclimatic variables (i.e., precipitation, temperature, and WD) are complex, nonlinear, and nonstationary; therefore, these conventional methods may give an inaccurate or unreasonable diagnosis for the change trend of long-term hydroclimatic variables [4]. With the rapid development of signal detection technology, the ensemble empirical mode decomposition (EEMD), a new time-series signal-processing method, overcoming the problem of mode mixing, was developed based on EMD, which is more suitable for hydroclimatic trend analysis [9]. In recent years, EEMD-based hydroclimatic change detection studies have obtained certain achievements of theoretic and realistic significance. Among them, Tao et al. [2] analyzed the extreme temperatures (maximum and minimum) in Xinjiang and revealed nonlinear change trends of them. Bai et al. [4] investigated multiscale WD responses to climate change of the headwater region of Kaidu River in Xinjiang and revealed a positive correlation between precipitation and temperature at different time scales, especially at the interdecadal scale. Qian and Zhou [10] analyzed the multidecadal variations of drought in North China with EEMD and presented their temporal and spatial characteristics of their changes. However, to the best of our knowledge, little effort has been made to apply EEMD in the investigation of hydroclimatic trends in the ungauged basins in alpine regions where the discharge variations in these regions pose a great impact on water resource partitions for the downstream regions [6].

As a typical mountainous watershed with water supply mainly from the Chemayungdung Glacier in the northern slope of the Himalayas, the YZRB is well known as the highest river basin in the world, with the mean elevation exceeding 4600 m a.s.l. [11, 12]. The complex topography with high elevation of the YZRB, especially in the upper reach of the basin, has hampered the settlement of hydrometric gauges and stations [13]. Due to the potential Indian monsoon activity, the annual precipitation gradually decreased from greater than 2000 mm over the down-reach (east) to less than 300 mm on the upper reach (west) of the YZRB [14]. The stream discharge in the upper reach of the YZRB comes primarily from melting water of glaciers, snow, and frozen soils. Therefore, climate change poses a significant impact on water resource partitioning of middle streams and downstreams of this region [15]. In the context of global warming, studies of climate change impacts on WD fluctuations in the upper reach of the YZRB are of both practical and scientific significance. Therefore, the present study focused on the investigation of the nonlinear trend of climate change and its effects on discharge fluctuations at different time scales over the upper reach of the YZRB with EEMD.

The objective of this study was achieved by conducting analysis in four phases: (1) detecting the oscillations and nonlinear trend of climate change (precipitation and temperature) at different time scales; (2) characterizing the dynamics and fluctuations of annual discharge, as well as their statistical significance at different time scales; (3) quantifying the impacts of climate change on discharge fluctuations; and (4) exploring the mechanism that temperature affects annual discharge fluctuations. The remainder of this paper is structured as follows: Section 2 introduces the study area and datasets used. Section 3 describes the methodology adopted for the purpose of this study. Results and discussion as well as some uncertainty analyses are given in Section 4, followed by the conclusion being presented in Section 5 at the end.

2. Study Area and Data Utilization

2.1. Study Area

As shown in Figure 1, the upper reach of the YZRB was defined as the drainage area controlled by the Lazi hydrological station [15], geographically located between 82°00′–87°36′E and 28°46′–31°03′N. The total area of the upper reach of the YZRB is approximately 50,134 km2, accounting for 19.5 percent of the YZRB in China. The elevation ranges from 3992 to 7031 m a.s.l., with an average elevation of roughly 5090.5 m a.s.l. Being controlled by the combined effects of the Indian monsoon and westerly circulation, the upper reach of the YZRB is characterized by a semiarid climate of the cold temperate zone [5, 16] where a number of glaciers and an extensive area of seasonal frozen soils and permafrost developed in this region. According to the Second Chinese Glacier Inventory (SCGI) [17], the glacier-covered area accounts for about 1.7% of the total area in the upper reach of the YZRB. Based on the cloud-free daily snow-cover products provided by Huang et al. [18], the average snow-cover duration was found to be approximately 50.4 days. Seasonal frozen soils and permafrost predominated the study area with an areal ratio of 1.1 : 1. Due to the attenuation of the Indian monsoon activity from downstream to upstream and given the unique topography of the Yarlung Zangbo Suture Zone, the multiyear averaged precipitation in the region is below 300 mm, with the mean temperature of approximately −3.0°C [19].

The meteorological station network in the upper reach of the YZRB, organized by National Meteorological Information Center (NMIC) and China Meteorological Administration (CMA), available to public access, is only composed of a scarce and unevenly distributed set of stations [13, 19]. An insufficient number of meteorological stations could result in high uncertainty in climate change analysis. To overcome this issue, high-resolution gridded precipitation and temperature datasets (Figure 1), produced by NMIC and CMA using an intensive meteorological station network, were employed in this study [13].

2.2. Data Utilization

China daily gridded Precipitation Analysis Product and China daily gridded Temperature Analysis Product (CPAP and CTAP, respectively) are produced and routinely calibrated by NMIC and CMA, based on the observations from approximately 2472 national meteorological stations [13, 20]. All the observations used in CTAP and CPAP have undergone a rigorous quality control (QC) check in three steps of extrema detection, internal consistency check, and spatial consistency check before being interpolated [12, 21]. Both CTAP and CPAP consist of daily observations for precipitation and temperature with a 0.25° × 0.25° grid size. The CPAP has been widely used as a reference to validate the accuracy of satellite-based precipitation products [2225] and as a reference to investigate climate change [13, 26]. Also, it was pointed out that the CTAPs, including maximum (Tmax), minimum (Tmin), and mean (Tmean) temperatures, have enough accuracy for analyzing climate change at the basin scale [12, 13]. The CPAP and CTAP during 1961–2009 were selected to obtain annual total precipitation and temperature for each grid in this study. And the basin-scale annual average precipitation (Pavg) and temperature (Tavg) for all grids in the upper reach of the YZRB were calculated.

The annual WD observations at the national hydrological station of Lazi (Figure 1) during 1961–2009 were also collected in our analysis. The data have undergone a strict QC check. With these observations, we produced the annual WD anomaly data series during the same period by calculating the difference between the original WD series and its average value.

3. Methodology

3.1. EEMD

In this study, the EEMD [9, 27] was employed to extract intrinsic mode functions (hereafter referred to as IMFs) and the nonlinear trend component (hereafter referred to as RES) from the nonlinear and nonstationary hydroclimatic data series during 1961–2009. The IMFs revealed the oscillation characteristics from high frequency (HF, less than a 10-year period) to low frequency (LF, not less than a 10-year period) at different time scales, while RES represented the nonlinear trend in the original time series [12]. This method consisted of self-adaptation mechanisms and the temporal locality, which was suitable for time-frequency analysis of nonlinear and nonstationary variables with long-term variations [2]. For a given time series in EEMD, x(t) can be expressed as follows:where is the number of IMFs. The jth IMF can be defined in polar coordinates as follows:where is the time-dependent amplitude; represents the time-dependent frequency; and denotes the trend component, namely, the RES. IMFs are different from Fourier modes, for which both and are time-independent. The following two properties define an IMF: (1) each IMF has exactly one zero crossing between two consecutive local extremes, namely, a sequence of maxima and minima, and (2) the local mean of each IMF is zero [28]. In this study, the significance of the IMF components was tested by a posteriori test method reported in [9, 27]. First, the distribution of the energy concerning the period in the form of the spectral function should be obtained. The IMF’s energy density () can be expressed using the following equation:where is the length of an . Next, the mean period of the () can be defined as follows:where is the number of peaks for each IMF. There is a simple equation that relates the average energy density and the mean oscillation period as follows [9]:

As shown in Figure 3, the relation between (the x-axis) and (the y-axis) can be expressed by a straight line with the slope of −1 [4]. Theoretically, the IMF component of the white noise series is supposed to be distributed on the line; however, due to little deviation within the actual application, the confidence interval for the energy spectrum distribution of white noise is presented as follows [29]:

Through decomposition, if the energy of an IMF component exceeds the upper limit of a given confidence interval (i.e. α = 0.05, 0.10, 0.80, and 0.50), it indicates that the periodic oscillation is significant at the corresponding significance level [4, 30]. Furthermore, the higher significance level means the more actual physical meaning is contained in the IMF component; otherwise, the more white noise is contained in that component [29].

One of the fundamental approaches to overcome the mode mixing issue in the EEMD is to add a white Gaussian noise series, with a probability density function (PDF) equal to that of the normal distribution, into original series x(t). It is noteworthy that the scaling factor, controlling amplitude of noise, should be neither too small nor too large; here, its optimum value ranges from 0.2 to 0.3, and the ensemble number is set to 100 [4]. In this study, we conducted the decomposition of hydroclimatic data series with the scaling factor ranging from 0.2 to 0.3 but found little differences in IMFs and RES. Therefore, the white Gaussian noise series with the amplitude of 0.2 times the standard deviation was considered in EEMD we employed.

3.2. Calculation of Temperature Extremes

To explore how the temperature extremes affect WD in the context of global warming, we computed four temperature extremes (TNn, TNm, TXm, and TXx) from CTAP datasets in the RClimDex software (v1.1; http://etccdi.pacificclimate.org/software.shtml). Another two temperature extremes (TMn and TMx) were calculated using the IDL codes developed in the present study. The details of temperature extremes used in this study, including ID, name, definitions, and units, are listed in Table 1.

4. Results and Discussion

4.1. Analysis of Climate Change

The average annual precipitation (Pavg) and temperature (Tavg) over the upper reach of the YZRB during 1961–2009 were, respectively, 299.73 mm and −3.03°C. Anomaly series of Pavg and Tavg are illustrated in Figure 2. Increasing and decreasing alternative trends in annual precipitation variations were found during the period, accompanying with a distinct turning point in the mid-1980s (Figure 2(a)). The precipitation was relatively higher in 1970–1985 than other periods and was followed by several fluctuations over 1986–2009.

As shown in Figure 2(b), the temperature in the upper reach of the YZRB tended to be increasing all over the study period, with the turning point in the early 1990s. Before this time, the temperature rose slowly with stronger fluctuations but with a stable, increasing trend with small fluctuations afterward. The average temperatures for 1961–1990 and 1991–2009 were −3.42°C and −2.41°C, respectively, with the maximum (occurred in 2006) to minimum (occurred in 1966) mean annual temperature ratio of 3.55 : 1. Referring to the 5-year moving average of Pavg and Tavg anomaly series during 1961–2009 (Figure 2), nonlinear and nonstationary variations in precipitation and temperature can be recognized clearly. Therefore, the EEMD, as a nonlinear method of extracting trends and period information from a given time series, was suitable to be applied to decompose Pavg and Tavg over the study period.

After the decomposition, we obtained four IMFs (IMF1–4) and one trend component (RES) for each time series (Pavg and Tavg). These IMFs revealed oscillation characteristics of Pavg and Tavg from HF to LF at different time scales, and the RES reflected the trend of the original time series over time. The periodic significances of IMFs for Pavg and Tavg are illustrated in Figures 3(a) and 3(b), respectively, in which ln Tk and ln Ek represent the natural logarithm of the original period and the energy spectral density for each IMF, correspondingly. According to Figure 3(a), the oscillation periods of IMF1 and IMF3 were not able to reach the 90% confidence limits, indicating that precipitation over the upper reach of the YZRB had a weak quasi-3-year period and a weak quasi-10-year period. However, IMF2 and IMF4 fall above 95% confidence limits, implying that the precipitation in the studied region has a significant quasi-5-year period and a significant quasi-39-year period. For the periodic significance of Tavg (Figure 3(b)), a weak quasi-3-year period appeared in IMF1 and a weak quasi-42-year period in IMF4, while two apparent quasi-5-year and quasi-12-year periods appeared in IMF2 and IMF3, respectively. Similar to temperature, precipitation had stronger HF periodic oscillations (IMF1-2) during 1961–2009 than LF periodic oscillations (IMF3-4). These multiscale oscillations of precipitation and temperature reflected not only the periodic variations of the climate system under external forcing but also the nonlinear feedback of the climate system [30].

The detailed quasiperiods and their corresponding variance contribution rate of each IMF and RES to Pavg and Tavg series are listed in Tables 2 and 3, respectively. According to the results, Pavg and Tavg had similar oscillation periods in IMFs, especially IMF1–3, with some differences in significance levels (Figure 3). The quasi-3-year oscillation period of Pavg and Tavg might be explained by the tropospheric biennial oscillation (TBO) with a 2- to 3-year period, which is the basic characteristic of interannual variations of atmospheric circulation [31]. In this study, a weak TBO of precipitation and temperature variations was identified over the upper reach of the YZRB. The quasi-5-year oscillation period of Pavg and Tavg may be related to the El Niño-Southern Oscillation (ENSO) and North Atlantic Oscillation (NAO) [32, 33]. It can be concluded that the ENSO and NAO pose a significant impact on precipitation variation but had an insignificant impact on temperature variation over the upper reach of the YZRB. The quasi-11-year (Pavg) and quasi-12-year (Tavg) oscillation periods were probably dependent on solar activity (e.g., roughly quasi-11-year period of the sunspot numbers) [34]. Detailed explanation of the mechanism on how multiscale oscillations of precipitation and temperature formed requires further study of the large-scale oscillations in the future. Relating variance contribution rates of each IMF and RES to precipitation and temperature series, the HF oscillations (IMF1-2) were stronger than LF oscillations (IMF3-4), especially for the temperature. However, LF oscillations cannot be ignored to maintain the total energy of the original time series. Also, HF oscillations had more noises than LF oscillations, since all IMFs were not statistically significant for temperature.

According to the theory of EEMD [9], the interannual precipitation and temperature series were reconstructed by summing IMF1, IMF2, and RES, while their interdecadal series were obtained by summing IMF3, IMF4, and RES. Interannual and interdecadal fluctuations as well as nonlinear change trends of Pavg and Tavg anomaly series are exhibited in Figures 4(a) and 4(b), respectively. From Figure 4(a), we can see that the oscillations of the interannual precipitation series were nearly consistent with the oscillations of the original precipitation anomaly series during the 1961–2009 period except from 1984 to 1995, which might be attributed to the external noises contained in the HF oscillations [4], especially IMF1. The interdecadal precipitation series had positive- and negative-phase alterations during 1961–2009, representing wet and dry periods over the study period. Moreover, it can be inferred that the precipitation tends to be decreasing in general in the future. The RES presented an increasing-decreasing change trend during 1961–2009, with a turning point in the mid-1980s. As for temperature (Figure 4(b)), the reconstructed interannual and interdecadal variation series perfectly captured the small-scale and large-scale oscillations of the original anomaly series. According to the interdecadal variation series, the temperature tends to be continually increasing in the future. Also, the RES indicated a temperature increase over 1961–2009, with an accelerated change rate after the early 1990s.

4.2. Multiscale Impacts of Climate Change on WD Fluctuations

In the context of global warming, climate change usually results in changes in hydrological cycles and causes repartitioning of water resources in time and space, especially in the mountain areas (Deng et al. [6] and Bai et al. [4]). To analyze the impacts of precipitation and temperature on water resources, four IMFs and one RES of WD series were obtained using EEMD in the upper reach of the YZRB. As the significance test results of the four IMFs presented in Figure 5, WD had relatively significant HF oscillation periods than LF oscillation periods. However, neither HF nor LF oscillation periods were statistically significant. Furthermore, the WD had quasi-3-year and quasi-6-year periods falling into confidence limits between 80% and 90%, while relatively weaker quasi-11-year and quasi-28-year periods falling into confidence limits between 50% and 80% were also identified. Therefore, WD series had similar oscillation periods to precipitation and temperature. To investigate the information contained in each IMF and RES for original WD series, we computed their variance contribution rates as presented in Table 4. The variance contribution rates of HF components (IMF1-2) were higher than those of LF components (IMF3-4), with the RES component of the maximum variance contribution rate about 39.2%.

To investigate the WD variations at different time scales, similar to precipitation and temperature, we generated the interannual (the sum of IMF1, IMF2, and RES) and interdecadal (the sum of IMF3, IMF4, and RES) WD series during 1961–2009, and the results are presented in Figure 6 against the original WD anomaly series as well as RES. The reconstructed interannual WD series portrayed the variations of the original WD anomaly series over the study period, implying that the interannual variations play an essential role in explaining the differences in the initial WD series over the upper reach of the YZRB. However, the interannual change series was not highly consistent with the original WD series since the early 1990s to a greater extent because of the systematic external noises contained in the HF oscillations [4]. The reconstructed interdecadal WD series depicted the alternating abundance and shortage stages over the study period. The interdecadal series, which is obtained by superimposing large-scale oscillations on a decreasing-increasing hydrological trend, was a secular and complicated process of WD variations. Therefore, according to the fluctuations of the interdecadal WD series, it can be inferred that WD might be decreasing in the future. In comparison with the original WD anomaly series, the RES can depict perfectly the decreasing-increasing change trend, with the turning point in the early 1990s.

To quantify the impacts of climate change on WD variations, we calculated the multiscale correlation coefficients between precipitation, temperature, and annual WD at different times (Table 5) and the corresponding significance levels. The precipitation and WD have a nonsignificant correlation at all different time scales, which implied that precipitation had a slight influence on WD variations probably due to the semiarid climate with less precipitation but mainly in the solid state of the area. As for temperature, significant correlations with WD at the IA vs. IA, IA vs. ID, and ID vs. ID scales were found, especially at the ID vs. ID scale with the significance level of 0.01 (Table 5). The correlation between temperature and annual WD at the large scale (ID vs. ID) was stronger than that at the small scale (IA vs. IA), with the correlation coefficients of 0.39 and 0.35, respectively. Therefore, we can conclude that annual WD varied much more sensitively and quickly in response to temperature changes than to precipitation changes over the semiarid, cold upper reach of the YZRB, China.

Studies of multiscale impacts of climate change on WD variations revealed a certain level of understanding on hydrological responses to climate changes over the semiarid, cold upper reach of the YZRB; however, it has been pointed out that climate extremes have more direct and significant effects on the hydrologic regime than normal states [12, 35]. Given the insignificant impact of precipitation on WD over the upper reach of the YZRB, the effect of temperature extremes on WD was thus investigated further in this study. Similar to procedures as applied to normal states of climate changes and WD series, temperature extremes including TNn, TNm, TMn, TMx, TXm, and TXx were decomposed by using EEMD, and interannual and interdecadal series of these temperature extremes were also reconstructed in the similar way as did for climate change indices. Finally, the correlation coefficients between temperature extremes and WD at two primary time scales (i.e., IA vs. IA and ID vs. ID) were obtained by the multiscale correlation analysis method (Table 6). Accordingly, the high-temperature indices (TMx, TXm, and TXx) had a relatively greater impact on WD than low-temperature indices (i.e., TNn, TNm, and TMn) at both interannual and interdecadal scales, and even the correction was not statistically significant for low-temperature indices. This implied that WD variations can be better explained by high-temperature indices over the upper reach of the YZRB, which coincided with the common sense that the high temperature extremes determine the hydrological regime of the upper reach of the YZRB characterizing less rainfall but abundant snow and glaciers [5, 36]. Therefore, melting water from snowing, snow cover, glaciers, and frozen soils predominated WD of this region [15, 37], and high temperatures played an important role in alpine basins (Deng et al. [6] and Liu et al. [19]). At the interannual scale, TXm most closely correlated with annual WD with a correlation coefficient of 0.54, followed by TMx and then TXx. The TXx, as a proxy for high temperature, can result in an increase of evapotranspiration and a decrease of WD [6, 12]. Therefore, the TXm was considered as the most representative temperature extreme for investigating the impacts of temperature on WD variations at both interannual and interdecadal scales in the upper reach of the YZRB.

4.3. Discussion

In this study, based on the secular hydroclimatic datasets observed during 1961–2009 in the upper reach of the YZRB, oscillations and nonlinear trend of climate (precipitation and temperature) changes and the WD variations at different time scales were detected by using EEMD. All investigated hydroclimatic variables, i.e., precipitation, temperature, and WD, were found to be varied nonlinearly with clear multiscale oscillations characterizing great differences in the oscillation periods, corresponding significance levels, and variance contribution rates. Among them, precipitation posed a weak impact on WD variations, while temperature played a significant role in WD fluctuations. At the interannual scale, TXm highly correlated with annual WD, followed by TMx and TXx. The TXx, with increasing evapotranspiration and reducing WD, did not show a significant correlation with WD at both interannual and interdecadal scales. Significant correlations between partial temperature indices (i.e., Tavg, TXm, and TXx) and annual WD were detected at both interannual and interdecadal scales, which indicated that a synchronous change existed between them. Based on above-described studies, interannual and interdecadal variations as well as a nonlinear change trend of hydroclimatic variables were constructed. According to the reconstructed interdecadal variation series, the temperature tends to be continually increasing in the foreseeable future. Also, the RES indicated an increasing tread of temperature over 1961–2009, with an accelerated change rate after the early 1990s. However, the reconstructed interdecadal series of WD were rather secular and complicated, and the WD might be increasing in a certain period of time near the future but tends to be decreasing in the foreseeable future. Uncertainties do exist in the prediction of hydroclimatic variations in the results of the present study. One source of uncertainties introduced in the analysis of hydroclimatic variations of the study area is due to inaccessibility of the hydroclimatic data to public after 2009 because of the data management policies, which might bring negative impacts on the prediction of future changes in precipitation, temperature, and annual WD. Another one came from the limited hydroclimatic data available in the upper reach of the YZRB, which resulted in that all the analyses were constrained on the impact of precipitation and temperature on WD variations without taking spatial heterogeneity of those variables at the basin or subbasin scales into account. However, according to Guo et al. [32], there might be notable spatial differences in multiscale oscillations and change trends over the study period in the study area. Although some reasons were provided for why TXx cannot be regarded as the most robust temperature index to account for WD variations, quantification on the impacts of evapotranspiration on WD variations was not attempted due to the lack of accurate observation data.

The extensively developed seasonal frozen soil and permafrost contribute a large portion of WD in the upper reach of the YZRB. In the context of global warming, both seasonal frozen soils and permafrost have undergone a noticeable degradation, leading to a significant increase in thickness of the active layer [38, 39]. Conversely, the increasing active layer thickness of frozen soil posed a great impact on WD variations and meanwhile involved the associated unknown responses to climate change, such as higher evaporation, greater water storage in the soil profile, and less surface runoff yield. Although climate change has been taken as the main driving force for the changes of the hydrologic system, the effects of frozen soil degradation on WD variations yet remained as an important issue poorly understood [40]. Therefore, further studies are required to explore the quantitate effects of frozen soil degradation on WD variations in the upper reach of the YZRB, China.

Attempts have been made to detect the turning points for precipitation, temperature, and annual WD series accurately with the Mann–Kendall test [32]. However, no significant mutations were found, especially in precipitation and annual WD series, given that their secular variations are nonlinear and nonstationary over the upper reach of the YZRB. To investigate the temporal effect of climate changes on annual WD variations, the approximate turning points were identified (e.g., mid-1980s and early 1990s) using RES [4, 12], which was believed to be sufficient for the present study with the analysis of secular hydroclimatic variations.

5. Conclusions

This study investigated nonlinear and nonstationary climate change (averaged and extreme states) and associated impacts on WD variations using EEMD over the period 1961–2009 in the upper reach of YZRB. The precipitation, temperature, and annual WD presented obvious and similar multiscale oscillations, such as HF (quasi-3-year and quasi-5/6-year periods) and LF (quasi-11/12-year period) oscillations; however, there were some differences in the longest oscillation period of them; the longest oscillation period of annual WD (quasi-28-year period) was smaller than that of precipitation (quasi-39-year period) and temperature (quasi-42-year period). Moreover, the significance levels and variance contribution rates of these multiscale oscillations for precipitation, temperature, and annual WD also had great differences. Overall, these multiscale oscillations and trend component indicated that the variations of precipitation, temperature, and annual WD are nonlinear and nonstationary in the upper reach of YZRB, and the EEMD is a sophisticated change detection method for estimating hydroclimatic trends at the basin scale.

With respect to the impacts of climate changes on WD variations, precipitation posed a weak impact on WD variations at all time scales (i.e., IA vs. IA, IA vs. ID, ID vs. IA, and ID vs. ID), while temperature played a significant role in WD variations in the upper reach of YZRB at interannual and interdecadal scales, with the correlation coefficients of 0.35 and 0.39, respectively. Furthermore, the high-temperature indices (TMx, TXm, and TXx) had stronger impacts on annual WD than low-temperature indices (TNn, TNm, and TMn) at both interannual and interdecadal scales. TXm was the most significant extreme temperature index among all temperature extremes, while TXx was not a dominant temperature index for reducing WD. Moreover, the Tavg, TXm, and TXx had significant correlations with annual WD at both interannual and interdecadal scales, indicating that a synchronous change existed between them over the upper reach of YZRB.

Overall, uncertainties do likely exist in this study. However, these results offer an opportunity to better understand multiscale oscillations and change trend of secular hydroclimatic series and the impacts of climate change on annual WD at different time scales in the upper reach of YZRB, which will help us develop better water resource management to account for climate change impact.

Data Availability

The meterological data used to support the findings of this study are available from the corresponding author upon request. The watershed discharge data used to support the findings of this study have not been made available because of the data management policy.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

The authors contributed equally to this paper.

Acknowledgments

This study was supported by the National Key R&D Program of China under grants 2016YFA0602302 and 2016YFB0502502. The authors would like to thank all data providers and Dr. Masoud Jafari Shalamzari from the Institute of Remote Sensing and Digital Earth, CAS, for polishing the language of this paper.