Spatiotemporal Variation of Drought and Associated Multi-Scale Response to Climate Change over the Yarlung Zangbo River Basin of Qinghai–Tibet Plateau, China

: Drought is one of the most widespread and threatening natural disasters in the world, which has terrible impacts on agricultural irrigation and production, ecological environment, and socioeconomic development. As a critical ecologically fragile area located in southwest China, the Yarlung Zangbo River (YZR) basin is sensitive and vulnerable to climate change and human activities. Hence, this study focused on the YZR basin and attempted to investigate the spatiotemporal


Introduction
Drought is one of the most threatening natural disasters in the world, which has terrible impacts on agricultural irrigation and production, ecological environment, and socioeconomic development [1], and it is generally categorized into four types: meteorological, agricultural, hydrological, and socioeconomic [2].Meteorological drought (represented as drought in this study) is a prerequisite to identify other drought categories, which is characterized as the occurrence of below-average precipitation over a long period of time and across a wide area [3,4].Under the background of global warming, the frequency and intensity of drought have increased significantly [5].Therefore, it is of great importance for drought monitoring and warning to understand the spatiotemporal variation of drought and associated multi-scale response to climate change.
Influenced by multiple climate factors, drought often displays nonlinear, nonstationary complex processes with periodic oscillations [6][7][8].However, previous research mostly focused on linear variations and neglected the nonlinear responses to climate change [9][10][11].Hence, investigating the nonlinear and nonstationary process is essential to understand long-term drought variations.As a novel signal processing method, the ensemble empirical mode decomposition (EEMD), having advantages of robust self-adaptability and local variation characteristics, was developed to solve nonlinear and nonstationary problems [12,13].The intrinsic mode functions (IMFs) and the trend component of climate factors are separated from the original signal based on the EEMD method [14,15].Hence, the EEMD method is one of the most popular techniques to extract periodic information and changing trends, which has been applied to analyze long-term nonlinear variations in climate research [16,17].
The scPDSI (self-calibrating Palmer drought severity index), based on evapotranspiration estimation using the FAO-56 Penman-Monteith equation, has been considered as the best drought index in China [18][19][20], especially for monitoring and warning drought events.Firstly, it has a solid physical foundation by considering precipitation, temperature, and the locally-available water content to assess the precipitation deficit/surplus [21].Whereas, the SPI (standardized precipitation index) only focuses on the precipitation [22], and the SPEI (standardized precipitation evapotranspiration index) ignores soil moisture [23].Secondly, compared with the empirical Thornthwaite (Th) equation that only describes the mean daily temperature and latitude, the FAO-56 Penman-Monteith equation incorporates the energy balance and aerodynamic theory [24].Thirdly, it is more suitable for global wide drought monitoring as it automatically calibrates itself at any location using dynamically computed values based on the framework of the PDSI model [25].In addition, the seasonal snowpack dynamic variation is considered in the water balance model, which provide a more accurate measure of the availability of moisture for snowy areas when the snowpack melts [26].This feature is particularly suitable for high altitude regions covered with glaciers and snow, such as the Yarlung Zangbo River basin in the Qinghai-Tibet Plateau.Finally, the high accuracy scPDSI dataset produced by the CRU (climate research unit) has been successfully applied in different regions over the world [27][28][29].
Long-term variation of drought is closely associated with climate variables, such as precipitation, temperature, potential evapotranspiration, etc.In general, less precipitation and higher temperature will lead to more soil moisture consumed by evapotranspiration and furthermore result in drought events.Although comprehensive research concerning the response mechanism of drought to climate variables has been extensively carried out [30][31][32], it is still unclear what the response mechanism between extreme temperatures (i.e., maximum, minimum temperature and diurnal temperature range) and drought variation process, especially in high-altitude regions featured with the fragile ecological environment.In addition, the multi-scale response mechanism of drought under global climate change should be considered more.
Since the twentieth century, severe drought has brought great economic and social losses over China, and the dry area has increased by 3.72% per decade throughout China [33].Chen and Sun (2015) further detected drought variation characteristics and considered that drought has been frequent and intensive across China [34].Li et al. (2009) pointed out that significant decreasing soil moisture may induce severe drought since the 1950s [35].The Qinghai-Tibet Plateau, located in southwest China, also known as the Third Pole of the Earth and the Water Tower of Asia, is more sensitive to global climate change [36].Currently, it is one of the ideal and perfect regions to study global climate change [37].Recent investigations have indicated that drought in the Qinghai-Tibet Plateau not only exerted important influence to local agricultural production and social development [38,39], but also gave feedback to global climate because it is an elevated heat source and sink that drive global circulation [40].The Yarlung Zangbo River (YZR) basin (82 ) is a key region for the precipitation generation mechanism over the Qinghai-Tibet Plateau.In addition, the YZR contains enormous water vapor transportation channels which carry abundant moisture from the Indian Ocean to the inner region of the Qinghai-Tibet Plateau [41].Thus, the YZR basin is of overwhelming importance to agricultural development, sustainable environment, and economic prosperity of the Qinghai-Tibet Plateau [42].Furthermore, the YZR basin is also the political, economic, and cultural center of Tibet Autonomous Region, China.The cultivated land, mainly distributed in the river valley of the midstream, accounts for about 62.89% of the area of the Tibet Autonomous Region.However, the statistically-significant warming and intensive drought were observed in the YZR basin [43].Li et al. (2010) and Song et al. (2011) concluded that the mean temperature showed an evident increasing trend [44,45], which would result in severe drought in the YZR basin.Li et al. (2015) demonstrated that under the context of global warming, the duration and magnitude of drought had gradually worsened [46].Desertification in the YZR basin has been a critical issue during the past few decades due to the significant warming [47].The available water resources for agricultural production, economic prosperity, and society development, significantly affected by warming and drought, will become more challenging in the future.Previous investigations more focused on the linear variation of drought, while non-linear and non-stationary drought process and associated multi-scale response to climate change are still indistinct.In addition, inadequate and irregular observation networks in the YZR basin have been restricting the ability to study basin-wide climate change.Therefore, to solve the problems mentioned above, the primary objectives of this study are as follows: (1) To explore the spatiotemporal variation of drought in the YZR basin based on scPDSI; (2) to identify periodic cycles at different time scales and reconstruct time series of the scPDSI based on the EEMD method; and (3) to investigative the multi-scale response mechanism of drought under global climate change.The results of this study will contribute to our understanding of drought in the YZR basin and provide useful information for ecological environment protection and socioeconomic development.

