# The Spatial-Temporal Variation Characteristics of Natural Vegetation Drought in the Yangtze River Source Region, China

^{1}

^{2}

^{*}

*Int. J. Environ. Res. Public Health*

**2021**,

*18*(4), 1613; https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph18041613

Next Article in Journal

Next Article in Special Issue

Next Article in Special Issue

Previous Article in Journal

Previous Article in Special Issue

Previous Article in Special Issue

Department of Geographical Sciences, Faculty of Resources and Environmental Science, Hubei University, Wuhan 430062, China

Changjiang River Scientific Research Institute, Changjiang Water Resources Commission of the Ministry of Water Resources of China, Wuhan 430010, China

Author to whom correspondence should be addressed.

Academic Editor: Elena Maria Otazo Sánchez

Received: 13 November 2020 / Revised: 26 January 2021 / Accepted: 30 January 2021 / Published: 8 February 2021

(This article belongs to the Special Issue Global Warming & Water)

In the context of climate change, ecosystem in Yangtze River Source Region (YRSR) is under threat from severe droughts. This study introduced a new natural vegetation drought index, standardized supply-demand water index (SSDI), and identified natural vegetation drought events and parameters (e.g., duration, severity, peak, and coverage area) based on run theory. Then the drought-prone regions were investigated via 2-dimensional joint copula. The results indicate that (1) compared with traditional meteorological drought index, the SSDI is reliable and can reflect the comprehensive characteristics of the ecological drought information more easily and effectively; (2) the YRSR had witnessed the most severe drought episodes in the periods of late-1970s, mid-1980s, and mid-1990s, but the SSDI showed a wetting trend since the mid-2000s. Additionally, droughts in the Southern YRSR were relatively more severe with longer drought duration; (3) in most areas of Togton River Basin and Dam River Basin, the severe ecological drought events occurred more frequently; (4) drought duration and severity in the YRSR were more susceptible to temperature when the temperature rise was above 1.0 °C. The average drought duration and severity increased by 20.7% and 32.6% with a temperature rise of 1 °C. Investigating and evaluating drought characteristics, causes, and drought index effectiveness provide essential information for balanced water resource allocation, utilization, and drought prevention. Understanding these spatial-temporal characteristics of drought and return period was useful for drought risk assessment and sustainable development of water resources.

Drought is one of the widespread climatic extreme events that is characterized by prolonged water deficit. Both social economic systems and ecological environment systems are vulnerable to drought [1,2,3]. According to statistics, about 8% of global vegetated land area has been affected by extreme drought, causing GPP (Gross Primary Productivity) reduction of 1.5 PgC (petagrams of carbon) during 2001–2011 [4]. It has been estimated that the cereal losses due to droughts reached up to 1820 million Mg in worldwide in the past four decades [5]. Globally, more than 11 million deaths have been caused by drought disasters and more than 2 billion people have been affected since the 1990s [6]. Furthermore, the drought frequency and intensity in the mid to late 21st century were projected to increase in the context of climate change [6,7], which would lead to more drought problems in the future. Thus, it is important to assess drought characteristics in the planning and management of water resources. In recent research, a plenty of drought indices have been developed to identify different types of droughts, such as meteorological drought, hydrologic drought, agricultural drought, and socioeconomic drought [8,9]. The most widely used drought indices include standardized precipitation index (SPI) [10,11], Standardized Precipitation-Evapotranspiration Index (SPEI) [12], Palmer drought severity index (PDSI) [13], soil moisture drought index [14], and so on. In addition, a run theory for identifying drought parameters concerning drought duration, magnitude, and intensity was proposed by Yevjevich [15]. This approach has been applied in several drought models and analyses [16,17,18,19,20,21].

However, most drought indices use a single hydro-meteorological variable (e.g., precipitation, runoff) or in combination with other meteorological elements, loosely connected to ground conditions [22]. Although many drought indices are proposed for meteorological drought, hydrologic drought, and agricultural drought, few have been used for terrestrial natural vegetation. The Vegetation Condition Index (VCI) [23] is a kind of remote-sensing index and can be used for vegetation drought assessment. It takes Normalized Difference Vegetation Index (NDVI) as input, which has advantages in describing the vegetation biomass in drought condition. However, it is difficult to be used for drought forecasting compared with the indices with hydro-meteorological variables as input. A natural vegetation drought index with simple hydro-meteorological variables is needed for drought modeling and forecasting. Weng et al. proposed a Generalized Drought Assessment Index (GDAI) for assessing drought events [24]. This index considers water supply and water demand of water resources systems. Similarly, a natural vegetation drought index can be built with water supply and water demand during the growth period. Additionally, these two major inputs can be evaluated by hydro-meteorological variables, making it easy to model and forecast vegetation drought.

The Yangtze River Source Region (YRSR), located in the center of the Tibet Plateau, is one of areas that are sensitive and vulnerable to climate. Since the mid-20th century, the hydro-meteorological variables and water resources showed considerable changes in the YRSR under global warming [25,26]. The temperature, precipitation, actual evapotranspiration, and runoff of the YRSR increased by 0.34 °C/decade, 11.4 mm/decade, 7.6 mm/decade, and 3.3 mm/decade, respectively, in the last 60 years [27]. Furthermore, the annual average temperature and precipitation are projected to increase by 2.2 °C and 9.8% in the next 30 years based on CMIP5 (Coupled Model Intercomparison Project Phase 5) Climate Models [28]. Both the water supply and water demand of natural vegetation (especially grass and forest) has been influenced in these warmer and wetter climate conditions. Investigating natural vegetation drought characteristics in the context of climate change is helpful for freshwater planning and management. Thus, we addressed the following key issues: (1) introduce a new natural vegetation drought index considering the difference between water supply and demand during the growth period; (2) reveal the spatial-temporal variation characterization of natural vegetation drought in the YRSR; (3) investigate the drought-prone regions by bivariate probability and return period.

Yangtze River Source Region (Latitude: 32°25′ E and 35°53′ E; Longitude: 89°43′ E to 97°19′ E), located in the western Tibetan plateau, covers an area of 159,065 km^{2} (Figure 1). The elevation ranges from 6456 m in the West to 3512 mm in the East, with an average of 4779 m. The average annual precipitation is approximately 327.4 mm, which is concentrated from June to September. The average annual temperature is about −5.5–3.0 °C from northeast to southwest, with an average of 2.9 °C [27,29]. The drought index is 3.67 in the YRSR, which means the climate is very dry. The land cover in the YRSR consists primarily of grasslands (63.8%) and unused land (29.4%). Water area and forest land are less dominant land cover types, accounting for 6.5% and 0.3% of total area of YRSR. Due to the cold and dry climatic conditions, the eco-hydrology process of the whole region is sensitive to climate change.

