Abstract

Understanding the evolution and propagation of different drought types is crucial to reduce drought hazards in arid and semiarid regions. Here, Standardized Precipitation Index (SPI), Streamflow Drought Index (SDI), and Vegetation Condition Index (VCI) were used to investigate the spatiotemporal variation of different drought types and correlations between Pre (Pre-R)/post (Pos-R)-reservoir. Results showed that the average peak/intensity/duration/severity of meteorological droughts (MD) were greater in the Pre-R than in the Pos-R period in the upstream Heihe River Basin (UHRB), while there was little change between the Pre-R and Pos-R periods in the midstream Heihe River Basin (MHRB). The average peak/intensity/duration/severity of hydrological drought (HD) decreased in the mainstream for Yingluoxia (Ylx) but increased for Zhengyixia (Zyx) station in the Pos-R period. Propagation time decreased by 3 months (negative effect) in Ylx and increased by 8 months (positive effect) in Zyx compared with the Pre-R period. In the Pos-R period, propagation time increased (1–3 months) for tributaries (positive effect). Propagation times for the mainstream and tributaries varied for different seasons and time periods. Pearson’s correlation coefficient values were lower at short timescales (1–3 months) but higher at long timescales for the Pos-R period in Ylx and Zyx for SDI-1 with different timescales of SPI. The SDI and SPI had no lag in the UHRB and MHRB. However, VCI with SPI had a significant lag correlation at short timescales in the UHRB (lag 6 months) and MHRB (lag 4 months), and the VCI with SDI had a significant lag correlation for 1 month in the MHRB. The propagation time from MD to HD has been reduced for Pos-R in the UHRB. There was a positive effect (prolonged MD propagation HD time) in Pos-R but still faces serious drought stress in the MHRB.

1. Introduction

Drought is recognized as the world’s most costly and pressing natural hazard, giving rise to significant losses in the fields of the economy, ecology, and environment (e.g., crop losses, degradation and desertification, urban water supply shortages, and forest fires) [14]. The damage from droughts is expected to increase in severity [5]. Compared with other natural hazards, the spatial extent of drought is very large and its time of influence is commonly much longer. Under global warming, more frequent and severe droughts have been projected in the 21st century [68], particularly in the mid-latitudes [9, 10]. Severe and extreme droughts have become more frequent in the eastern and central portions of northern China since the late 1990s [11, 12].

In general, drought can be classified into four types: meteorological, hydrological, agricultural, and socioeconomic [13]. Meteorological droughts (MD) are the water shortages caused by an imbalance in precipitation and evaporation; precipitation is commonly used for MD analysis [14, 15]. Hydrological droughts (HD) are related to a period with inadequate surface and subsurface water resources for established water uses in a given water resource management system. Streamflow data are widely applied for HD analysis [14, 16]. Because the influence of drought on terrestrial ecosystems is becoming increasingly acute, the vegetation response to drought is a crucial topic in the domain of climate research [1719]. The Normalized Difference Vegetation Index (NDVI) is a good measure to estimate green biomass, leaf area index, and patterns of productivity and has been widely used to assess vegetation degradation, ecosystem features, and the physiological drought conditions of vegetation [20, 21].

Drought indices are the best ways to monitor drought and drought events at present. Lloyd-Hughes [22] counted over 100 drought indicators that have been developed for different types of drought [23]. A commonly used index is the Standardized Precipitation Index (SPI) [24]. Its calculation is simple because it only needs precipitation as an input, and it can be calculated over different timescales, which allows SPI to monitor both short-term droughts, such as agricultural droughts, and long-term droughts, such as HD [3, 2527]. For HD index, the Streamflow Drought Index (SDI) proposed by Nalbantis and Tsakiris [28] overcomes the problems of predicting drought onset and duration using cumulative streamflow volumes and the areal extent of drought using a spatially integrated streamflow at a basin outlet; it is also widely used in different regions of the world [2833]. The Vegetation Condition Index (VCI) is a normalization of the NDVI, which filters out the contribution of local geographic resources to the spatial variability of the NDVI [17, 34]. It reflects agricultural (vegetation) drought accurately in many parts of the world and has been applied in daily drought monitoring by the National Oceanic and Atmospheric Administration of the United States (NOAA) and the National Satellite Meteorological Centre of China [35].

Many researchers in related fields have focused on drought evolution, drought prediction, and developing reliable drought indices using observed and simulated data [4, 3639]. Recently, studying the propagation time of VCI and HD based on MD using Pearson’s or Spearman’s correlation coefficients has become a research focus [8, 10, 18, 19, 21, 23, 4046]. These two kinds of droughts can reflect the different stages of drought development. Commonly, MD develops and ends relatively quickly, while HD is the result of MD [4]. The corresponding propagation time depends on local landscape conditions [46]. Therefore, investigating the propagation time from MD to HD [4] and from MD to VCI along with its potential influencing factors is necessary for establishing an effective monitoring and warning system for HD and VCI based on MD.

However, little research has been done into the propagation of MD to HD and VCI and HD propagation to VCI, in particular in subarid and arid watersheds, which is one of the aims of this paper. The investigation of VCI and HD propagation time from MD at different timescales under the influence of global warming and human activities (e.g., construction of reservoirs or dams, land use, and land cover change) will help to elucidate the impacts of droughts on terrestrial ecosystems and inform planning and optimizing water resource allocations during droughts. It is important to identify the drought propagation process and mechanisms to help establish a drought early warning system. The primary objectives of this study are as follows: (1) to inspect the spatiotemporal characteristics of MD and HD duration, intensity, severity, and peak, while detecting the spatiotemporal characteristics of VCI; (2) to investigate VCI and HD propagation times based on MD at different timescales and VCI propagation times based on HD at different timescales; and (3) to explore the potential factors influencing drought characteristics and propagation times.