Study Area
The YZR basin, one of the largest river basins in southwest China, is sited in the southern edge of the Qinghai-Tibet Plateau that is located at a longitude of 82 • 0 E-97 • 1 E and latitude of 27 • 8 N-31 • 2 N (Figure 1).The YZR flows through the south of Qinghai-Tibet Plateau from west to east, with a total length of ~2000 km and an entire area of ~240,000 km 2 in China.Originating from the northern foothills of the Himalayas Glacial on the Qinghai-Tibet Plateau, it serves as the most important freshwater source for downstream areas.In this study, the YZR basin was divided into three sub-areas with the ArcSWAT tool (upstream, midstream, and downstream).The climate in the YZR basin, strongly influenced by the southwest monsoon, is complicated.65%-80% of the annual precipitation amount falls in June to September, while precipitation mostly comprised of snow in winter and early spring is relatively little, indicating the precipitation has an obvious uneven intra-annual variation [48].Precipitation in this region is characterized by a large spatial heterogeneity with a decreasing trend from the downstream to upstream areas.The annual mean precipitation of the cold and dry upstream is less than 300 mm, the annual mean precipitation in the midstream is 300-600 mm, whereas the annual precipitation amount of the warm and humid downstream is more than 2000 mm.The YZR basin is covered with extensive glaciers and snow cover.The area of glaciers reaches 4225 km 2 accounting for about 2.1% of the YZR basin and nearly 30% of the area is covered by perennial snow and seasonal snow [49,50].Furthermore, the YZR basin has shown a remarkable increase in temperature over the past decades, which is consistent with the warming trend across the whole Tibetan Plateau [51][52][53], and increasing runoff of the YZR basin was noted in the recent years, which can be attributed to the melting of the glacier and snow coverage over the basin [54].