The daily discharge data measured from 1960 to 2016 at Zhimenda hydrological station were provided by Qinghai Hydrology Bureau. The daily meteorological data, including precipitation, highest temperature, lowest temperature, average temperature, sunshine hours, wind speed, relative humidity over the past 57 years (1960 to 2016) at 12 meteorological stations were obtained from National Meteorological Information Center (http://data.cma.cn/). The distribution of these selected hydro-meteorological stations was shown in Figure 1.

Land use/Land cover (LULC) data with a spatial resolution of 1 × 1 km for the years 1990, 2000, and 2010 was obtained from Resource and Environment Data Cloud Platform (RESDC) (http://www.resdc.cn/). This dataset is generated by remote sensing images with manual visual interpretation method. The quality has been controlled and integration has been checked by RESDC.

The digital soil map and properties of dominant soil types were collected from the China soil database (http://www.soil.csdb.cn/), which is provided by Institute of Soil Science, Chinese Academy of Sciences.

The DEM data with a resolution of 90 m was downloaded from the CGIARCSI GeoPortal (https://srtm.csi.cgiar.org/).

The SSDI describes the degree of deviation between dry and wet conditions based on the difference between water supply and demand [30]. The algorithm of this index is similar to SPEI, but it defines drought in a broad term, considering water resources shortage instead of climatological factors. Taking crop evapotranspiration and effective precipitation as input, the SSDI has been used to identify the agricultural drought [30]. In the same way, it can be used to quantify the natural vegetation drought conditions. For the vegetation drought estimation in YRSR, this study made several assumptions. Taking into account only perennial or seasonal vegetation and their corresponding phenological stages. As there are hardly any rainfed crops in the Yangtze River Source Region, rainfed crops were not included. According to the previous studies, the green water is important for the natural grass and forest vegetation [31]. Thus, the green water flow (GWF) was selected as the water supply in this study. While, the water demand can be represented by ecological water requirement (EWR) of vegetation. With these two major variables, the SSDI for ecological drought evaluation is briefly described as follows:

The water surplus or deficit (WD_{i}) is defined as Equation (1)
where, GWF_{i} and EWR_{i} refer to green water flow and natural vegetation water requirement at each station or subregion for the ith month respectively. The GWF can be quantified by hydrological model. In this study, the ArcSWAT 2012 (https://swat.tamu.edu/software/arcswat) (Blackland Research & Extension Center, Temple, TX, USA) is used for model setup and parameterization and GWF estimation [32,33]. In the previous study, we already built a SWAT (Soil and Water Assessment Tool) model for YRSR and simulated the flows of green water, as shown in Yuan et al. [28]. Vegetation type, climate, and soil moisture can all determine the vegetation’s EWR, which was calculated by quota valuation mode in this research.
where ET_{0} represents the potential evapotranspiration of vegetation (mm), which can be calculated by Penman–Monteith formula. Kc represents the plant water potential coefficient; Ks represents the soil moisture coefficient. The detailed mathematical procedure can be found in Zhao et al. [34].

$$W{D}_{i}=GW{F}_{i}-EW{R}_{i}$$

$$EWR={K}_{s}\times {K}_{c}\times E{T}_{0}$$

Accumulated WD value in a window of k-months (D^{k}) is calculated as Equation (3).
where, n is the analyzed month.

$$\begin{array}{c}{D}_{n}^{k}={\displaystyle {\displaystyle \sum}_{i=0}^{k-1}}\left(GW{F}_{n-i}-EW{R}_{n-i}\right)\hspace{1em}n\ge k\end{array}$$

According to SPEI calculation procedure, a three-parameter log-logistic distribution is selected to consider the D^{k} series:
where, F(x) is the probability distribution function (PDF) of D^{k}; α, β, γ are the scale, shape, and origin parameters, respectively, which can be estimated by L-moment method [35].

$$F\left(x\right)={\left[1+{\left(\frac{\alpha}{x-\gamma}\right)}^{\beta}\right]}^{-1}$$

Finally, the obtained values of F(x) are converted into the standardized values by Equation (5)
where, P = 1 − F(x); according to SPEI calculation procedure, the constants are C_{0} = 2.515517, C_{1} = 0.802853, C_{2} = 0.010328, d_{1} = 1.432788, d_{2} = 0.189269, and d_{3} = 0.001308. Table 1 describes the categorization of drought magnitude by SSDI.

$$\begin{array}{l}if\hspace{1em}P\le 0.5\\ W=\sqrt{-2ln\left(P\right)}\\ SSDI=W-\frac{{c}_{0}+{c}_{1}W+{c}_{2}{W}_{2}}{1-{d}_{1}W+{d}_{2}{W}_{2}+{d}_{3}{W}_{3}}\\ else\\ W=\sqrt{-2ln\left(1-P\right)}\\ SSDI=\frac{{c}_{0}+{c}_{1}W+{c}_{2}{W}_{2}}{1-{d}_{1}W+{d}_{2}{W}_{2}+{d}_{3}{W}_{3}}-W\end{array}$$

The ‘run theory’ was proposed for drought characteristics definition by Yevjevich [15], which has been widely used for drought event identification [9,24]. In this study, we suggest that a drought event starts when the value of the SSDI value reaches −1.0 or less and ends when the SSDI value is more than −1.0 (Figure 2). Then, the characteristics of each drought event can be described by the parameters as listed below:

- 1.
- Drought duration (DD) is defined as the time between the initiation and termination of a drought event, which is expressed in months in this study.$$D{D}_{n}=T{t}_{n}-T{i}_{n}$$
_{n}refers to drought duration for the nth drought event. Ti_{n}and Tt_{n}are initiation time and termination time, respectively.

The average duration can be calculated as follows:
where DD_{avg} refers average duration during a given period. DD_{n} refers to drought duration for the nth drought event. ND is the number of drought events during a given period.

$$D{D}_{avg}=\frac{{{\displaystyle \sum}}_{n=1}^{ND}D{D}_{n}}{ND}$$

- 2.
- Drought severity (DS) is the sum of SSDIs during the drought duration.$$D{S}_{n}={\displaystyle \sum}_{i=1}^{D{D}_{n}}SSD{I}_{i}$$
_{n}refers to drought severity for the nth drought event.

The average severity can be calculated as follows:
where DS_{avg} refers to average severity during a given period. DS_{n} refers to drought severity for the nth drought event. ND is the number of drought events during a given period.

$$D{S}_{avg}=\frac{{{\displaystyle \sum}}_{n=1}^{ND}D{S}_{n}}{ND}$$

- 3.
- Drought peak (DP) is the maximum absolute value of SSDIs of a drought event.$$D{P}_{n}=\mathrm{max}\left(\left|SSD{I}_{T{i}_{n}}\right|,\left|SSD{I}_{T{i}_{n}+1}\right|,\cdots ,\left|SSD{I}_{T{t}_{n}-1}\right|,\left|SSD{I}_{T{t}_{n}}\right|\right)$$
_{n}refers to drought peak for the nth drought event. Ti_{n}and Tt_{n}are initiation time and termination time, respectively.

- 4.
- Drought coverage area (DA) is the region affected by the drought, which is calculated as follows:$$D{A}_{n}=\frac{A\left(n\right)}{A}\times 100\%$$
_{n}refers to the ratio of average drought affected area. A(i) is the area of region experiencing drought conditions. A is the area of the study region. In this research, the parameter A represents the area of grassland and forest ecosystem in the YRSR.

To identify the interannual trends of drought parameters (e.g., SSDI values and coverage area), the linear regression method was adopted to eliminate the increase or decrease rate [36], which can be calculated as follows:

Linear regression method was utilized to analyze the interannual trends of drought parameters (e.g., SSDI values and coverage area) [36]. The equation is listed as follows:
where θ_{slope} is the linear slope of the time series variable, which can be used to characterize the increase or decrease rate during a given study period; n is the number of years (here n = 15); X_{i} is the drought parameter for the ith year (I = 1,2,…n).

$${\theta}_{Slope}=\frac{n\times {{\displaystyle \sum}}_{i=1}^{n}\left(i\times {X}_{i}\right)-{{\displaystyle \sum}}_{i=1}^{n}i{{\displaystyle \sum}}_{i=1}^{n}{X}_{i}}{n\times {{\displaystyle \sum}}_{i=1}^{n}{i}^{2}-{({{\displaystyle \sum}}_{i=1}^{n}i)}^{2}}$$

The Mann–Kendall test is widely used in the field of hydrology and climatology [37]. The Mann–Kendall statistic, Z_{MK}, can be calculated by the following equation:
where, Statistic S can be calculated through the following equation:
where, j > k; n is the number of years, X_{j} and X_{k} are the annual values in the years j and k, respectively. The equation sgn(X_{j} − X_{k}) can be calculated as follows:
where, Z_{MK} follows a standard normal distribution, if |Z| > Z_{1}_{−α∕2}, where α denotes the significance level, then the trend is significant.

$${Z}_{MK}=\left\{\begin{array}{c}\frac{S-1}{\sqrt{VAR\left(S\right)}}\hspace{1em}if\hspace{1em}S>0\\ 0\hspace{1em}if\hspace{1em}S=0\\ \frac{S+1}{\sqrt{VAR\left(S\right)}}\hspace{1em}if\hspace{1em}S<0\end{array}\right.$$

$$S={\displaystyle \sum}_{k=1}^{n-1}{\displaystyle \sum}_{j=k+1}^{n}sgn\left({X}_{j}-{X}_{k}\right)$$

$$sgn\left({X}_{j}-{X}_{k}\right)=\left\{\begin{array}{c}1\hspace{1em}if\hspace{1em}{X}_{j}-{X}_{k}>0\\ 0\hspace{1em}if\hspace{1em}{X}_{j}-{X}_{k}=0\\ -1\hspace{1em}if\hspace{1em}{X}_{j}-{X}_{k}<0\end{array}\right.$$

The results of the statistical analysis above can obtained by MATLAB (The MathWorks, Natick, MA, USA).

Previous research has recommended a larger set of marginal distributions for fitting drought variables [38,39,40]. In this study, the exponential, Gamma, Weibull, Generalized Pareto, Generalized Extreme Value distributions were selected as candidates to fit the DD, DS, and DP. The cumulative distribution functions are presented in Equations (16)–(20).
where Γ(⋅) is the gamma function. F(x) is the marginal distribution and x is the value of drought variable. α, β, λ are scale, shape, and location parameters, respectively, which can be estimated by maximum likelihood method.

$$\mathrm{Exponential}:F\left(x\right)=1-exp\left(-\frac{x-\beta}{\alpha}\right)$$

$$\mathrm{Gamma}:F\left(x\right)={{\displaystyle \int}}_{0}^{x}\frac{{x}^{\alpha -1}}{{\beta}^{\alpha}\Gamma \left(\alpha \right)}exp\left(-\frac{x}{\beta}\right)dx$$

$$\mathrm{Weibull}:F\left(x\right)=1-exp{\left(-\frac{x}{\alpha}\right)}^{\beta}$$

$$\mathrm{Generalized}\mathrm{Pareto}:F\left(x\right)=1-{\left[1-\beta \left(\frac{x-\lambda}{\alpha}\right)\right]}^{\frac{1}{\beta}}$$

$$\mathrm{Generalized}\mathrm{Extreme}\mathrm{Value}\begin{array}{c}F\left(x\right)=exp\left\{-{\left[1-k\left(\frac{x-\lambda}{\alpha}\right)\right]}^{\frac{1}{\beta}}\right\}\hspace{1em}\left(\beta \ne 0\right)\end{array}$$

The fitted distributions were compared against empirical distribution estimated by Equation (21) [41].
where P_{m} is the empirical cumulative probability of the mth value. m represents mth smallest value in the dataset arranged in ascending order. n is the total number.

$${P}_{m}=\frac{m-0.44}{n+0.12}$$

Bivariate models based on copula functions were employed for probabilistic modelling and return period estimation in this study. A copula function C can link marginal distributions to their bivariate probability distribution as follows:
where F_{XY}(x,y) is the bivariate probability distribution. F_{X}(·) and F_{Y}(·) are marginal distributions. X and Y are correlated variables.

$$P\left(X\le x,Y\le y\right)={F}_{XY}\left(x,y\right)=C\left({F}_{X}\left(x\right),{F}_{Y}\left(y\right)\right)$$

Among all types of copulas, the Archimedean copulas with only one parameter derived by Genest and MacKay [43] have been widely used for analyzing drought characteristics [33,44,45]. In this research, three families of 2-dimensional Archimedean copulas, Clayton, Frank and Gumbel, were used to construct bivariate distributions for DD vs. DS, DD vs. DP and DS vs. DP. Expressions of the mentioned copulas are presented as follows:
where C(u,v) is the copula bivariate distribution function for marginal cumulative probabilities u and v. θ is the parameter which can be estimated by maximum likelihood method.

$$\mathrm{Clayton}:\begin{array}{c}C\left(u,v\right)=max\left({\left({u}^{-\theta}+{v}^{-\theta}-1\right)}^{-\frac{1}{\theta}},0\right)\hspace{1em}\theta \in \left(0,\infty \right)\end{array}$$

$$\mathrm{Frank}:\begin{array}{c}C\left(u,v\right)=-\frac{1}{\theta}ln\left(1+\frac{\left({e}^{-\theta u}-1\right)\left({e}^{-\theta v}-1\right)}{{e}^{-\theta}-1}\right)\hspace{1em}\theta \in (-\infty ,0)\cap (0,+\infty )\end{array}$$

$$\mathrm{Gumbel}:\begin{array}{c}C\left(u,v\right)=exp\left(-{\left({\left(-lnu\right)}^{\theta}+{\left(-lnv\right)}^{\theta}\right)}^{\frac{1}{\theta}}\right)\hspace{1em}\theta \in \left.1,\infty \right)\end{array}$$

The empirical copula function is introduced to estimate goodness of fit [46], which is written as:
where C_{emp}(u,v) is the empirical copula function. R_{i} and S_{i} are the ranks of the ith observed data. I(A) is the indicator function of event A which is presented as Equation (27).

$${C}_{\mathrm{emp}}\left(u,v\right)=\frac{1}{n}{{\displaystyle \sum}}_{i=1}^{n}I\left(\frac{{R}_{i}}{n+1}\le u,\frac{{S}_{i}}{n+1}\le v\right)$$

$$I\left(A\right)=\left\{\begin{array}{l}1\hspace{1em}A\hspace{1em}is\hspace{1em}True\\ 0\hspace{1em}A\hspace{1em}is\hspace{1em}False\end{array}\right.$$

The fitted distributions were compared against empirical distribution and then the K-S tests and SED were used to test goodness of fit and select optimal copula function.

The bivariate probability and return period for DD vs. DS, DD vs. DP and DS vs. DP can be defined by two cases, namely Case “∩” and Case “∪”. Taking pair of DD vs. DS as an example, the two cases can be described as follow:

(1) Case “∩”: (DD > d) and (DS > s)

(2) Case “∪”: (DD > d) or (DS > s)

The bivariate probabilities in these two cases can be estimated by Equations (28) and (29) [38].
where, P_{D}_{∩P} and P_{D}_{∪P} are bivariate probabilities for Case “∩” and Case “∪”, respectively. F_{D}(·) and F_{S}(·) are the marginal distributions for drought duration and severity, respectively. C(·) is the copula bivariate distribution for drought duration and severity.

$${P}_{D\cap S}=P\left(\begin{array}{ccc}DD\ge d& and& DS\ge s\end{array}\right)=1-{F}_{D}\left(d\right)-{F}_{S}\left(s\right)+C\left({F}_{D}\left(d\right),{F}_{S}\left(s\right)\right)$$

$${P}_{D\cup S}=P\left(\begin{array}{ccc}DD\ge d& or& DS\ge s\end{array}\right)=1-C\left({F}_{D}\left(d\right),{F}_{S}\left(s\right)\right)$$

With the copula-based distribution, the bivariate return period in the mentioned two cases is denoted as [18,38]:
where, T_{D}_{∩P} and T_{D}_{∪P} are bivariate probabilities for Case “∩” and Case “∪”, respectively. ζ =NY/ND, NY refers to the period of SSDI time series (years), and ND is the number of drought events in NY years [39].

$${T}_{D\cap S}=\frac{\zeta}{P\left(\begin{array}{ccc}DD\ge d& and& DS\ge s\end{array}\right)}=\frac{\zeta}{1-{F}_{D}\left(d\right)-{F}_{S}\left(s\right)+C\left({F}_{D}\left(d\right),{F}_{S}\left(s\right)\right)}$$

$${T}_{D\cup S}=\frac{\zeta}{P\left(\begin{array}{ccc}DD\ge d& or& DS\ge s\end{array}\right)}=\frac{\zeta}{1-C\left({F}_{D}\left(d\right),{F}_{S}\left(s\right)\right)}$$

One of the major motivations of this research is to develop a reliable and reasonable method that can be used to assess drought in forest and grass ecosystems objectively. The Normalized Difference Vegetation Index (NDVI) values can depict the vegetation biomass, which is an ancillary reflection of drought conditions [47]. In view of this dynamic characteristic, we compared the SSDI with NDVI to discover whether the new index constructed in this study can capture the regional drought. Figure 3 illustrates the relationship between SSDI/SPEI and NDVI during 2000–2015. It can be found that NDVI has much higher Pearson’s correlation coefficient with SSDI (r_{NDVI}_{, SSDI}) than does SPEI (r_{NDVI}_{, SPEI}). For example, the sub-basins with r_{NDVI}_{, SSDI} ≥ 0.4 account for 63.4% of the total, while the sub-basins with r_{NDVI}_{, SPEI} ≥ 0.4 account for only 25.6% (Figure 3). This indicates that the SSDI is superior to the SPEI for ecological drought assessment in the YRSR because it can capture the NDVI variation.

The interannual variability of the SSDI at different timescales is shown in Figure 4. We found that the SSDI series at different regions and the whole YRSR exhibit a wetting trend during 1960–2016. Taking the SSDI-12m of December as an example, the SSDI value increased by 0.0154 ± 0.0063 per year in the study period. It also can be seen that persistent droughts occurred frequently in 1977–1979 and 1984–1997, and periods of late-1970s, mid-1980s, and mid-1990s had witnessed the most severe drought episodes. A clear reversal of dry/wet condition in 2006 can be detected in Figure 3, which indicates that drought has been gradually relieved since mid-2000s. This change has happened largely because of the effects of climate change. A couple of studies have reported that the precipitation and actual evapotranspiration in the YRSR have been increased in the last decade [27,28,48]. In other words, the water supply for grassland and forest ecosystem would increase in this context, causing drought events to become less frequent in many regions of the YRSR.

To investigate the temporal evolution of dry/wet condition in YRSR, the linear regression analysis method and nonparametric Mann–Kendall test method were applied to analyze the trends of SSDI from 1960 to 2016. The dry/wet conditions of year, spring, summer, autumn, and winter are represented by SSDI-12m of December, SSDI-3m of May, SSDI-3m of August, SSDI-3m of November, and SSDI-3m of February (Figure 5). For more than four-fifths of the sub-basins across YRSR, the SSDI-12m of December presented increasing trends from 1960 to 2016 which indicate that most parts of YRSR have experienced wetting trends in the study period. In addition, the significant increasing trends of SSDI-12m (α = 0.05) have been detected in nearly half of the sub-basins, especially in northern YRSR (Figure 5a). The SSDI-3m of May, SSDI-3m of August, SSDI-3m of November showed a similar trend as the SSDI-12m, which illustrated that the dry conditions in spring, summer and autumn had been alleviated (Figure 5b–d). However, the change of SSDI values in winter was much different from that in other seasons. The dry condition in winter appeared to be more severe in Togton River Basin and Dam River Basin. However, the trend did not pass the significance test through the Mann–Kendall method (Figure 5e).

Percentages of drought affected areas for seasonal and yearly scales were calculated with Equation (11) and the results are presented in Figure 6. It can be easily identified that the widespread drought events occurred frequently in the 1990s. The average percentages of drought affected areas (SSDI-12m < −1) were found to be 26.2% for this period, more than 1.6 times the average value for the 1960–2016 period. The notable severe ecological drought years were 1979, 1995, 1994, and 1960, with more than 75% of areas being affected (Figure 6a). The drought affected areas in spring, summer, and autumn have decreased since the early 2000s, but in winter have increased, especially in the most recent 10 years (Figure 6b–e).

According to the drought event identification and characterization based on Run theory and SSDI at 6-month timescale, the spatial patterns of drought duration, severity, peak, maximum severity, and drought count in 1960–2016 are illustrated in Figure 7. The spatial distribution of average drought duration showed that the Southern YRSR (e.g., Dam River Basin, Nyaqa River Basin) suffered from droughts of a longer duration compared with the northern YRSR (e.g., Qumar River Basin, Beilu River Basin). However, the peak of SSDI in the northern YRSR was generally higher than that in Southern YRSR. The spatial distribution in average drought severity agreed well with that in drought duration, namely, high in the south and low in the north. In addition, the spatial distribution of maximum severity and number of drought events based on SSDI at 6-month timescale was computed in this research. It is clear that the higher maximum severity (S_{max}) recorded by SSDI was within the northern YRSR, with S_{max} more than 25. Average drought count was significantly higher in Togton River Basin and Dam River Basin, where precipitation is less and temperature is lower.

The time scale of SSDI is flexible, e.g., a 3-month SSDI can be used as a short-term drought index while a 12-month SSDI can be used as an intermediate-term one. Figure 8 compared the SSDI value for the time scales of 3, 6, and 12 months. It can be found that the time series of SSDI-12M are more ‘smoothed’ than that of SSDI-3M and SSDI-6M. That is, the dry/wet conditions identified by smaller time scale SSDI took place more alternately. The scatter plots for SSDI-3M and SSDI-6M/12M are shown in Figure 9. Drought grades identified by SSDI-6M/12M were generally lower than that identified by SSDI-3M. The difference between SSDI-3M and SSDI-6M was relatively smaller with the correlation coefficient nearly to 0.85. However, the difference between SSDI-3M and SSDI-12M was obvious. The characteristics of drought events derived from multiscale SSDI were compared in Figure 10. It can be concluded that the drought characteristics are affected by the time scales. The drought event appears to be longer and intensified but low frequency at higher timescales. Taking the entire YRSR as an example, the median drought duration and severity were 2.0 months and 2.4, respectively, at 3-month scale, while 8.0 months and 10.5, respectively, at 12-month scale. These results show that the choice of time scales to calculate SSDI value can introduce uncertainties in drought characteristics evaluated.

Previous studies have shown that the temperature has been increasing significantly in the YRSR [49,50,51]. In the context of warming climate, both the green water flow and ecological water requirement would increase, which may lead to complex changes in the balance between water supply and demand. Additionally, then drought characteristics are likely to be changed accordingly. To investigate the sensitivity of drought characteristics to the temperature in the YRSR, a series of hypothetical temperature increase scenarios, changing from 0.5 to 2.5 °C in increments of 0.5 °C, were simulated and analyzed. Generally, the drought duration and severity in the whole YRSR would increase along with the increase of temperature (Figure 11). In the hypothetical scenarios with ΔT ≤ 1 °C, the average drought duration and severity were less susceptible to changes in temperature. The resulting slope indicates that a 1 °C increase in temperature corresponds to increases in average drought duration and severity by 4.2% and 6.9%, respectively (dashed red lines in Figure 11). It is clear that the drought characteristics were more susceptible to the temperature increase in the hypothetical scenarios with ΔT > 1 °C. The results showed increases of 20.7% and 32.6% for drought duration and severity respectively when temperature was increased by 1 °C (dashed red lines in Figure 12). According to the previous study, the annual average temperature in the YRSR was projected to increase by 2.2 °C in the next 30 years based on CMIP5 Climate Models [28]. In addition, that means the ecological droughts with longer duration and larger severity may happen more frequently in the future. However, this trend needs to be further explored because the precipitation was also projected to increase. In addition, the temperature increase would also have an effect on growth period of grass or forest in the YRSR [51,52,53]. This change in phenology would alter the water demand process and drought characteristics indirectly, which should be further investigated.

The exponential, Gamma, Weibull, Generalized Pareto, and Generalized Extreme Value distributions were chosen as appropriate marginal distributions for the DD, DS, and DP in each subregion and the entire YRSR. The best fit distributions for each drought variable were selected. The theoretical probability and empirical probability of selected distribution for drought variables were plotted in Figure 12. It can be seen from Figure 12 that the associations of theoretical and empirical probabilities are close to the 1:1 line, which proved the reliability of marginal distribution selection.

With the optimal marginal distribution of drought variables chosen above, Copula function for modeling bivariate distribution could be structured. Three kinds of Archimedean copulas, namely Clayton, Frank and Gumbel, were compared in this research and the estimated parameters for best-fitted copulas are listed in Table 2. In case of DD vs. DS, the Frank Copula in all six subregions and the entire YRSR was the best function with the lowest SED value. Details of the optimal copulas for DD vs. DP and DS vs. DP in subregions can be found in Table 2. The comparison plot between theoretical joint cumulative probability and empirical joint cumulative probability is shown in Figure 13. It could be observed that the theoretical value is almost equal to the empirical one. This result proved that the selected copulas and estimated parameters are reliable. It should be noted that Frank Copula appears to be proper and robust in modeling bivariate distribution of drought variables in many of the subregions of YRSR.

Using the selected optimal copulas in I–V subregions and the entire YRSR, the bivariate probability of DD vs. DS, DD vs. DP and DS vs. DP in both Case “∩” and Case “∪” can be estimated and the results are presented in Figure 14. For example, when the DD and DS exceed 3 months and 3, respectively (above moderate grade), the P_{D}_{∩S} was 0.459, 0.434, 0.453, 0.484, 0.508 and the P_{D}_{∪S} was 0.590, 0.580, 0.577, 0.613, 0.645 in subregions I–V, respectively. For the drought event with DD > 3 and DS > 4.5(above severe grade), the P_{D}_{∩S} in the five subregions was 0.435, 0.411, 0.422, 0.461, 0.489 and P_{D}_{∪S} was 0.478, 0.465, 0.467, 0.504, 0.537. It could be found that the joint probabilities are different among the five subregions, with significantly larger P_{D}_{∩S} and P_{D}_{∪S} in Dam River Basin and Qumar River Basin, implying that the drought risk is higher in these regions.

The bivariate return periods (T) of DD vs. DS, DD vs. DP and DS vs. DP for subregions I–V and the YRSR were estimated according to Equations (30) and (31), and results are illustrated in Figure 15. It can be found that various pairs of drought variables can result in the same return period. With the 3D-surface figure, the associated return period in Case ∩ and Case ∪ for a specific drought variable can be extracted. For example, DD = 6 months and DS = 6 result in T_{D}_{∩S} = 11.5 year and T_{D}_{∪S} = 6.2 years for the YRSR. The two kinds of return period are compared, and the “∪” return periods are shorter than the “∩” return periods in all the sub-basins and the entire YRSR. Further, the univariate and bivariate drought return periods for different variation patterns were compared (Table 3). Apparently, the univariate return period is longer than “∪” return periods, while shorter than the “∩” return periods. With the DD, DP and DS increasing, both the “∪” and “∩” return periods increased. In addition, the difference between “∪” and “∩” return periods would also increase as the values of drought variables increase. Taking the entire YRSR as an example, the DD, DS and DP were about 11.1, 17.2 and 1.8, respectively, in a situation whereby the univariate return period was 50 years. For this particular event, the T_{D}_{∪S}, T_{D}_{∪S} and T_{P}_{∪S} were about 34.1, 28.9 and 29.6 years while the T_{D}_{∩S}, T_{D}_{∩S} and T_{P}_{∩S} were about 93.3, 186.6 and 161.5 years. The results are important for drought risk assessment and would be helpful in planning and management of water resources systems under severe or extreme drought scenarios.

The average DD, DS and DP for the YRSR are 3.46 months, 5.37 and 1.54, respectively. For these particular drought characteristics, the bivariate probabilities of the 119 sub-basins divided for spatial analysis [28] were estimated and mapped in Figure 16. A noticeable spatial variability of P_{D}_{∪S} and P_{D}_{∩S} could be observed obviously. The bivariate probabilities of DD vs. DS increased from the north to the south, implying that the droughts in southern YRSR are relatively more severe and have longer duration. In the case of DD > 3.46 and DS > 5.37, the bivariate probabilities of P_{D}_{∪S} and P_{D}_{∩S} were more than 0.40 and more than 0.35 in the Dam River Basin. While in the Qumar River Basin, P_{D}_{∪S} and P_{D}_{∩S} were evaluated as 0.30–0.45 and 0.25–0.40, respectively. Comparing with P_{D}_{∪S} and P_{D}_{∩S}, the spatial difference of P_{D}_{∪P}, P_{D}_{∩P}, P_{S}_{∪P}, and P_{S}_{∩P} among YRSR is not obvious. P_{D}_{∪P} and P_{S}_{∪P} in nearly all of the sub-basins in YRSR were more than 0.45. The P_{D}_{∩P} and P_{S}_{∩P} in most parts of YRSR ranged from 0.30 to 0.33.

The spatial distribution of bivariate return periods is similar with that of probabilities (Figure 17). The short return periods are often associated with high probabilities and dispersed in southern YRSR, especially for DD vs. DS. In most areas of Togton River Basin and Dam River Basin, the T_{D}_{∪S} and T_{D}_{∩S} remained in the variation range from 2.33–2.49 years and 2.67–2.85 years. While, those two kinds of return periods in Qumar River Basin were around 2.6 and 3.0 years, respectively. These results suggest that severe ecological drought events are more likely to occur in the Togton River Basin and Dam River Basin.

This study introduced a new natural vegetation drought index, standardized supply-demand water index (SSDI), and identified spatial-temporal variability of natural vegetation drought characteristics drought events. Then the drought-prone regions were investigated. The primary conclusions are as follows:

- The time-series of SSDI and Standardized Precipitation and Evapotranspiration Drought Index (SPEI) with Normalized Difference Vegetation Index (NDVI) were compared in this study. There exists a higher correlation between constructed SSDI and NDVI. This result indicated that the constructed SSDI was reliable and can reflect the comprehensive characteristics of the ecological drought information more easily and effectively.
- The YRSR had witnessed the most severe drought episodes in the periods of late-1970s, mid-1980s and mid-1990s, but the SSDI showed a wetting trend since the mid-2000s, mainly because of a warmer and wetter climate in the most recent 10 years. However, the climate change has different effects on the dry condition at seasonal scales. The drought affected areas in spring, summer and autumn have decreased since 2000 while this area in winter has increased. The drought duration and severity showed a spatial variation among different regions in the YRSR. Generally, droughts in the Southern YRSR were relatively more severe with longer drought duration, implying that the Southern YRSR was an area that had been facing challenging drought conditions. The average drought duration and severity in the YRSR would be less susceptible to changes in temperature when the increase temperature was above 1.0 °C. However, the characteristics would be more susceptible to temperature in the YRSR when the increase temperature were above 1.0 °C. The average drought duration and severity is shown to increase by 20.7% and 32.6% with a 1 °C increase in temperature for the hypothetical scenarios with ΔT > 1 °C.
- The return periods of five sub-basins and the entire YRSR for case ‘‘∩” were longer than those in case ‘‘∪” and their spatial trends are highly consistent. High return periods were found in Qumar River Basin. While, low return periods were found in most areas of Togton River Basin and Dam River Basin, implying that severe ecological drought events occurred more frequently.

The study investigated the spatial-temporal variation characterization of ecological drought in the Yangtze River Source Region, which provides a comprehensive understanding of regional drought events in the YRSR and highlights the bivariate drought probabilities and return periods. However, two aspects should be improved in further study. First, the choice of time scales to calculate SSDI value can introduce uncertainties in drought characteristics evaluated. Therefore, future study should compare multiscale SSDI indices and select the optimal time scale for drought assessment in the YRSR. Second, temperature is one of the most crucial and direct driving factors of the changes in growth period. Additionally, the changes in growth period for grass and forest could affect the process of ecological water requirement. This dynamic change should be quantified in the algorithm of SSDI. Furthermore, the current paper selected the YRSR as a case study and the methodological development of SSDI and drought identification can be used in other catchments or regions as well.

Z.Y. worked on forming ideas of this paper; J.Y. and Z.Y. worked together in calculating and writing of this manuscript; T.L. was responsible for the preliminary data processing; Z.Y. provided supervision during the whole process. All authors read and approved the final manuscript.

This research was funded by [National Natural Science Foundation of China] grant number [41890821, 52079008, 51909080]; [National Key Research and Development Project] grant number [2017YFC1502404]; [National Public Research Institutes for Basic R&D Operating Expenses Special Project] grant number [CKSF2019292/SH+SZ, CKSF2017061/SZ]; [Young and middle-aged talent project of Hubei Province Education Department] grant number [Q20181001]; [Open Research Fund of Changjiang River Scientific Research Institute, Changjiang Water Resources Commission of the Ministry of Water Resources of China] grant number [CKWV2018493/KY].

Not applicable.

Not applicable.

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

- Touma, D.; Ashfaq, M.; Nayak, M.; Kao, S.-C.; Diffenbaugh, N. A multi-model and multi-index evaluation of drought characteristics in the 21st century. J. Hydrol.
**2015**, 526, 196–207. [Google Scholar] [CrossRef] - Vicca, S.; Balzarolo, M.; Filella, I.; Granier, A.; Herbst, M.; Knohl, A.; Longdoz, B.; Mund, M.; Nagy, Z.; Pintér, K.; et al. Remotely-sensed detection of effects of extreme droughts on gross primary production. Sci. Rep.
**2016**, 6, 28269. [Google Scholar] [CrossRef] - Schwalm, C.R.; Anderegg, W.R.L.; Michalak, A.M.; Fisher, J.B.; Biondi, F.; Koch, G.; Litvak, M.; Ogle, K.; Shaw, J.D.; Wolf, A.; et al. Global patterns of drought recovery. Nature
**2017**, 548, 202–205. [Google Scholar] [CrossRef] - Du, L.; Mikle, N.; Zou, Z.; Huang, Y.; Shi, Z.; Jiang, L.; McCarthy, H.R.; Lou, Y. Global patterns of extreme drought-induced loss in land primary production: Identifying ecological extremes from rain-use efficiency. Sci. Total Environ.
**2018**, 628–629, 611–620. [Google Scholar] [CrossRef] - Lesk, C.; Rowhani, P.; Ramankutty, N. Influence of extreme weather disasters on global crop production. Nature
**2016**, 529, 84–87. [Google Scholar] [CrossRef] - United Nations International Strategy for Disaster Reduction Secretariat. Global Assessment Report on Disaster Risk Reduction. In Risk and Poverty in a Changing Climate, Invest Today for a Safer Tomorrow; UN, International Strategy for Disaster Reduction: Geneva, Switzerland, 2009. [Google Scholar]
- Stocker, T.F.; Qin, D.; Plattner, G.K.; Tignor, M.M.B.; Allen, S.K.; Boschung, J.; Nauels, A.; Xia, Y.; Bex, V.; Midgley, P.M. Climate Change 2013: The physical Science Basis. In Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2013. [Google Scholar]
- World Meteorological Organization. Drought and Agriculture. World Meteorological Organization; Technical Note 138; WMO Publication: Geneva, Switzerland, 1975. [Google Scholar]
- Mishra, A.K.; Singh, V.P. A review of drought concepts. J. Hydrol.
**2010**, 391, 202–216. [Google Scholar] [CrossRef] - McKee, T.B.; Doesken, N.J.; Kleist, J. The Relationship of Drought Frequency and Duration to Time Scales. In Proceedings of the 8th Conference on Applied Climatology, Anaheim, CA, USA, 7–22 January 1993. [Google Scholar]
- McKee, T.B.; Doesken, N.J.; Kleist, J. Drought Monitoring with Multiple Time Scales. In Proceedings of the 9th Conference on Applied Climatology, Dallas, TX, USA, 15–20 January 1995. [Google Scholar]
- Serrano, S.M.V.; Beguería, S.; Moreno, J.I.L. A multiscalar drought index sensitive to global warming: The standardized precipitation evapotranspiration index. J. Clim.
**2010**, 23, 1696–1718. [Google Scholar] [CrossRef] - Palmer, W.C. Meteorologic Drought. In Research Paper No. 45; US Department of Commerce Weather Bureau: Washington, DC, USA, 1965. [Google Scholar]
- Hollinger, S.E.; Isard, S.A.; Welford, M.R. A New Soil Moisture Drought Index for Predicting Crop Yields. In Applied Climatology; American Meteor Society: Anaheim, CA, USA, 1993. [Google Scholar]
- Yevjevich, V. An Objective Approach to Definitions and Investigations of Continental Hydrologic Droughts. In Hydrology Papers 23; Colorado State University Publication: Fort Collins, CO, USA, 1967. [Google Scholar]
- Loaiciga, H.A.; Leipnik, R.B. Stochastic renewal model of low-flow streamflow sequences. Stoch. Hydrol. Hydraul.
**1996**, 10, 65–85. [Google Scholar] [CrossRef] - Chung, C.H.; Salas, J.D. Drought occurrences probabilities and risks of dependent hydrologic processes. J. Hydrol. Eng.
**2000**, 5, 259–268. [Google Scholar] [CrossRef] - Salas, J.D.; Fu, C.; Cancelliere, A.; Dustin, D.; Bode, D.; Pineda, A.; Vincent, E. Characterizing the severity and risk of drought in the Poudre River. Colorado. J. Water Resour. Plan. Manag.
**2005**, 131, 383–393. [Google Scholar] [CrossRef] - Mishra, A.K.; Desai, V.R.; Singh, V.P. Drought forecasting using a hybrid stochastic and neural network model. J. Hydrol. Eng.
**2007**, 12, 626–638. [Google Scholar] [CrossRef] - Serinaldi, F.; Bonaccorso, B.; Cancelliere, A.; Grimaldi, S. Probabilistic characterization of drought properties through copulas. Phys. Chem. Earth
**2009**, 34, 596–605. [Google Scholar] [CrossRef] - Jiang, R.; Xie, J.; He, H.; Luo, J.; Zhu, J. Use of four drought indices for evaluating drought characteristics under climate change in Shaanxi, China: 1951-2012. Nat. Hazards
**2015**, 75, 2885–2903. [Google Scholar] [CrossRef] - Zargar, A.; Sadiq, R.; Naser, B.; Khan, F.I. A review of drought indices. Environ. Rev.
**2011**, 19, 333–349. [Google Scholar] [CrossRef] - Kogan, F.N. Remote sensing of weather impacts on vegetation in non-homogeneous areas. Int. J. Remote Sens.
**1990**, 11, 1405–1419. [Google Scholar] [CrossRef] - Weng, B.; Yan, D.; Wang, H.; Liu, J.; Yang, Z.; Qin, T.; Yin, J. Drought assessment in the Dongliao River basin: Traditional approaches vs. generalized drought assessment index based on water resources systems. Nat. Hazards Earth Syst. Sci.
**2015**, 15, 1889–1906. [Google Scholar] [CrossRef] - Yao, Z.; Liu, Z.; Huang, H.; Liu, G.; Wu, S. Statistical estimation of the impacts of glaciers and climate change on river runoff in the headwaters of the Yangtze River. Quat. Int.
**2014**, 336, 89–97. [Google Scholar] [CrossRef] - Wang, R.; Yao, Z.; Liu, Z.; Wu, S.; Jiang, L.; Wang, L. Snow cover variability and snowmelt in a high-altitude ungauged catchment. Hydrol. Process.
**2015**, 29, 3665–3676. [Google Scholar] [CrossRef] - Du, Y.; Berndtsson, R.; An, D.; Zhang, L.; Hao, Z.; Yuan, F. Hydrologic Response of Climate Change in the Source Region of the Yangtze River, Based on Water Balance Analysis. Water
**2017**, 9, 115. [Google Scholar] [CrossRef] - Yuan, Z.; Xu, J.; Wang, Y. Historical and future changes of blue water and green water resources in the Yangtze River source region, China. Theor. Appl. Climatol.
**2019**, 138, 1035–1047. [Google Scholar] [CrossRef] - Chen, J. Water cycle mechanism in the source region of Yangtze River. J. Yangtze River Sci. Res. Inst.
**2013**, 30, 1–5. [Google Scholar] - Yuan, Z.; Xu, J.; Chen, J.; Huo, J.; Yu, Y.; Locher, P.; Xu, B. Drought assessment and projection under climate change: A case study in the middle and lower Jinsha River Basin. Adv. Meteorol.
**2017**, 5757238. [Google Scholar] [CrossRef] - Schuol, J.; Abbaspour, K.C.; Yang, H. Modeling blue and green water availability in Africa. Water Resour. Res.
**2008**, 44, 212–221. [Google Scholar] [CrossRef] - Zhang, W.; Zha, X.; Li, J.; Liang, W.; Ma, Y.; Fan, D.; Li, S. Spatialtemporal change of blue water and green water resources in the headwater of Yellow River Basin, China. Water Resour. Manag.
**2014**, 28, 4715–4732. [Google Scholar] [CrossRef] - Abbaspour, K.C.; Rouholahnejad, E.; Vaghefi, S.; Srinivasan, R.; Yang, H.; Kløve, B. A continental-scale hydrology and water quality model for Europe: Calibration and uncertainty of a high-resolution large-scale SWAT model. J. Hydrol.
**2015**, 524, 733–752. [Google Scholar] [CrossRef] - Zhao, F.; Li, H.; Li, C.; Cai, Y.; Wang, X.; Liu, Q. Analyzing the influence of landscape pattern change on ecological water requirements in an arid/semiarid region of China. J. Hydrol.
**2019**, 578, 124098. [Google Scholar] [CrossRef] - Ahmad, M.I.; Sinclair, C.D.; Werritty, A. Log-logistic flood frequency analysis. J. Hydrol.
**1988**, 98, 205–224. [Google Scholar] [CrossRef] - Yin, J.; Yan, D.; Yang, Z.; Yuan, Z.; Yuan, Y.; Wang, H.; Shi, X. Research on Historical and Future Spatial-Temporal Variability of Precipitation in China. Adv. Meteorol.
**2016**, 2016, 9137201. [Google Scholar] [CrossRef] - Tabari, H.; Talaee, P.H.; Ezani, A.; Some’e, B. Shift changes andmonotonic trends in autocorrelated temperature series over Iran. Theor. Appl. Climatol.
**2012**, 109, 95–108. [Google Scholar] [CrossRef] - Shiau, J.T. Fitting drought duration and severity with two-dimensional copulas. Water Resour. Manag.
**2006**, 20, 795–815. [Google Scholar] [CrossRef] - Ayantobo, O.O.; Li, Y.; Song, S.; Javed, T.; Yao, N. Probabilistic modelling of drought events in China via 2-dimensional joint copula. J. Hydrol.
**2018**, 559, 373–391. [Google Scholar] [CrossRef] - Guo, Y.; Huang, S.; Huang, Q.; Wang, H.; Wang, L.; Fang, W. Copulas-based bivariate socioeconomic drought dynamic risk assessment in a changing environment. J. Hydrol.
**2019**, 575, 1052–1064. [Google Scholar] [CrossRef] - Gringorten, I.I. A plotting rule for extreme probability paper. J. Geophys. Res.
**1963**, 68, 813–814. [Google Scholar] [CrossRef] - Weng, B.; Zhang, P.; Li, S. Drought risk assessment in China with different spatial scales. Arab. J. Geosci.
**2015**, 8, 10193–10202. [Google Scholar] [CrossRef] - Genest, C.; MacKay, J. The joy of copulas: Bivariate distributions with uniform marginals. Am. Stat.
**1986**, 40, 280–283. [Google Scholar] - Wang, R.; Zhao, C.; Zhang, J.; Guo, E.; Li, D.; Alu, S.; Ha, S.; Dong, Z. Bivariate copula function-based spatial-temporal characteristics analysis of drought in Anhui Province, China. Meteorol. Atmos. Phys.
**2019**, 131, 1341–1355. [Google Scholar] [CrossRef] - Ribeiro, A.F.; Russo, A.; Gouveia, C.M.; Páscoa, P. Copula-based agricultural drought risk of rainfed cropping systems. Agric. Water Manag.
**2019**, 223, 105689. [Google Scholar] [CrossRef] - Nelsen, R.B. An Introduction to Copulas. In Statistics; Springer: New York, NY, USA, 2006. [Google Scholar]
- McVicar, T.R.; Jupp, D.L.B. The current and potential operational uses of remote sensing to aid decisions on drought exceptional circumstances in Australia: A review. Agric. Syst.
**1998**, 57, 399–468. [Google Scholar] [CrossRef] - Liu, J.; Chen, J.; Xu, J.; Lin, Y.; Yuan, Z.; Zhou, M. Attribution of Runoff Variation in the Headwaters of the Yangtze River Based on the Budyko Hypothesis. Int. J. Environ. Res. Public Health
**2019**, 16, 2506. [Google Scholar] [CrossRef] - Kang, S.; Zhang, Y.; Qin, D.; Ren, J.; Zhang, Q.; Grigholm, B.; Mayewski, P.A. Recent temperature increase recorded in an ice core in the source region of Yangtze River. Chin. Sci. Bull.
**2007**, 52, 825–831. [Google Scholar] [CrossRef] - Chongyi, E.; Hu, H.; Xie, H.; Sun, Y. Temperature change and its elevation dependency in the source region of the Yangtze River and Yellow River. J. Geogr. Geol.
**2014**, 6, 124–131. [Google Scholar] - Piao, S.; Fang, J.; Zhou, L.; Ciais, P.; Zhu, B. Variations in satellite-derived phenology in China’s temperate vegetation. Glob. Chang. Biol.
**2006**, 12, 672–685. [Google Scholar] [CrossRef] - Ding, M.; Zhang, Y.; Sun, X.; Liu, L.; Wang, Z.; Bai, W. Spatialtemporal variation in alpine grassland phenology in the Qinghai-Tibetan Plateau from 1999 to 2009. Chin. Sci. Bull.
**2013**, 58, 396–405. [Google Scholar] [CrossRef] - Wang, S.; Zhang, B.; Yang, Q.; Chen, G.; Yang, B.; Lu, L.; Shen, M.; Peng, Y. Responses of net primary productivity to phenological dynamics in the Tibetan Plateau, China. Agric. For. Meteorol.
**2017**, 232, 235–246. [Google Scholar] [CrossRef]

Drought Grade | Range of SSDI Value |
---|---|

Near normal | 1.00 > SSDI > −1.00 |

Moderate drought | −1.00 ≥ SSDI > −1.50 |

Severe drought | −1.50 ≥ SSDI > −2.00 |

Extreme drought | −2.00 ≥ SSDI |

Subregion | DD vs. DS | DD vs. DP | DS vs. DP | |||
---|---|---|---|---|---|---|

Copula | Parameter (θ) | Copula | Parameter (θ) | Copula | Parameter (θ) | |

I | Frank | 33.094 | Frank | 8.192 | Frank | 10.202 |

II | Frank | 25.537 | Frank | 9.184 | Frank | 11.220 |

III | Frank | 31.090 | Frank | 7.205 | Frank | 8.897 |

IV | Frank | 32.708 | Frank | 9.804 | Frank | 11.906 |

V | Frank | 18.411 | Clayton | 2.042 | Clayton | 3.509 |

YRSR | Frank | 27.058 | Frank | 9.024 | Frank | 10.889 |

Region | T | DD | DS | DP | DD vs. DS Return Period | DD vs. DP Return Period | DS vs. DP Return Period | |||
---|---|---|---|---|---|---|---|---|---|---|

Case∪ | Case ∩ | Case∪ | Case ∩ | Case∪ | Case ∩ | |||||

Subregion I | 5 | 3.6 | 5.3 | 0.2 | 4.7 | 5.3 | 4.1 | 6.3 | 4.3 | 6.0 |

10 | 6.3 | 9.2 | 0.6 | 9.0 | 11.2 | 7.3 | 15.7 | 7.6 | 14.5 | |

20 | 9.0 | 13.2 | 1.0 | 16.5 | 25.4 | 12.9 | 44.1 | 13.4 | 39.2 | |

50 | 12.5 | 18.4 | 2.0 | 35.1 | 87.1 | 28.5 | 204.2 | 29.2 | 173.7 | |

100 | 15.2 | 22.3 | 3.1 | 62.3 | 253.9 | 53.7 | 719.7 | 54.6 | 598.3 | |

Subregion II | 5 | 3.2 | 4.9 | 0.2 | 4.7 | 5.4 | 4.2 | 6.1 | 4.3 | 5.9 |

10 | 5.7 | 8.7 | 0.5 | 8.8 | 11.5 | 7.5 | 14.9 | 7.8 | 13.9 | |

20 | 8.2 | 12.6 | 1.0 | 16.0 | 26.7 | 13.3 | 40.6 | 13.7 | 36.7 | |

50 | 11.5 | 17.6 | 1.9 | 33.8 | 96.1 | 29.0 | 182.6 | 29.7 | 158.2 | |

100 | 14.0 | 21.4 | 3.0 | 60.4 | 289.5 | 54.3 | 633.5 | 55.1 | 536.3 | |

Subregion III | 5 | 2.8 | 4.1 | 0.1 | 4.8 | 5.2 | 4.2 | 6.2 | 4.3 | 6.0 |

10 | 5.5 | 7.8 | 0.4 | 9.1 | 11.0 | 7.4 | 15.4 | 7.7 | 14.3 | |

20 | 8.1 | 11.6 | 0.7 | 16.9 | 24.5 | 13.1 | 42.6 | 13.6 | 38.1 | |

50 | 11.6 | 16.6 | 1.2 | 36.0 | 81.6 | 28.7 | 194.9 | 29.4 | 167.1 | |

100 | 14.2 | 20.4 | 1.7 | 63.7 | 232.0 | 54.0 | 682.7 | 54.8 | 571.7 | |

Subregion IV | 5 | 3.7 | 5.5 | 0.2 | 4.8 | 5.3 | 4.3 | 6.0 | 4.4 | 5.8 |

10 | 6.6 | 9.8 | 0.5 | 9.1 | 11.2 | 7.6 | 14.5 | 7.9 | 13.6 | |

20 | 9.5 | 14.0 | 0.8 | 16.6 | 25.1 | 13.4 | 39.2 | 13.9 | 35.7 | |

50 | 13.3 | 19.6 | 1.3 | 35.4 | 85.4 | 29.2 | 174.1 | 29.9 | 151.8 | |

100 | 16.1 | 23.9 | 1.8 | 62.7 | 246.9 | 54.5 | 599.6 | 55.4 | 511.0 | |

Subregion V | 5 | 3.1 | 4.8 | 0.1 | 4.8 | 5.3 | 4.3 | 5.9 | 4.5 | 5.6 |

10 | 6.2 | 9.5 | 0.6 | 9.1 | 11.1 | 7.7 | 14.2 | 8.2 | 12.8 | |

20 | 9.3 | 14.3 | 1.0 | 16.8 | 24.6 | 13.6 | 38.0 | 14.5 | 32.1 | |

50 | 13.3 | 20.5 | 1.7 | 35.9 | 82.2 | 29.4 | 166.1 | 31.0 | 129.6 | |

100 | 16.4 | 25.3 | 2.4 | 63.6 | 234.5 | 54.8 | 567.8 | 56.7 | 422.9 | |

YRSR | 5 | 3.1 | 4.8 | 0.2 | 4.7 | 5.3 | 4.2 | 6.1 | 4.3 | 5.9 |

10 | 5.5 | 8.5 | 0.6 | 8.9 | 11.4 | 7.5 | 15.0 | 7.8 | 14.0 | |

20 | 7.9 | 12.3 | 1.0 | 16.1 | 26.3 | 13.2 | 41.0 | 13.7 | 37.2 | |

50 | 11.1 | 17.2 | 1.8 | 34.1 | 93.3 | 28.9 | 185.0 | 29.6 | 161.5 | |

100 | 13.5 | 20.9 | 2.7 | 60.9 | 278.6 | 54.2 | 643.0 | 55.0 | 549.6 |

Note: The first column is the given return period of single variable DD, DS and DP. The second to the fourth columns are the value of DD, DS and DP with given return period, which can be obtained by inverse function of Equations (16)–(20). The sixth to eleventh column are the combined and co-occurrence return period, which can be obtained by Equations (30)–(31).

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. |

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).