2. Materials and Methods

2.1. Study Area

The Heihe River Basin (HRB) is the second largest inland river basin in the arid zone of northwest China. The HRB flows through Qinghai Province, Gansu Province, and the Inner Mongolia Autonomous Region. The HRB is located along the land Silk Roads (“One Belt, One Road”) between 97.1°E–102.0°E and 37.7°N–42.7°N with a total area of approximately 1.43 × 105 km2 [47, 48]. The HRB includes three sections from south to north: the upstream Heihe River Basin (UHRB) from the Qilian Mountains to Yingluoxia station (outlet of the mountains) belongs to the subhumid and semiarid temperate continental monsoon climate zone; the midstream Heihe River Basin (MHRB) running from Yingluoxia station to Zhengyixia station belongs to the arid climate zone; and downstream HRB terminates in the Juyan Lakes (east and west branches, respectively) [49]. From upstream to downstream, the elevation decreases, the water availability decreases, and the landscape changes from glaciers and alpine biomes in the UHRB to steppes and agricultural ecosystems in the MHRB to riparian ecosystems and vast areas of desert in the downstream area [47]. In this study, the UHRB and MHRB were selected as the study area. Figure 1 shows the study area and the hydrometeorological stations. Analyses were also conducted at the sub-basin level in the UHRB and MHRB. In the MHRB, the annual average precipitation, potential evapotranspiration, and streamflow were 226 mm, 924 mm, and 88 mm, respectively. Over 90% of the population, grain production, and major industries are concentrated in the MHRB. Approximately 84% of the total available water was consumed for irrigation. On average, a decline in groundwater of about 1.86 m is attributable to water consumption for irrigation of the farmland area in Heihe River from 1981 to 2010 [50] with demand constantly increasing [51, 52]. There are a series of reservoirs (one large reservoir, 9 medium size reservoirs, and 89 small reservoirs) and a relatively complete irrigation system consisting of more than 893 main canals and branch canals [53, 54]. Land use and land cover have also significantly changed in the HRB. Urban, cultivated, and forest land have increased, while water area and grassland have become increasingly degraded.

2.2. Data