Data
Previous research has typically considered in-situ observations to analyze the spatiotemporal variation of climate factors [55][56][57].However, the observational networks are irregularly distributed, which reduces the representativeness of the spatial variability.In addition, most stations lack long-term and continuous observation data [58,59].The YZR basin is a typical high-altitude region with an inadequate and irregularly observation network [60], which has only 15 national meteorological stations mainly concentrated in the midstream and downstream, while its drainage area is over 240,000 km 2 .
To overcome the problems that exist in sparse observational networks, the CRU data have been used extensively in recent climate research to analyze the historical climate factors of the world [61][62][63][64].The CRU data was developed and has been subsequently updated, improved, and maintained with support from a number of funders, principally the UK's Natural Environment Research Council (NERC) and the US Department of Energy.Long-term support is currently provided by the UK National Centre for Atmospheric Science (NCAS), a NERC collaborative center.Based on surface meteorological stations, monthly values of climate variables firstly were constructed for the reference period by applying quality control techniques, interpolation methodologies, and cross-validation.Then, the time series of climate variables were filled by interpolating the observed data.In addition, sometimes one variable was used to estimate another: for example, sunshine duration and daily temperature range were used to estimate cloudiness.Detailed information on the procedures used to build this database are given in New et al. [65,66].The CRU data has a spatial and temporal resolution of 0.5 • × 0.5 • and one month, respectively.The CRU data was validated and applied across the Qinghai-Tibet Plateau, indicating a reasonable representation of spatiotemporal variations of precipitation and surface air temperature [67,68].
In this study, precipitation (PRE), mean temperature (TEM), potential evapotranspiration (PET), maximum (TMX) and minimum temperature (TMN), and diurnal temperature range (DTR) were extracted from the CRU TS 3.26 data.The PET data were calculated based on the FAO-56 Penman-Monteith equation as recommended by the World's Food and Agriculture Organization (FAO).
The scPDSI data were produced based on CRU TS 3.26.It was calculated using the time series of precipitation, temperature, and fixed parameters related to the soil/surface characteristics at each location, which made this scPDSI data more appropriate for drought monitoring and warning [69].
The range of values of the scPDSI was between −4 and +4.Negative scPDSI values indicated dry condition; positive scPDSI values indicated wet condition.The classification of the PDSI values is summarized in the below charts (Table 1).

Linear Regression Method
The linear regression method was applied to detect the changing rate of scPDSI and climate variables in the YZR basin from 1956 to 2015.The regression model was as follows: where y is the scPDSI and climate variables; t is the year; a is the slope, indicating the increasing rate when a > 0 or decreasing rate when a < 0. Specially, when variable y is scPDSI, a > 0 indicates the wetting tendency, a < 0 indicates the drying tendency; b is the intercept of the regression.
Locally-weighted scatter point smoothing (LOWESS), a nonparametric method, has been also used to smooth the time series and visualize changing trends in climate factors [70][71][72].With the LOWESS, each smoothed value was determined by neighboring data points within the span and locally-weighted polynomial regression [73].

Mann-Kendall Significance Test
As an effective and practical statistical method, the nonparametric Mann-Kendall significance test, widely applied in the research field of climate change, mainly aims to identify trend significance of discharge, precipitation, temperature, and other climate factors [74][75][76].The Mann-Kendall test statistic Z c was estimated as follows: To a series of X t = (x 1 , x 2 , . . ., x n ), where with where x i and x j are the sequential data values; n is the length of the time series; t is the extent of any given time.The null hypothesis H 0 is rejected if , where α is the significance level of the test.
In this study, the trend was significant when |Z c | > 1.96 corresponding to the significance level of 0.05.In addition, Sen's slope coupled with Mann-Kendall test was used to carry out trend analysis in this study.β, the Sen's slope, was calculated as follows:

Moving t-test
The principle of the moving t-test to identify the abrupt point is to test whether the mean values of two sub-samples change significantly at the significance level of 0.05 [77,78].Two sub-datasets are manually separated by setting a datum point, but generally the lengths of two subsamples are equivalent.The statistical value of the moving t-test method was calculated as follows: where x 1 , x 2 are the two sets of subsamples; n 1 , n 2 are the lengths of the subsamples; x 1 , x 2 and s 1 , s 2 are mean values and standard deviations of the two sub-datasets, respectively.When t > t 1−α/2 (i.e., α was the significance level of the student test and equals 0.05 in this study), it can be considered that an abrupt change takes place at the significance level of 0.05.

Ensemble Empirical Mode Decomposition
The ensemble empirical mode decomposition (EEMD) with robust self-adaptability and local variation characteristics, can be applied to solve nonlinear and nonstationary problems.Therefore, we applied EEMD method to decompose the scPDSI and climate factors from 1956 to 2015 into different intrinsic mode functions (IMF i ) and a trend component (T).Each IMF component, consisting of a signal and white noise of finite amplitude, was defined as the mean of an ensemble of trials [79,80].The natural waveforms were applied in the EEMD method compared with other decomposition methods, and the physical information with specific time scales did not change when the new data was added [81,82].The original time series x(t) was decomposed according to the formula below.

Significance Test of IMF Components
To determine the oscillation scale of each IMF component, we examined a more detailed distribution of the energy for the period in the form of a spectral function.The energy density of IMF components (E i ) can be defined as follows: where N is the length of the IMF component.The white noise sequence was tested using the Monte Carlo method.A simple equation that related the energy density (E i ) and the averaged period (T i ) was obtained: If we plot ln(T i ) as the X-axis and ln(E i ) as the Y-axis, the relation between the energy density and the averaged period can be expressed by a straight line whose slope is −1.
In theory, the IMF components of the white noise series should be distributed on the line; however, the actual application produced a little deviation, and the confidence interval for the energy spectrum distribution of white noise was; thus, presented as follows: where α is the significance level.

Variation Contribution Rate
The variance contribution rate (W) is a measure of the effects of the frequency of the fluctuation and amplitude at different scales on the original signal characteristics [83,84].The variance contribution rate of the IMF i can be calculated as follows: where var(IMF i ) and var(T) are, respectively, the variances of the IMF i and the trend component.W i is the variation contribution rate of the IMF i .

Temporal Variation Analysis
By averaging scPDSI values over all pixels, the linear regression and LOWESS method were applied to identify the long-term variation trend of scPDSI anomaly from 1956 to 2015 in the YZR basin.The results are shown in Figure 2. The scPDSI anomaly showed a significant increasing trend with the slope value of 0.024/yr (R = 0.347, Z c = 2.54), indicating the YZR basin was, overall, becoming wetter during the study period.Yuan et al. (2017) used the ration of potential evapotranspiration to precipitation to detect spatial and temporal variations of the wet-dry condition over China, and concluded that the degrees of wetness in the Qinghai-Tibet Plateau had substantially increased from 1961 to 2015 [85].Owing to the better capability of the LOWESS method to capture the intrinsic variation trend by smoothing time series, an obvious climate transition period from wet to dry during the mid 1990s was detected (brown line in Figure 2).Especially, the scPDSI showed a significantly decreasing trend after the transition period, indicating that the YZR basin turned drier.Zhang et al. (2019) came to the same conclusions, showing that soil moisture kept a persistently decreasing trend in the Qinghai-Tibet Plateau since the twenty-first century [86], which indicates a drier condition.Though the LOWESS curve shows the value is dropping, the end of 2015 was still wetter than the beginning of 1956, which restated that the YZR basin showed a basin-wide wetting tendency.In conclusion, although the YZR basin showed an overall wetting tendency, it experienced a "wetting-to-drying transition" process during the mid 1990s, as shown by the results of the LOWESS method.To further investigate and validate the transition characteristics of the dry-wet condition inferred by the LOWESS curve, the moving t-test method was utilized to identify the occurrence time of the breakpoint.Figure 3 shows the moving t-test results of the annual scPDSI anomaly.There were two breakpoints located over the 95% confidence interval.One located in 1975 and the other located in 1994.The result of the transition point in 1975 was close with the opinion that the occurrence frequency of drought showed an obvious change around the 1980s in China [87].The transition point of 1994 was consistent with results detected by the LOWESS curve, illustrating a transition period from wet to dry occurred in the middle of the 1990s.