Daily precipitation datasets from 1961 to 2012 were collected by five meteorological stations (Table 1) and derived from the China Meteorological Data Sharing Service System of National Meteorological Information Centre V3.0 (http://www.nmic.gov.cn/). Rigorous quality control was conducted by the China National Meteorological Information Centre before the data were released. Monthly observed streamflow data, which were also of high quality, for eight hydrological stations (Table 2) were collected from the Hydrological Bureau of Gansu Province. Yingluoxia (Ylx) and Zhengyixia (Zyx) stations were located in the main stream of the HRB. The monthly NDVI (1 km spatial resolution) dataset was the satellite remote sensing data of SPOT/VEGETATION NDVI and MODIS from 1999 to 2012. The data were obtained from the resource and environmental science data center of the Chinese Academy of Sciences (http://www.resdc.cn), and LULC data in the UHRB (2000 and 2010) and MHRB (1975, 1985, 1995, 2000, and 2010) were obtained from the Cold and Arid Regions Science Data Centre (http://westdc.westgis.ac.cn/). The historical disaster records were derived from the China Meteorological Disaster Yearbook (Gansu volume). Zhang et al. [49] used Pettitt tests to determine abrupt change points in the annual streamflow series from 1945 to 2012 and from 1957 to 2012 at the Ylx and Zyx stations to ascertain the impact of the reservoir activities in the HRB. The results indicated only one significant abrupt change point in 1979 at Ylx station and in 1984 at Zyx station [49]. The impact of human activities is more evident after the Chinese Government initiated the Ecological Water Diversion Project (EWDP) in mainstream Heihe River for 2000 [54]. Therefore, in this paper, we focused on the mainstream Ylx and Zyx stations, between the Pre-R and Pos-R periods, before and after EWDP. The abrupt change points for the tributaries for six stations are shown in Table 2 (Pettitt tests).

2.3. Methods
2.3.1. Standardized Precipitation Index (SPI)

The SPI index was developed by McKee et al. [24]. Within a certain geographic area, the precipitation usually fluctuates regularly. If the precipitation is less than the average annual precipitation, a drought may occur. If precipitation exceeds the annual average, flooding may occur. To calculate the SPI, a frequency distribution function is first constructed from a series of long-term precipitation observations. In this paper, a three-parameter log-logistic distribution was used because the frequencies of precipitation accumulated at different timescales are well modeled using these statistical distributions. Because SPI is based on the cumulative probability of a given timescale, here the total amounts of precipitation in the current month and previous i months (i = 1, 2, 3, ......) were used to calculate SPI on a timescale of i + 1 month. In our study, we used the SPI (1, 3, 6, 9, ......, 48 months) to analyze the characteristics of drought change in the HRB. For drought classifications, please see Table 3.

2.3.2. Streamflow Drought Index (SDI)

The SDI was developed by Nalbantis and Tsakiris [28] to characterize hydrological drought based on developing the concepts of the SPI. It is worth noting that the procedure should be performed separately for each month and can be calculated for any timescale (e.g., 1, 3, 6, 9, and 12 months). The total streamflow in a given month j and year i depends on the chosen timescale k and can be calculated as [55]where Vi-1,l and Vi,l denote streamflow volumes in years i − 1 and i, respectively. Due to the three-parameter log-logistic distribution used to calculate the SPI, we also chose this probability distribution to standardize the observed streamflow data used to get the SDI to facilitate comparison in the HRB. It is emphasized here that positive SDI values reflect wet conditions, while negative values indicate a hydrological drought.

2.3.3. Vegetation Condition Index (VCI)

Kogan and Sullivan [56] proposed the VCI index in 1993. It is here assumed that the change of the VCI index is only affected by weather factors; in this way, the ecosystem factors and extreme weather factors are separated, the impact of extreme weather on vegetation is evaluated, and agricultural drought is monitored. A smaller VCI index correlates to a greater drought. The VCI index is calculated as follows:where NDVI, NDVImax, and NDVImin are monthly NDVI, multiyear maximum NDVI, and multiyear minimum NDVI (in this study, 14 years), respectively, for each grid cell. The VCI changes from 0 to 100 corresponding to changes in vegetation condition from extremely unfavorable to optimal. In the case of an extremely dry month, the vegetation condition is poor, and the VCI is close to or equal to zero. A VCI of 50 reflects fair vegetation conditions. At optimal vegetation conditions, the VCI is close to 100 [34]. According to t-test, the trend was divided into extremely significant increase (slope > 0, p < 0.01), significant increase (slope > 0, 0.01 ≤ p < 0.05), insignificant change (p ≥ 0.05), extremely significant decrease (slope < 0, p < 0.01), and significant decrease (slope < 0, 0.01 ≤ p < 0.05).

2.4. Drought Variables and Run Theory

According to McKee et al. [24] and Spinoni et al. [57], the derived drought variables based on run theory [58] follow certain definitions (Figure 2). When the monthly SPI/SDI values were below −1, the corresponding month was potentially identified as a drought event (e.g., the four drought events for D0, D1, d0, and d2 shown in Figure 2). Thresholds of −1 (moderate drought), −1.5 (severe drought), and −2 (extreme drought) were used to identify drought events and drought characteristics. Drought duration was defined as the number of months from the first month in which the indicator goes lower than −1 to the last month with a negative value before the indicator returns to positive values. If the duration of a drought event only contained one month where SPI/SDI < −1, this month was regarded as a single drought event (e.g., D1 with a magnitude S1). Drought severity was defined as the sum of the monthly absolute values of the index when the index was ≤−1.0 (−1.5, −2). A drought event may contain a few consecutive months with negative SPI/SDI (e.g., duration D0 and its magnitude S0). If the duration of a drought event contained two branches, such as d0 and d2, and the interruption period d1 between d0 and d2 was less than 6 months in which 0 < SPI/SDI < 1 [59, 60], these months were still regarded as a single drought event (D2 = d0 + d1 + d2). The corresponding magnitude was then defined as S2 = s1 + s2 [59, 60]. Intensity stands for the number of months in which the drought indicator was lower than −1 (−1.5, −2). Drought peak referred to the month in the “drought event” with the lowest value of the indicator [61].

2.5. Correlation Analysis

Pearson’s correlation coefficient (PC) was calculated between SDI-1 (one month SDI) and SPI (1, 3, 6, 9, ....., 48 months), monthly (April–October) VCI, and different timescales of SPI and SDI (1, 3, 6, 9, and 12 months) series indicates the characteristics of interannual propagation at 95% and 99% confidence levels. The PC ranges from −1 to 1; a PC > 0 indicates a positive correlation and a PC < 0 indicates a negative correlation. In this study, the PC method was also used to detect the intra-annual propagation time of the SDI-1 to the different timescales of SPI and monthly VCI to the different timescales of SDI and SPI conditions. For instance, if the maximum correlation in the March SDI and the March series of SPI at the timescales of 1, 3, …, 48 months is recorded at the 3-month timescale, this means that the SDI in March is mostly determined by the cumulative water balance occurring from January to March [19, 21, 23, 41]. Meanwhile, to determine whether there is a lag between the SPI (accumulation periods of 1–48 months) and SDI-1, cross-correlations were calculated for SDI-1 series which were lagged by 0 to 6 months after the SPI series. In this case, the SPI accumulation period with the strongest correlation with SDI-1 was denoted as the lagged SPI-n (1, 3, 6, 9, ....., 48 months) [23]. The same interpretation applies to the lag correlation of the VCI with different timescales SPI and VCI with different timescales SDI.

3. Results

3.1. MD and HD for the Mainstream Characteristics between Different Timescales and Thresholds

The distribution of drought characteristics based on thresholds of −1, −1.5, and −2 for different timescales at different periods, including Pre-R (1965–abrupt change year), Pos-R (abrupt change year +1–2012), and total period (1965–2012), was evaluated in the UHRB and MHRB using boxplots. The drought characteristics of SPI-1 and SDI-1 in the MHRB (Table S1) were almost consistent compared with historical records and, for different combinations of timescales, thresholds and periods (Table S2) are shown in the supporting information. The drought peak and intensity (Figures S1-S2), duration, and severity for different timescales at different thresholds are shown in Figure 3 for MD and in Figure 4 for HD. More drought events were identified at shorter accumulation periods (1–12 months for MD and 1–6 months for HD) based on thresholds equal to −1.0, −1.5, and −2.0 in the UHRB (Ylx) and MHRB (Zyx). For MD, the average intensity, duration, and severity of fluctuations increased with cumulative timescale increases between Pre-R and Pos-R periods. The changes were more obvious in Pos-R compared with Pre-R periods for all thresholds in the UHRB. The average peak, intensity, duration, and severity of the MD did not change much for Pre-R, whereas Pos-R fluctuations increased with cumulative timescale increases using thresholds of −1.0, −1.5, and −2.0. The average peak, intensity, duration, and severity of value changed little in the MHRB. With regard to HD, the duration and severity were larger for HD than for MD because a short-term MD may not develop into an HD [42] and as the number of events decreases, duration lengthens, and severity worsens. The average intensity (Figure S2), duration, and severity of the HD were significantly different in Ylx and Zyx for different periods based on the thresholds −1.0, −1.5, and −2.0. In Ylx, the characteristics of HD were similar to MD in that, compared to the Pre-R period, the average intensity, duration, and severity of the HD decreased significantly, especially for SDI-3 and SDI-12. However, in Zyx, the HD were quite inconsistent with MD, and the average intensity, duration, and severity of the HD increased significantly during the Pos-R period in Zyx, especially for SDI-1, SDI-6, and SDI-12.

3.2. Comparison of MD Characteristics for SPI-1 and HD for SDI-1

The characteristics of MD and HD based on SPI-1 and SDI-1 using thresholds of −1, −1.5, and −2 in the UHRB and MHRB are shown in Tables 4-5 and Figure 5. More MD events occurred at Yng station (names and abbreviations are shown in Table 1) in the UHRB (Figure 5) and at Sd station in the MHRB. More HD events occurred at Zyx (Figure 5). It can be seen from Tables 4-5 that the average, maximum, and minimum peak/intensity/duration/severity in Zyx, Yng, and Ql stations were identified as longer and more severe based on thresholds of −1.0, −1.5, and −2.0 for MD and HD. HD characteristics based on these thresholds for the average peak, maximum peak, and minimum peak in Sh, respectively, were more obvious, and longer average, maximum, and minimum intensities were found in Ylx, Ljq, and Sh (Table 5). Longer, more severe HD occurred in Ylx, Ljq, Sh, and Wfc.

3.3. Variation in VCI

Vegetation condition in the extremely dry months of April, May, and October (Figure 6) was poor, especially for alpine vegetation, cultivated vegetation, grassland, and marshy grassland. From June to September, the vegetation was in optimal condition for coniferous forests, shrubs, and marshy grassland in the UHRB and for coniferous forests and cultivated vegetation in the MHRB. It is worth noting that the vegetation condition was poor for alpine vegetation over the whole growing season.

We further analyzed the trend and significance level of the growing season for VCI. Figure 7 shows that, over the entire growing season, alpine vegetation (UHRB) and grassland (MHRB) showed a significant decreasing trend, while cultivated vegetation (MHRB), coniferous forest, shrubs, and grassland (UHRB) showed a significant increasing trend.

3.4. Correlation of VCI, HD, and MD
3.4.1. Interannual Propagation Time between VCI, HD, and MD at Different Timescales in Different Periods

Figures 8-9 describe the relationship between SDI-1 and VCI-1 with all SPI accumulation periods (1–48 months) and VCI with all SDI accumulation periods (1–12 months) for different periods (Pre-R and Pos-R, before and after EWDP, and total period) in the mainstream for Ylx (SPI calculated using the arithmetic mean of precipitation from Yng and Ql stations), Zyx (SPI calculated using the arithmetic mean of precipitation from Gt, Sd, and Zy stations), and tributaries (SPI calculation is based on the area flowing through or the adjacent area because all meteorological stations were not available, e.g., Sh-Sd, Sss-Ql, Lyb-Yng, Ljq-Sd, Wfc-Ql, and Syk-Zy). Figure 8 indicates that the propagation times were 1, 1, 3, 1, 1, and 1 in the Pre-R period, respectively; 1, 3, 9, 3, 6, and 6 in the Pos-R period, respectively; and 1, 3, 6, 1, 6, and 6 in the total period for six tributaries, respectively. Figure 9(a) shows that the SDI-1 and SPI had the highest PC (PC = 0.27, ) in the Pre-R period for the 9-month period, and, in the Pos-R period, the highest PC (PC = 0.46, ) was observed for the 6-month period. Figure 9(a) also shows the highest PC for SDI-1 and SPI (PC = 0.38, ) was for the 6-month period before the EWDP period, and, after the EWDP period, the highest PC (PC = 0.48, ) was observed for the 6-month period in Ylx. Figure 9(b) reveals the highest PC (PC = 0.19, ) was found for the 1-month period Pre-R and the 9-month period (PC = 0.24, ) Pos-R. Figure 8(b) also demonstrates the highest PC (PC = 0.17, ) was found for the 3 months before the EWDP period and over 36 months (PC = 0.35, ) after the EWDP period in Zyx. These results demonstrate that the changes in the streamflow were more sensitive to precipitation in the Pre-R period for short timescales (1–3 months), and Pos-R propagation time increased (2–6 months) for tributaries, while the propagation time decreased by 3 months in Ylx and increased by 8 months in Zyx compared to the Pre-R period. Thus, reservoir operation increased the propagation time of the streamflow variables to the SPI during the Pos-R period in Zyx, whereas the propagation time of the streamflow variables to the SPI was decreased during the Pos-R period in Ylx. There was no change in Ylx before and after the EWDP, while after the EWDP, the propagation time of the HD to the MD increased on a 33-month timescale compared with before the EWDP in Zyx. Figure 9(c) shows that the VCI and SPI strongest negative PC was found for the 6-month (PC = −0.09) period in the UHRB, and the strongest positive PC was found for the 36-month (PC = 0.13) period in the MHRB. Figure 9(d) shows that the VCI and SDI highest PC was found for the 6-month period (PC = 0.23, ) in the UHRB and the 3-month period (PC = 0.43, ) in the MHRB. The results indicated that correlation between VCI and MD was not significant in the UHRB and MHRB. However, the correlation between VCI and HD was significant in the UHRB (Ylx) and MHRB (Zyx), while the VCI propagation time to HD was found at short timescales.

3.4.2. Intra-Annual and Lag Correlations between VCI, HD, and MD at Multiple Scales in Different Periods

The SDI-1 with SPI at different timescales was investigated in relation to the propagation from MD to HD. Figure 10 shows the propagation times from MD to HD had significant seasonal characteristics for tributaries in different periods. In spring (March–May), the drought propagation time was relatively short (1–12 months) in Sss, Ljq, and Wfc, whereas the propagation time was relatively long (1–48 months) in Sh, Lyb, and Syk. In summer (June–August), the drought propagation time was relatively long, except for Wfc. Drought propagation times in autumn (September–November) were relatively short (1–6 months) in Sss, Lyb, and Syk and relatively long (1–48 months) in Sh, Ljq, and Wfc. Drought propagation times in winter (December–February) were similar to those in summer and were relatively long (1–48 months), except for Lyb in the Pre-R period. Drought propagation times in autumn and winter were relatively long, whereas in spring and summer, they were relatively short in the Pos-R period. The total period was similar to the Pos-R period. To summarize, in the Pos-R period, the drought propagation times in winter did not change and remained long, whereas drought propagation times were increased in autumn. Propagation time was longer for Ss, Lyb, and Syk and shorter for Sss, Ljq, and Wfc in spring; and propagation times in summer were shorter, except for Wfc and Ss. It can be clearly seen from Figure 11 that the intra-annual correlations between the monthly SDI-1 and SPI at multiple scales and in different periods have strong seasonal characteristics in Ylx and Zyx. The propagation times in spring, summer, autumn, and winter were approximately 3, 9, 6 and 12 months in the Ylx for the Pre-R period (Figure 11(a)), respectively. The propagation times in the spring, summer, autumn and winter were approximately 6, 9, 1, and 6 months in Ylx for the Pos-R period (Figure 11(b)), respectively. The propagation times in spring, summer, autumn, and winter were 6 months in Ylx for the total period (Figure 11(c)). The propagation times in spring, summer, autumn, and winter were approximately 6, 1, 1, and 15 months in Zyx for the Pre-R period (Figure 11(d)), respectively. The propagation times in spring, summer, autumn, and winter were 9, 12, 36, and 9 months in Zyx for the Pos-R period (Figure 11(e)), respectively. The propagation times in spring, summer, autumn, and winter were approximately 9, 9, 3, and 6 months in Zyx for the total period (Figure 11(f)), respectively. These results show that the propagation time in the Pos-R period increased in spring, decreased in autumn and winter, and did not change in summer in Ylx compared to the Pre-R period, whereas in Zyx it increased in spring, summer, and autumn, and it decreased in winter. Compared with the Pre-R period, the PC was lower on a 1-month timescale in Ylx (Figure 11(b)) and on a 1–3-month timescale in Zyx (Figure 11(e)). However, the PC between SDI-1 and different timescale of SPI in the Pos-R period was higher than that in the Pre-R period for long-term timescales (>12 months) in Ylx (except for January and June) and Zyx (except for June, July, November, and December). It can be seen from Figures S3(a)-S3(b) that, after the EWDP, the PC of SDI-1 and the timescale of SPI were higher in spring and autumn and lower in summer (except for August) and winter than before the EWDP in Ylx. However, Figures S3(c)-S3(d) show that, after the EWDP, the PC of SDI and SPI was higher in spring, summer, autumn, and winter, specifically in spring, in Zyx. In different periods (Pre-R and Pos-R and before and after EWDP) and different timescales, the PC of the Pos-R period and after the EWDP was higher in Ylx (except for 1, 36, and 48 months) and Zyx, especially in Zyx than in the Pre-R period and before the EWDP. Figure 12(a) shows that the higher PC between VCI and SPI in April at 6 months was 0.63 (), whereas the lower PC in May at 1 month and 48 months was −0.06. The highest PC between VCI and SDI in May at 12 months was 0.78 (), whereas the lowest PC in March at 3 months was −0.69 () in the UHRB (Figure 12(b)). Figure 12(c) shows that the highest PC between VCI and SPI was 0.84 () in July at 9 months, whereas the lowest PC was −0.19 in October at 1 month. The highest PC between VCI and SDI was observed at 0.71 in May at 6 months, and the lowest PC was −0.23 in March at 6 months in the MHRB (Figure 12(d)).

We further analyzed whether there was a lag between the different timescales of SPI with SDI-1. The PC between SDI-1 lagged by 0–6 months, with SPI accumulation periods of 1–48 months for different timescales at different periods (Pre-R, Pos-R, and total), and the strongest correlation was detected at a lag of 0 months (i.e., no lag) in Ylx and Zyx (Figure 13). We also found no lag before and after the EWDP (Figure S5), whereas there was a higher PC at 3–12 months before the EWDP and after the EWDP at 6–15 months in Ylx. The highest PC was at 1–12 months before the EWDP and after EWDP at 9–48 months in Zyx. We also analyzed the lag correlation of the SPI with VCI-1 and SDI (accumulation periods of 1–12 months) with VCI-1 (Figure 14), and PC was calculated for the VCI-1 series, which were lagged by 0 to 6 months after the SPI and SDI series. The results demonstrated that the highest PC (PC = 0.21, ) was found at a 6-month lag between VCI and SPI, while the strongest correlation was at a lag of 0 months (i.e., no lag) between VCI and SDI in the UHRB (Ylx) during 1999–2012. The highest PC (PC = 0.26, ) was observed at the 4-month lag scale between VCI and SPI, and the strongest PC (PC = 0.47, ) was found at a lag of 1 month between VCI and SDI during 1999–2012 in the MHRB (Zyx). As mentioned above, it can be concluded that there was no lag between SDI and SPI; however, VCI with SPI had significant lag correlations in the short term in the UHRB and MHRB. Additionally, the VCI with SDI had a significant 1-month lag correlation in the short term in the MHRB (Zyx).

4. Discussion

Our results indicated that MD were more serious in Pos-R than in Pre-R in the UHRB, there was little change in the MHRB, and MD were more serious in the UHRB than in the MHRB for duration and severity based on thresholds of −1.0, −1.5, and −2.0 for the total period. These results are inconsistent with [44], probably because the precipitation data types and SPI probability distribution function were different. In our study, precipitation data types were based on meteorological stations, and the SPI probability distribution function was a three-parameter log-logistic distribution. However, Ma et al. [44] used 0.5° resolution precipitation data and their SPI probability distribution function was an empirical Gringorten plotting position. The Pos-R period had a significant positive impact on the evolution of different timescales based on thresholds of −1.0, −1.5, and −2.0 in Ylx, which was consistent with Wu et al. [60], who showed that Pos-R period in the HD had reduced duration and severity compared with Pre-R period. However, there were more serious duration and severity, specifically at 6 and 12 months, than in the Pre-R period in Zyx. For the total period (Table S2), HD events were more frequent at Zyx than at Ylx for different timescales based on thresholds of −1.0, −1.5, and −2.0. These results were consistent with Ma et al. [44]. Figure 1(c) and Table 6 show that the main vegetation type in the UHRB was grassland, and there was an increasing trend to cultivated land for 2000–2010. The UHRB had fewer reservoirs and less arable land than the MHRB, and the drought duration and severity were more likely to be related to climate change. For example, Ma et al. [44], Qiu et al. [52], Gao and Zhang [62], Yang et al. [63], Cong et al. [64], and Cheng et al. [65] showed that an increase in streamflow was caused by climate change (increased precipitation) factors rather than human activities in the UHRB. At Zyx, the duration and severity of drought were greater than those in the Pre-R period. More extreme MD events occurred in 1985, 1992, 2000, and 2005 (i.e., Pos-R) [66], and the decrease in precipitation may have led to a decrease in streamflow and thus more serious HD events compared with the Pre-R period (Table S2). The area of cultivated land increased, and unused land and grassland decreased from 1975 to 2010. The total irrigation volume in the area increased by 0.57 × 108 m3/year in 1990–2010 in the MHRB [54]. Along with the rapid increase in the number of pumping wells, groundwater irrigation has increased 2.11 × 108 m3/year [54]. It is clear that human activities have changed the land surface and altered hydrological processes, including evapotranspiration, infiltration, surface runoff, and storage of water and; in this way, they affect the development of drought and thereby exacerbate drought [67]. Zhu et al. [68] also pointed out that human activities, such as the expansion of irrigation, rapid population growth, and socioeconomic development, have deeply modified hydrological processes in the MHRB. Therefore, the abovementioned climate change and human activities may affect changes in drought characteristics between the different periods in the UHRB and MHRB.

Previous research has examined the propagation time and lags of HD to MD in different regions of the world [4, 23, 43, 44]. In this study, we found that the PC at Ylx was larger than that at Zyx. The correlation with climate change and human activities changed from positive to negative at Zyx. Ma et al. [44] found that climate change was inclined to increase streamflow and propagation time, contributing from −57% to 63% in the UHRB, whereas human activities played a dominant role in water consumption with a contribution rate greater than −89% to further alter HD characteristics and propagation time in the MHRB. In this study, we found the MD propagation to HD was different in different catchments and had a significant seasonal discrepancy at different timescales for different periods. The reservoirs increased the propagation time of the hydrological variables to MD during the Pos-R period at Zyx, and the results were in line with Wu et al. [60]. The PC values between SDI and SPI in the Pos-R period were lower on a short-term scale (1–3 months) than those in the Pre-R period, whereas they were higher for the long-term timescale (>12 months) at Ylx and Zyx, which was also consistent with Wu et al. [60]. The propagation time in the Pos-R period increased in spring, decreased in autumn and winter, and did not change in summer at Ylx compared with the Pre-R period, while it increased in spring, summer, and autumn and decreased in winter at Zyx. The possible causes of the uneven spatial variations of precipitation are temperature and melting of snow and glaciers in the Qilian Mountains when the temperature is rising, which may lead to the differences in HD variations [44, 62]. The propagation time for tributaries would also vary in different periods in relation to natural and social environmental factors [4, 23, 69]. For instance, Deng et al. [70] demonstrated that the 33 tributaries in the MHRB no longer joined the mainstream after 1980s, and they gradually disappeared and formed independent irrigation oases. These may diverge in terms of different climate and human activity impacts. The lack of lag in the SDI and SPI was in line with Barker et al. [23] at Ylx and Zyx.

The high correlation between monthly NDVI and the MD index at different timescales is an indicator that can describe the impact of drought on vegetation, and the month of highest correlation means the greatest sensitivity of vegetation to drought [10, 41]. Our results exhibited a negative correlation on short timescales (≤9 months) between VCI and SPI, consistent with Li et al. [19]. Meanwhile increasing drought timescales (≥12 months) suggested positive correlations in the UHRB and MHRB (Figure 8(c)). Niu et al. [45] also found that variations in the Enhanced Vegetation Index (EVI) generally followed the evolution of soil moisture in the top 1 m of the soil profile. There were both positive and negative time lags. For example, during a 2001 drought event, it was observed that the Soil Moisture Anomaly Index reached a negative peak value. However, the EVI did not show a consistent response, indicating it was not affected by natural drought events because it was maintained by irrigation activities in cropland in the MHRB. The strongest PC values between VCI and SPI were for 6-month and 36-month timescales for the growing season in the UHRB and MHRB, respectively. These results were inconsistent with Vicente-Serrano et al. [41], who found that, in arid and humid regions, vegetation responded to MD at short timescales as well as in semiarid and subhumid regions at long timescales. The main reason may be the relationship with different vegetation types, anthropogenic activity, and climate change, although the UHRB is a semiarid region. In Figure 1(c), there is little coniferous forest, and the region is mainly grassland, shrub, and alpine vegetation. However, cultivated vegetation dominates in the MHRB, and there are 24 irrigation districts (Figure S5) and 71 reservoirs (four under construction) with thousands of canals and over 6000 pumping wells for irrigation [71]. The total irrigation volume in the area has increased in the MHRB [52]. Thus, crop improvements and water-saving irrigation technologies would increase water use efficiency and drought tolerance capacity [70]. Our findings also differed from Zhang et al. [18], who found significant annual high PC values in arid and semiarid regions, and the corresponding drought timescales were 3–6 months (January–December). This may be related to the different timescales (April–October). We further found a significant positive and high PC value (PC ≥ 0.53, ) between VCI and SPI on an intra-annual scale, implying that water availability is the critical factor for vegetation’s various spatiotemporal activities [52] in the UHRB and MHRB, especially in the MHRB (PC ≥ 0.70, ). Our results also revealed that, with an increasing drought timescale, the PC values between VCI and SPI were higher in the MHRB, specifically in August (the PC values fluctuated from 0.43 to 0.82 for 1–48 months) and September (the PC values fluctuated from 0.07 to 0.83 for 1–48 months). These results suggested that croplands can suffer from extreme and prolonged drought conditions [10] in the MHRB. Furthermore, different vegetation types have different propagation times for MD; the grassland and cultivated vegetation at a 3-month timescale and a 12-month timescale for the shrubland, coniferous forest, and broadleaf forest [8, 21]. Figures 14 and 6 demonstrate that alpine vegetation and grassland were poor (especially in April and May) and there was an extremely significant decrease for the growing season, whereas cultivated vegetation was good, and there was an extremely significant increase for the growing season in the UHRB and MHRB. Additionally, VCI with SPI had a significant lag correlation in the short term in the UHRB (6-month lag) and MHRB (4-month lag), and the VCI with SDI had a significant 1-month lag correlation in the MHRB, similar to Lin et al. [35]. Additionally, the correlation between a single timescale of the HD and VCI with the different timescales of the MD did not fully reflect their relationships. The reason is that the streamflow and vegetation propagation times with precipitation differ at different timescales [10, 18, 19, 21, 72, 73].

Our study had some limitations. The quantification of drought is a difficult task and drought is intrinsically multiscalar. There is no unique physical variable we can measure to quantify drought intensity. Vegetation propagation time to droughts is still an open scientific problem due to the complexity of droughts and limited knowledge of physiological processes. However, understanding the relationship between these mechanisms and the characteristics of droughts is crucial for improving our knowledge of vegetation vulnerability to climate fluctuations and climate change [41]. Meanwhile, with global warming, the propagation time of VCI and HD, VCI and MD, and HD and MD to climate change also has changed [5], which poses new challenges that need further research. Furthermore, increasing human activity has dramatically changed natural droughts; to manage drought effectively, we need to acknowledge that human influence is as integral to drought as natural climate variability. The complex interactions between natural and human processes need to be considered [67]. In addition, there are still some issues worthy of in-depth discussion. Often, people are considered as passive recipients at the end of the propagation cascade from MD via soil moisture drought to HD [74] but, in fact, people actively influence drought propagation [67]. Anthropogenic changes to the land surface alter hydrological processes, including evapotranspiration, infiltration, surface streamflow, and storage of water [67, 75]. Understanding how irrigation [76] (the effect of long-term irrigation in the HRB on the evolution of drought [77, 78]), land use/land cover change [79], and construction of reservoirs or dams [42, 63, 64, 80, 81] impact MD, HD, and VCI characteristics and MD propagation to HD and VCI need further research. Trnka et al. [82] presented a set of 60 priority questions and considered interdisciplinary characteristics to optimize future drought research, primarily covering the following drought-related topics: monitoring, impacts, drought forecasts, climatology, adaptation, and planning. Therefore, in future drought research, it is necessary to comprehensively consider priority questions and interdisciplinary characteristics, especially for irrigated agricultural areas in arid regions.

5. Conclusions

Drought affects land surface dynamics, and the strongest impact of drought on vegetation occurs in arid and semiarid areas. Therefore, understanding the evolutionary characteristics and propagation of different drought types may provide useful information about appropriate adaptation and mitigation strategies against the effects of drought on agricultural production and vegetation production and even alleviate the losses caused by drought. Our results highlight the following points:(1)More drought events occurred at shorter accumulation periods based on thresholds = −1.0, −1.5, and −2.0 for MD and HD. Compared to the Pre-R period, the average peak, intensity, duration, and severity of HD decreased at Ylx, whereas the effects were more serious at Zyx. The PC values were lower at short-term scales (1–3 months) than at long-term scales for the Pos-R period at Ylx and Zyx. The PC values for MD and HD after the EWDP were higher than those before the EWDP in the UHRB and MHRB, especially in the MHRB. Human activities at different timescales (Pre-R period and Pos-R period and before EWDP and after EWDP) may affect the correlation between drought and the timescale of MD lag to HD.(2)The propagation time decreased 3 months at Ylx and increased 8 months at Zyx compared to the Pre-R period and between VCI and SPI at 6 months in the UHRB and 9 months in the MHRB. However, VCI with SDI at 6 months at Ylx and 3 months at Zyx did significantly affect the growing season. On an intra-annual scale, for Ylx, Zyx, and six tributaries, propagation times were different in the Pre-R and Pos-R periods, and there were seasonal differences between Pos-R and Pre-R periods.(3)The SDI and SPI showed no lag at Ylx and Zyx. However, VCI with MD had a significant lag correlation at the short-term scale in the UHRB (6-month lag) and MHRB (4-month lag), while the VCI with HD has significant 1-month lag correlation at Zyx.(4)Overall, the average peak/intensity/duration/severity of MD and HD are weakening and the propagation time from MD to HD is also reduced for Pos-R period compared with Pre-R period in the UHRB. Constructing reservoirs has had a positive effect, prolonging MD propagation to HD times, but the drought characteristics (average peak/intensity/duration/severity) have increased from the perspective of climatology. Therefore, positive drought prevention measures are necessary to consider the characteristics of propagation between MD and HD as well as seasonal differences in the MHRB.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest in this paper.

Acknowledgments

This research was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant nos. XDA19070502, XDA20100104, and XDA19040500), the National Natural Science Foundation of China (Grant nos. 41571516 and 41690144), and the Fundamental Research Funds for the Central Universities (Grant no. 2019jbkyjd013).