Periodic Cycle Analysis
EEMD method was utilized to decompose the scPDSI anomaly from 1956 to 2015 in the YZR basin, and; furthermore, to identify the periodic cycles by using the significance test.As shown in Figure 4, four IMFs and one trend component were decomposed, which contained the periodic changes and nonlinear feedback in the climate system [88].According to Section 2.3.5, the significance test results of different IMF components are plotted in Figure 5, where the energy density values and inherent periodic cycles are shown with the logarithmic form at the longitudinal axis and the horizontal axis, respectively.The periodic cycles of IMF components are listed in Table 2.The IMFs, with their specific physical meaning, reflected particular inherent oscillation in the original scPDSI anomaly series [89].Figure 5 and Table 2 indicate that the scPDSI anomaly had the quasi-3-year (IMF1), quasi-9-year (IMF2), quasi-15-year (IMF3), and quasi-56-year (IMF4) cycles.In addition, the IMFs decomposed by the original scPDSI anomaly all fell above the 90% significance level line, indicating each IMF component had an important inherent physical meaning (Figure 5 and Table 2).The trend component showed a nonlinear and nonstationary characteristic during 1956-2015, which exhibited an upward trend before the middle 1990s and a downward trend after 1990s, which is similar to the result in Section 3.1.Combined with the result of LOWESS and the trend component decomposed by EEMD, although the YZR basin experienced an overall wetting process from 1956 to 2015, a transition period from wet to dry occurred in the mid 1990s, implying a significant drying period after 2000, which should be taken into consideration seriously for drought monitoring and warning.According to the Tibet Natural Disaster Statistics from the Third Pole Environment Database (http://www.tpedatabase.cn),Nyingchi and Shigatse counties, located in the YZR basin, indeed faced continuous drought in 2001 and 2002, which brought great challenges to the local ecological environment, agricultural production, and social development [90].
Considering the amplitude of the oscillation, IMF1 showed a higher fluctuating amplitude in the 1970s and 2000s.Meanwhile, two distinct abrupt points were also detected by using the moving t-test method in the same period (Figure 3).This can indicate that IMF1 with the quasi-3-year cycle had a good capability to represent the drought abrupt characteristic.The periodic cycle with quasi-3-year in IMF1 was also consistent with that of precipitation records in the Tibetan Plateau [91].IMF2 with the quasi-9-year cycle showed a higher amplitude during the 1960s and 1980s, which might be influenced by the sunspot activity characterized by a 11-year periodic cycle [92].Due to less than 10 years of periodic cycles, IMF1 and IMF2 can be regarded as the interannual variations of the scPDSI.IMF3 and IMF4, with the lower amplitude and longer periodic cycle, could be used to reflect the interdecadal variation characteristics in drought, which have been considered to have a relationship with large-scale atmospheric circulations [93][94][95].Thus, the interdecadal variation indicated by IMF3 and IMF4 is a good reference for the long-term variation pattern of drought, which cannot be neglected in drought monitoring and warning.

Contribution Rates of IMFs and the Trend Component
To investigative the multi-scale variation characteristics of the scPDSI in the YZR basin, the inter-annual and inter-decadal scales were reconstructed based on the IMFs with different periodic cycles and the trend component.IMF1 and IMF2, with quasi-3-year and quasi-9-year cycles, respectively, reflected the inter-annual variation features in the original scPDSI time series.Hence, the inter-annual scale of scPDSI was reconstructed by the summation of IMF1, IMF2, and the trend component.Similarly, the inter-decadal scale of scPDSI was reconstructed by summing IMF3 (quasi-15-year), IMF4 (quasi-56-year), and the trend component.The results are shown in Figure 6.The reconstructed inter-annual scale (pink line in Figure 6) was consistent with the original scPDSI anomaly (black line in Figure 6) with a high Pearson correlation coefficient (R > 0.85), manifesting that the inter-annual variation was the dominant part in the original scPDSI anomaly.However, the Pearson correlation coefficient between the inter-decadal and the original scale was only 0.63, implying that the reconstructed inter-decadal scPDSI variation could not adequately portray the characteristic of the original scPDSI anomaly.This might be attributed to small-scale oscillations excluded from the reconstructed inter-decadal scPDSI anomaly variation [96].However, the reconstructed inter-decadal scale could capture the long-term variation process of the original scPDSI anomaly, and, for example, the YZR basin in the 1960s, 1980s, and 2000s exhibited a drying process with the decreasing scPDSI anomaly at inter-decadal scale, while the scPDSI showed an increasing tendency during the other periods, indicating the YZR basin became wetter.The effect of the fluctuated frequency and amplitude at each scale on the general characteristics of the original data can be expressed as the variance contribution rate.Table 2 shows the variance contribution rate of each component for the scPDSI anomaly.IMF2 had the maximum variance contribution (28.9%), followed by IMF1 (26.4%),IMF3 (21.3%), and the trend component (15.6%).The minimum variance contribution was found in IMF4 with the value of 7.8%.What is more, the trend component contributed up to 15.6% of the variance, indicating that the overall mean annual scPDSI in the YZR basin from 1956 to 2015 showed a nonlinear rise with wetting tendency, which was more obvious before the mid 1990s.Hence, it can be concluded that the inter-annual variation had a contribution of 55.3%, whereas the inter-decadal variation had only a contribution of 29.1%.

Spatial Variation of Drought
The mean annual value of the scPDSI anomaly from 1956 to 2015 was calculated for each grid.Meanwhile, the standard deviation (~0.60), calculated by using the mean annual values over the whole basin, was used as the classification criterion of drought in the YZR basin during the study period.The spatial pattern of the mean annual scPDSI anomaly is shown in Figure 7.The mean annual values of scPDSI anomaly > 0.3 were defined as a wet condition, while the mean annual values of scPDSI anomaly < −0.3 were defined as a dry condition.The wet regions accounted for 26.7% of the YZR basin, and were mainly concentrated in the north of midstream and downstream basins, especially in the headwater basin of Lhasa river and Parlung Tsangpo river.However, the arid regions accounted for 33.4% of the YZR basin, which were mainly located in the Nyang Qu river basin of the upstream and Bomi county of the northeast downstream basin.The Sen's slope and Mann-Kendall significance test were applied at each grid to detect trend variations from 1956 to 2015 in the YZR basin.Figure 8 shows that more than 95% of the grids experienced a wetting trend, which was consistent with the results of the linear regression method (Figure 2).The regions with the significantly wetting trend, accounting for about 48% of the grids, were mainly sited in the arid upstream and midstream regions.However, there existed a drying trend in some parts of the humid downstream regions, especially in the southern part, which was consistent with the results obtained by Wu et.al [97], demonstrating a frequent risk of drought in the east YZR basin.It can be concluded that the YZR basin was undergoing a wetting process due to the significant wetting trend in most parts of the basin.In addition, the wetting tendency occurred in the arid and semi-arid regions and the drying tendency occurred in the humid southern part, indicating a DWDW (dry gets wetter, wet gets drier) pattern in the YZR basin.

Multi-Scale Response to Climate Change
The Earth's climate system is considerably complex.Moreover, nonlinear and nonstationary processes with different periodic oscillations and scales are often influencing climate factors, as also noted in this study.It is meaningful for regional economic development and ecological protection to investigate the nonlinear response of drought to climate factors.The EEMD method can reveal hidden intrinsic non-stationary oscillation structures, and it was applied to decompose the long-term climate factor variations in this study.Precipitation (PRE), mean temperature (TEM), potential evapotranspiration (PET), maximum temperature (TMX), minimum temperature (TMN), and diurnal temperature range (DTR) were selected to explore the multi-scale response of drought to climate change.According to the reconstructed principal, all six climate factors were reconstructed at inter-annual and inter-decadal scales (Figure S1).The correlation coefficients between the scPDSI anomaly and different climate factors are illustrated in Figure 9.
In terms of original time series, precipitation had the closest positive correlation with drought variation with a correlation coefficient of 0.76 (Figure 9a), which indicated that precipitation was considerably important in drought variation processes.Similar results were also found in other regions [98,99].Zhang et al., (2010) concluded that precipitation was the main factor causing drought in southwestern China [98].Liu et al., pointed out that the Tibetan Plateau exhibited a wetter state, which could be mainly influenced by increased precipitation [99].Furthermore, precipitation exerted important influences at inter-annual and inter-decadal scales with correlation coefficients of 0.69 and 0.82, respectively.
A previous study showed that PET was a decisive drought-induced climate factor in semi-arid and arid regions [100] (i.e., consistent increased PET will induce more severe drought with lower scPDSI value).However, there was a small difference in our study.Although PET played a negative role in drought variation, the correlation coefficient between scPDSI and PET was relatively weak (−0.32, Figure 9a).Increased precipitation intensity has been observed in the YZR basin [101], which might be partly responsible for the relatively non-significant negative correlation between the scPDSI and PET [102].In addition, ET regulates the water and energy cycle and will keep a relatively balanced status, which is more evident in high altitude regions [103].That is, wetter conditions with increasing scPDSI values can supply more water, which is beneficial for higher ET, and higher ET can conversely consume more moisture, which will result in drought.That may explain why ET and drought showed a non-significant correlation in the YZR basin.
Long-term scPDSI anomaly variation in the YZR basin was strong negatively correlated with DTR, by a correlation coefficient of −0.70 in the original time series (Figure 9a).Some studies concluded that changes of PET were dominated by temperature variables [104][105][106].Although DTR is an important indicator of temperature, it showed a non-significantly correlation with changes of PET in the YZR basin as their correlation coefficient was only 0.36.However, in comparison to the weak correlation coefficient (−0.32) between PET and scPDSI, DTR showed the directly negative influence on the drought variation with the correlation coefficient up to −0.70, which indicated that DTR played a more complicated role in long-term drought variation rather than simply influencing the PET.This phenomenon can be attributed to the special geographical environment in the YZR basin.Due to the widely distributed glaciers and snow cover, the snow melting process in the cryosphere cannot be neglected and might change the surface albedo, latent heat, and other variables [107], which will contribute most to the DTR variation, and; meanwhile, also have great impacts on drought variation.Hence, an increased cloud thickness will contribute mostly to the decreasing DTR [108,109].Clouds can improve long-wave radiation and reduce short-wave radiation, which will further alter the energy balance and exert some influences on drought variation.Plus, this changing process in energy balance will be more obvious in the cryosphere over the YZR basin distributed with extensive snow cover and glaciers.In general, more clouds will lead to abundant precipitation [110,111], which also explained the cause that there was a significantly negative correlation between DTR and PRE over the YZR basin in the inter-annual (−0.56, Figure 9b) and inter-decadal (−0.67, Figure 9c) time scales.Furthermore, the significant increasing TMN was the main reason for the decrease of DTR, and the correlation coefficient between DTR and TMN was significantly negative due to the fact that the increasing slope with 0.023 • C/yr of the TMN was about twice that of the TMX with the slope of 0.012 • C/yr (Figure S2).

Discussion
In this study, correlation coefficients between scPDSI and TEM, TMX, and TMN, were very low (0.05, −0.18, and 0.28 in Figure 9a, respectively), which indicated that drought variation in the YZR basin was not predominantly determined by the temperature changes [112].This seems inconsistent with the previous study that demonstrated that rising temperature was the main reason for global drought, which might be attributed to the spatial scale inconsistency between global and regional scale levels and the specific geographical environment in the YZR basin.As mentioned above, glaciers, snow cover, and permafrost are distributed in the YZR basin, and the melting process in the cryosphere may play a considerable important role in long-term drought variation.Additionally, the correlation coefficient between scPDSI and TMX was −0.18, indicating that increasing TMX will induce more severe drought.Whereas TEM and TMN showed positive but little impacts on scPDSI variation, with correlation coefficients of 0.05 and 0.28, respectively, which can be attributed to more shallow soil moisture as a consequence of the accelerated melting process in cryosphere induced by climate change [113,114].Li et al. (2019) used GLDAS data to investigate soil moisture variations from 1970 to 2009 and concluded that, under global warming, the significantly increasing surface air temperature greatly influences the depletion of soil moisture and becomes the dominant controlling factor in soil moisture variations [115].These further indicate the important role of the cryosphere components in environmental changes over the YZR basin.
Correlations among scPDSI and climate factors at the inter-annual and inter-decadal scales (Figure 9b,c, respectively) were similar with those in the original time series (Figure 9a).However, there were some differences at the other two time scales.At inter-annual and inter-decadal scales, PRE and DTR still showed pronounced impacts on the drought variation process.Whereas, compared with the original time series, correlation coefficients between scPDSI and PRE and DTR at the inter-annual scale decreased (i.e., changing from 0.76 and −0.70 to 0.69 and −0.63, respectively).While the correlation coefficients strengthened from 0.76 and −0.70 to 0.82 and −0.83, respectively, at the inter-decadal scale.These indicated that PRE and DTR played a more important role at the inter-decadal time series, which may be linked to the atmospheric circulation factors in the longer term [116].Except for PRE and DTR, TMN had the most significant influence on the drought variation at the inter-decadal time scale, which portrayed long-term impacts of the snow melting process under global warming.

Conclusions
In this study, the spatiotemporal variation characteristics of drought and associated multi-scale response to climate change were investigated, based on the high accuracy scPDSI data and CRU data, from 1956 to 2015 over the YZR basin, located on the southern edge of Qinghai-Tibet Plateau, China.The main conclusions were summarized as follows.
(1) The YZR basin experienced a basin-wide wetting process.However, there existed an obvious transition point in the mid 1990s, which indicated that the YZR basin exhibited first a wetting and then a drying process.The drying process after the mid 1990s was harmful to the ecological environment, agricultural productivity, and social development.
(2) Most regions in the YZR basin, located in the arid upstream and downstream areas, showed a significant wetting tendency, estimated by the Sen's slope and Mann-Kendall significance test.Whereas regions showing drying tendency were mainly distributed in the humid southeastern part, which was inconsistent with the DDWW pattern, emphasizing the complexity of the YZR basin.
(3) The results of EEMD decomposition with significance test manifested that the scPDSI in the YZR basin had quasi-3-, quasi-8-, quais-15-, and quasi-56-year cycles.The reconstructed inter-annual scale with the variance contribution of 55.3% has good capability of representing the drought abrupt change characteristic, while the inter-decadal scale with the variance contribution of 29.1% can portray the long-term drought variation process.
(4) The multi-scale response of drought to climate change indicated that changes of precipitation and diurnal temperature range were the dominat driving factors in the drought variation at different time scales.Compared with potential evapotranspiration, diurnal temperature change was a more important drought-induced climate factor, which indicated the snow melting process of the cryosphere cannot be neglected in long-term drought variation over the YZR basin distributed with the glaciers, snow, and permafrost.
Drought will bring devasting influence on agricultural production, ecological restoration, and socio-economic development for the Yarlung Zangbo river basin, which is the political, economic, and cultural center in the Tibet Autonomous Region.Therefore, this study investigated the spatiotemporal variation of drought and discussed its multi-scale response to climate change in this region, providing certain implications for drought identification, monitoring, and adaption management.However, there are still some deficiencies.In the future, better quality and higher-resolution gridded data should be applied to further describe the spatial heterogeneity in the YZR basin, with complex terrain and climate.In addition, a distinct dry-wet transition point was detected in the mid 1990s, indicating that the water and energy cycle may be altered, which can redistribute hydroclimate variables and result in great influence on the terrestrial ecosystem.These changes deserve to be investigated deeply in the YZR basin.

Figure 1 .
Figure 1.The geographical location of the Yarlung Zangbo River (YZR) basin in Qinghai-Tibet Plateau and in China.

Figure 2 .
Figure 2. The inter-annual variation of the scPDSI anomaly in the YZR basin.

Figure 3 .
Figure 3.The abrupt detection of scPDSI with the moving t-test method in the YZR basin.The dotted purple line is 95% confidence level.

Figure 4 .
Figure 4.The original scPDSI anomaly, intrinsic mode functions (IMFs), and trend component in the YZR basin.

Figure 5 .
Figure 5.The significance test for the IMFs of the scPDSI anomaly in the YZR basin.

Figure 6 .
Figure 6.The original, inter-annual, inter-decadal, and trend component of scPDSI time series.

Figure 7 .
Figure 7.The spatial pattern of the scPDSI anomaly in the YZR basin.

Figure 8 .
Figure 8.The spatial pattern of the Sen's slope of the scPDSI anomaly in the YZR basin.The red dot represents that the trend is significant at the 95% confidence level.

Figure 9 .
Figure 9.The correlation coefficient matrix in the original time series (a), inter-annual (b), and inter-decadal (c).

Funding:
This work was jointly supported by the National Key R&D Program of China (No. 2018YFC1508702) and the National Natural Science Foundation of China (Grant No. 91647202, 41890822 and 51509247).

Table 1 .
Classification of the scPDSI (self-calibrating Palmer drought severity index) values.

Table 2 .
The contribution rates of IMF components and trend component of the scPDSI anomaly.