Supplementary Materials

Figure S1: boxplots showing MD peak and intensity in the UHRB ((a) and (c)) and in the MHRB ((b) and (d)) based on SPI using thresholds of −1, −1.5, and −2 for different timescales at different periods (the dot and midline in the boxplots are the mean and median; the upper and lower boundaries of box are the interquartile ranges). Figure S2: boxplots showing HD peak and intensity at Ylx ((a) and (c)) and at Zyx ((b) and (d)) based on SDI using thresholds of −1, −1.5, and −2 for different timescales at different periods (the dot and midline in the boxplots are the mean and median; the upper and lower boundaries of box are the interquartile ranges). Figure S3: the growing season spatial distribution of VCI during 2000–2001 in the UHRB and MHRB. Figure S4: heat map showing interannual Pearson’s correlation coefficient of SPI accumulation periods of 1–48 months with SDI-1 in the UHRB (Ylx) and MHRB (Zyx) (note: the x-axis represents the correlation between the SDI-1 and SPI at different timescales and the y-axis represents the 12 months of the year; the up arrow represents that the correlation coefficients are increased in the period after EWDP in the UHRB (Ylx) and MHRB (Zyx) compared to the period before EWDP, or vice versa); (a) (1965–1999, PC ≥ 0.35, ), (b) (2000–2012, PC ≥ 0.57, ) notes in the Ylx and (c) (1965–1999, PC ≥ 0.35 or ≤ −0.38, ), (d) (2000–2012, PC ≥ 0.56, ) notes in the Zyx. Figure S5: heat maps showing PC between VCI lagged by 0–6 months and SPI (1–48)/SDI (1–12) before (1965–1999) and after (2000–2012) EWDP; (a) (1965–1999, PC ≥ 0.10, ), (b) (2000–2012, PC ≥ 0.16, ) in the UHRM (Ylx) and (c) (1965–1999, PC ≥ 0.11 or PC ≤ −0.11, ), (d) (2000–2012, PC ≥ 0.17, ) in the MHRB (Zyx). Figure S6: reservoirs under construction and construction (a) irrigation system and area (b) in the MHRB. Table S1: identification of MD (threshold = −1, SPI-1) and HD (Zyx) characteristics (threshold = −1, SDI-1) based on run theory with comparison with historical drought records in MHRB. Table S2: MD and HD drought characteristics for different timescales for different threshold for prereservoir period, postreservoir period, and total period in UHRB (Ylx) and MHRB (Zyx). (Supplementary Materials)