Hydrological Drought Risk Assessment Using a Multidimensional Copula Function Approach in Arid Inland Basins, China

The aim of this research was to use the standardized runoff index (SRI) with a three-month timescale (SRI-3) to analyze hydrological drought risk in two arid river basins characterized by different runoff regimes, Northwest China. Based on SRI-3, hydrological drought levels for different events were defined through run theory. The hydrological drought risk in the two study basins was then comprehensively assessed using a multidimensional copula function that considered the multivariable joint probability of hydrological drought duration, severity, intensity and peak. Results indicate that: (1) the risk of hydrological drought in the two basins between 1961–2018 periodically changed. There was a slight increase in risk within the Yarkant River Basin, while there was a clear decrease in risk within the Kaidu River Basin. The magnitude of drought in the two basins was relatively low; both basins were dominated by mild to moderate hydrological droughts; (2) the drought probabilities of the Yarkant River Basin and Kaidu River Basin from 1961 to 2018 exhibited a falling-rising-falling pattern and a rising-falling trend through time, respectively. These trends were correlated with changes in precipitation and the area of glacial ice, which presumably influenced the amount and source of runoff in the two basins. Hydrological drought risk in the Yarkant River Basin was higher than in the Kaidu River Basin; and (3) the return period of mild, moderate, severe and extreme drought events was 2 yrs, 8 yrs, 20 yrs, and 60 yrs in the Yarkant River Basin, respectively, and 2 yrs, 8 yrs, 23 yrs and 74 yrs in the Kaidu River Basin, respectively.


Introduction
Drought is a natural phenomenon caused by below-normal precipitation over an extended period [1][2][3]. In response to global climate change, drought has become a serious natural hazard, affecting both the ecological environment and human life [4,5]. Droughts are complex events best characterized by a series of properties including their frequency, duration and intensity [6,7]. Droughts can also take on a variety of unique forms depending on which part of the hydrological cycle they impact most strongly. For example, they can be classified into four categories, including meteorological droughts, agricultural droughts, hydrological droughts, and socio-economic droughts [8]. Among these categories, a hydrological drought, defined as the deficit in surface or subsurface water, is the most important, considering the high dependence of social activities (e.g., irrigation, industry and urban living) on water resources from catchment impoundments especially in arid regions. risk was determined and analyzed. The results of this study provide valuable information for the administration of water in arid regions coping with climate change.

Study Area
In this paper, two important tributaries of the Tarim River, the Yarkant River and Kaidu River, were selected as typical arid region research basins. The Yarkant River Basin (76 • 40 -79 • 25 E to 37 • 20 -40 • 10 N) is located on the northern slope of the Karakoram Mountains ( Figure 1) and is the longest tributary of the Tarim River, which is the largest inland river in China [26]. Within the study area, the average annual precipitation measures about 53.14 mm, whereas the annual potential evaporation is much larger, measuring about 2196 mm. The annual average temperature is about 3.9 • C, which is typical of an inland arid climate. Modern glaciers and permanent snow cover occur in the mountainous areas of the basin. Glaciers in the basin cover 5853 km 2 , accounting for 11.95% of the basin area, whereas snow cover accounts for about 15% of the basin area. The basin is mainly recharged by ice and snow meltwater, where glacial and snow meltwater accounts for about 75% of the total runoff; interannual variations in runoff are significant [27,28]. The total annual runoff is about 6.875 billion m 3 [27]. The Yarkant river basin irrigation district is the fourth agricultural irrigation district in China [29], which is located in the arid region of central Asia. This region is characterized by a typical temperate continental arid climate, and is the most important region for the production of high-quality grain, cotton and fruits in Xinjiang in particular, and China in general. The Yarkant River Basin is a desert oasis that it is located between the Taklimakan Desert and the Bulgari and Togak Desert. Its terrestrial, riparian, and aquatic ecosystems are fragile, and may be impacted the region's low and irregular rainfall, high temperatures and evaporation, and notable periods of drought.
Water 2020, 12, x FOR PEER REVIEW 4 of 24 types of watersheds, which provide water to the downstream Tarim River, has very important practical implications to the sustainable use of water resources, the sustainable development of society and economy, and the restoration and reconstruction of the aquatic and riparian ecological environment in the watershed.

Data
In this study, the annual meteorological data (temperature and precipitation) from the stations at Tashkurgan and Bayanbulak in the upper reaches of the Yarkant River and the Kaidu River, respectively, from 1961 to 2016 were obtained from China Meteorological Administration (http://data.cma.cn/) to analyze recent, decadal-scale changes in climate within the two basins. Monthly runoff data from 1961-2018 from the Kaqun hydrographic station in Yarkant River Basin and the Dashankou hydrographic station in the Kaidu River Basin were obtained from Xinjiang Tarim River Basin Management Bureau to calculate the SRI, characterize hydrologic drought conditions and assess hydrologic drought risk.  (Figure 1). The average annual precipitation measures about 276 mm, the average annual temperature of the basin was only −4.16 • C, and the average maximum snow depth is 12 cm, which has obvious characteristics of alpine climate. River runoff is composed of a mixture of glacier meltwater, snow meltwater and rainwater, where precipitation accounts for about 45% of the runoff in the Kaidu River [30]. By contrast, glaciers Water 2020, 12, 1888 4 of 22 in the basin cover 333 km 2 , accounting for 1.75% of the basin area, glacial meltwater only accounts for only 14.1%. The total annual runoff is about 3.515 billion m 3 [31]. The Kaidu River is also one of the tributaries of the Tarim River (Figure 1), and it is the source of China's largest inland freshwater lake, Bosten Lake, with a total length of more than 500 km and a drainage area of 22,000 km 2 .

Methodology
The sources of runoff in both of these two river basins include glacial and snow meltwater and precipitation; the runoff varies greatly within and between years, as the arid areas are particularly sensitive to climate change. However, the contribution of water from each source differs between the two rivers so that the degree to which they are affected by climate change is different. At the same time, they are both inland rivers within an arid region, which are subjected to droughts. Therefore, the assessment of hydrological drought risk in these two types of watersheds, which provide water to the downstream Tarim River, has very important practical implications to the sustainable use of water resources, the sustainable development of society and economy, and the restoration and reconstruction of the aquatic and riparian ecological environment in the watershed.

Data
In this study, the annual meteorological data (temperature and precipitation) from the stations at Tashkurgan and Bayanbulak in the upper reaches of the Yarkant River and the Kaidu River, respectively, from 1961 to 2016 were obtained from China Meteorological Administration (http://data.cma.cn/) to analyze recent, decadal-scale changes in climate within the two basins. Monthly runoff data from 1961-2018 from the Kaqun hydrographic station in Yarkant River Basin and the Dashankou hydrographic station in the Kaidu River Basin were obtained from Xinjiang Tarim River Basin Management Bureau to calculate the SRI, characterize hydrologic drought conditions and assess hydrologic drought risk.

Calculation of the SRI and Hydrological Drought Classification
The standardized runoff index (SRI) is a widely used index to assess hydrologic drought conditions because of its simplicity and the availability of the required data [16,32]. In this study, the SRI was applied as a hydrological drought indicator to analyze hydrological droughts in the study basins. The SRI was calculated using the following approach.
If it is assumed that runoff x during a certain time period satisfies the probability density function f (x) for the Gamma distribution, then f (x) can be determined using the following equation: where γ and β are the shape and scale parameters, respectively, and x > 0, γ > 0, and β > 0. The parameters γ and β can be calculated using the maximum likelihood method. Γ(β) is the gamma function. The cumulative probability of runoff x for a given time scale is as follows: After normalizing the probability of the Γ distribution, SRI can be calculated: Water 2020, 12, 1888

of 22
where Γ is the runoff distribution probability that is associated with the Γ function, and S is the probability density coefficient. When F > 0.5, S = 1; when F ≤ 0.5, S = −1; c 0 = 2.515517, c 1 = 0.802853, c 2 = 0.010328, d 1 = 1.432788, d 2 = 0.189269, d 3 = 0.001308; the values are empirical values. SRI has a characteristic of multiple time scales, which reflect hydrological drought phenomena in different time scales and can be represented by SRI-m, where m is month, for example, SRI-1, SRI-3, SRI-6, SRI-12, respectively.

SRI Value Drought Intensity
Severe drought ≤−2 Extreme drought

The Identification of Hydrological Drought Variables
Run Theory has been extensively employed to identify drought events. In this study, hydrological drought duration (Dd), hydrological drought severity (Ds), hydrological drought intensity (Di) and hydrological drought peak (Dp) were selected as the variables to characterize hydrological drought conditions needed for hydrological drought risk analysis. In accordance with run theory (Figure 2), each variable was identified as follows: Dd represents the duration of a hydrological drought event: the count of continuous months at which the value of the SRI is below the threshold X 0 .
Ds represents the severity of a hydrological drought event: the absolute sum of all SRIs during the drought period.
Di represents the intensity of a hydrological drought event, which is the ratio of drought severity to drought duration. Di represents the average drought water shortage during the drought event.
Dp represents the peak value of a hydrological drought event: the minimum value of the SRI during the drought period.
Using the above definitions, hydrological drought identification was carried out by calculating the SRI over a 3-month scale. The threshold was set to be −0.5.
Water 2020, 12, x FOR PEER REVIEW 6 of 24 Figure 2. Run theory analysis was used to determine the hydrological drought variables. Numbers along the x-axis represent identified drought events. Dp, drought peak; Ds, drought severity; Dd, drought duration.

Hydrological Drought Risk Probability Model
The most important step before the hydrological drought risk assessment is to carry out hydrological drought identification. In this paper, the characteristics of a hydrological drought event, including Dd, Ds, Di and Dp, were calculated on the basis of run theory, as shown in Figure 2. Run theory analysis was used to determine the hydrological drought variables. Numbers along the x-axis represent identified drought events. Dp, drought peak; Ds, drought severity; Dd, drought duration.

Hydrological Drought Risk Probability Model
The most important step before the hydrological drought risk assessment is to carry out hydrological drought identification. In this paper, the characteristics of a hydrological drought event, including Dd, Ds, Di and Dp, were calculated on the basis of run theory, as shown in Figure 2. In addition, the marginal distributions of these variables were fitted. Once the above had been accomplished, a copula function was adopted to build the hydrological drought risk probability model.

The Theory of Copula
The copula function was initially used by Sklar [34]. Copula function models can be used to construct a multidimensional joint distribution by using marginal distributions and a correlation framework [35]. A variety of copula functions can be used to set up a multidimensional joint distribution of drought variables. The Archimedean copula function is one of the most commonly used copula functions in hydrology [36,37]; this type of copula includes the Gumbel-Hougaard (GH), Clayton, and Frank functions. These three forms of copula have been commonly selected because their associated correlation models perform well. The correlation coefficient of the Gumbel copula is 2 − 2 1/θ , whereas the correlation coefficient of the lower tail is zero, which enhances its ability to describe the variation of two random variables with upper tail correlation. The Frank copula is a symmetry correlation function in the Archimedean copula family. The correlation coefficients of the upper and lower tails are equal and all zero. The lower tail coefficient of the Clayton copula is 2 1/θ , and the upper tail correlation coefficient is zero. Therefore, the characteristics of these three copula functions are representative and can explain the problem from different aspects.
(1) Parameter Estimation A nonparametric estimation method was adopted to compute the parameters of the copula in this study. This technique is primarily related to the copula parameter (θ) ( Table 2). The relationship between θ and τ (Kendall correlation coefficient) is represented by the following equation: After calculating τ from the measured data, the parameters of the joint distribution can be obtained.
(2) Verification and Evaluation In order to quantitatively evaluate the fitting error and select the appropriate copula function, the root mean square error (RMSE) and the Akaike information criterion (AIC) were used [38]: where m is the number of model parameters, n is the number of samples, P i represents the copula value of consecutive sample observations, and Pe i represents the corresponding multivariate empirical probability. AIC is a measure of the quality of the statistical model fit to the data. For a particular copula function, the smaller the AIC value of the objective function value, the better the copula function simulation. We defined RMSE as: where y i (θ) is the simulated two-variable joint probability value, y i represents the empirical observation, i is the serial number of the variable, and n is the total number of variables. The range of RMSE is [0, ∞]; RMSE is equal to 0 for a perfect model.

Correlation Analysis and Establishment of Marginal Distribution Function
In order to determine whether the drought variables were correlated and well-suited for establishing the joint distribution function, the correlations between the drought variables were analyzed using the Kendall rank, Spearman's rank-order, and Pearson product-moment correlation coefficients. The correlations were calculated using the following equations: where X and Y are random variables, the numerator represents the sample covariance while the denominator is the product of the sample standard deviations.
Spearman correlation coefficients : where di is the difference between the ranks of corresponding values X and Y, and n is the number of values in each data set (same for both sets).
Kendall correlation coefficients : where C is the number of the pairs that have consistency in X, Y, D is the number of the inconsistent pairs in X, Y. Twenty different marginal distributions, including the exponential, generalized Pareto, generalized extreme value, gamma, Weibull, and Gaussian distributions, among others, were fit to the parameters Dd, Ds, Di and Dp. The parameters were estimated using a maximum likelihood approach and the Chi-square test for the goodness-of-fit. The marginal distribution of each variable was determined by assuming that the threshold value should be as small as possible to preserve the largest sample.

Joint Probability Distribution and Hydrological Drought Risk Assessment
In combination with the definition used herein for a hydrological drought, the drought risk probability was defined as the joint probability of Dd, Ds, Di and Dp. The analysis calculates the marginal distribution for each drought variable and obtains the specific parameters of the function. The fitting of the function to the actual data was compared using a quantile-quantile (Q-Q) plot. A two-dimensional symmetrical Archimedean copula function was used to construct a 2D copula which was fit to the joint distribution functions of the four indicators (Dd, Ds, Di and Dp) by using an inversion of the Kendall's tau method in the two watersheds. Then, by comparing the RMSE and AIC values of the different copula functions, the most suitable copula function for each set of variables was selected. Construction of the 3D copula function was also based on the previously obtained univariate marginal distribution function and by applying the same calculation as used for the 2D models. The most suitable functions were then selected on the principle of minimizing the RMSE and AIC. The specific formulas are in Table 2. Table 2. Symmetric Archimedean copula functions.

The Return Period of a Drought Event
The return period refers to the length of time between the reoccurrence of a random variable over a specified period. The calculation of the hydrological drought return period can provide valuable information for the rational utilization of water [39]. Use of a single-variable and univariate return period often leads to an overestimation or underestimation of the risk rate for a given event. Thus, this study calculated the two-and three-variable return periods for hydrological drought risk analysis.
The approach assumes that X 1 , X 2 and X 3 represent the sequence of characteristic variables with dependent relationships to drought events. The marginal distribution functions were u 1 = F(x 1 ), u 2 = F(x 2 ), and u 3 = F(x 3 ), while the univariate return period of drought eigenvalues was expressed as: where L represents the time interval over which a drought occurs, and E(L) represents the average time interval between drought events. E(L) = N/n, where N is the length of the time series of drought analysis, and n is the number of times of drought occurrence. The two-dimensional joint return period can be expressed as: The two-dimensional co-occurrence return period can be expressed as: The three-dimensional joint return period can be expressed as: The three-dimensional co-occurrence return period can be expressed as: where T or represents the joint return period of the drought elements, T and represents the co-occurrence return period of drought variables, and ) represents the joint probability of a three-dimensional variable group.

Changes in Climate within the Study Area
Based on the observation data from the meteorological stations in the two river basins, the annual precipitation and temperature between 1961 and 2016 were analyzed using the M-K test method (Figures 3 and 4). The results show that from 1961 to 2016, the temperature in the Yarkant and Kaidu river basins increased at rates of 0.027 • C/yr and 0.021 • C/yr, respectively. Precipitation increased at a rate of 0.51 mm/yr and 1.19 mm/yr, respectively. In addition, to the M-K test results of precipitation and air temperature and mutation curves of temperature and precipitation showed an overall upward trend since 1961. However, the trend was not significant at the 0.05 level, indicating that the trend in climate change within the research basins was not significant (Figure 4).
Water 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/water Yarkant and Kaidu river basins increased at rates of 0.027 °C/yr and 0.021 °C/yr, respectively. Precipitation increased at a rate of 0.51 mm/yr and 1.19 mm/yr, respectively. In addition, to the M-K test results of precipitation and air temperature and mutation curves of temperature and precipitation showed an overall upward trend since 1961. However, the trend was not significant at the 0.05 level, indicating that the trend in climate change within the research basins was not significant ( Figure 4).

Changing Temporal Trend in the SRI
Based on the 58-year runoff record (1961-2018) from the Yarkant and Kaidu river basins, the average standardized runoff index of different time scales (SRI1-SRI12) for the two river basins were calculated ( Figures 5 and 6). The figures show that the SRI of the two river basins exhibits a tendency to rise. However, the fluctuations in the SRI of the two river basins are different. These differences reflect differences in hydrological drought conditions between the two adjacent rivers. The fluctuations of different time scales SRI are also different, on the time scale of 1-6 months (SRI-1 to SRI-6), although there were obvious positive and negative alternations in the SRI value in the basin, the SRI value was mainly negative before the year of abrupt change, and gradually increased after that. On the time scale of 6-12 months (SRI6-SRI12), the positive and negative differentiation of SRI values became more and more obvious after the mutation year in two basins.

Changing Temporal Trend in the SRI
Based on the 58-year runoff record (1961-2018) from the Yarkant and Kaidu river basins, the average standardized runoff index of different time scales (SRI1-SRI12) for the two river basins were calculated ( Figures 5 and 6). The figures show that the SRI of the two river basins exhibits a tendency to rise. However, the fluctuations in the SRI of the two river basins are different. These differences reflect differences in hydrological drought conditions between the two adjacent rivers. The fluctuations of different time scales SRI are also different, on the time scale of 1-6 months (SRI-1 to SRI-6), although there were obvious positive and negative alternations in the SRI value in the basin, the SRI value was mainly negative before the year of abrupt change, and gradually increased after that. On the time scale of 6-12 months (SRI6-SRI12), the positive and negative differentiation of SRI values became more and more obvious after the mutation year in two basins.
After performing the M-K tests on the multi-year, average SRI values of the rivers (Figure 4), it is apparent that the SRI of the Yarkant River Basin mutated in 1991 and 1993; the SRI in the Kaidu River Basin mutated in 1994. The SRI of the two basins before the mutation year were generally negative, while more positive values occurred after the mutation year. Changes in the trend for the Yarkant River Basin were more significant. Furtheremore, a comparison of the changing trends in the SRI and the climatic factors within the two river basins shows that the trends in SRI and precipitation within the Kaidu River Basin were highly correlated (r = 0.69), while there was no clear correlation between the trends in SRI and precipitation within the Yarkant River Basin.
Water 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/water basins were calculated ( Figures 5 and 6). The figures show that the SRI of the two river basins exhibits a tendency to rise. However, the fluctuations in the SRI of the two river basins are different. These differences reflect differences in hydrological drought conditions between the two adjacent rivers. The fluctuations of different time scales SRI are also different, on the time scale of 1-6 months (SRI-1 to SRI-6), although there were obvious positive and negative alternations in the SRI value in the basin, the SRI value was mainly negative before the year of abrupt change, and gradually increased after that. On the time scale of 6-12 months (SRI6-SRI12), the positive and negative differentiation of SRI values became more and more obvious after the mutation year in two basins.  changing trends in the SRI and the climatic factors within the two river basins shows that the trends in SRI and precipitation within the Kaidu River Basin were highly correlated (r = 0.69), while there was no clear correlation between the trends in SRI and precipitation within the Yarkant River Basin. The SRI values exhibit different time scales (Figure 7a,b). Changes observed for SRI-1 are characterized by short intervals and frequent alternations between positive and negative values, which may lead to a misestimation of the occurrence of drought events. In contrast, SRI-3 possessed three-month intervals, which more closely reflected seasonal changes in the river basin. SRI-12 reflects the regularity of drought events throughout the year. Figure 7a shows that the general trends in SRI-3 and SRI-12 within the Yarkant River Basin from 1961 to 2018 were similar. Both of the two basins exhibited an increase in the SRI around 1993; values after 1993 were more positive. However, SRI-3 defined more droughts, separated by wet intervals, than SRI-12, which was characterized by more frequent changes from positive to negative SRI values. Figure 7b shows that the trends illustrated by SRI-3 and SRI-12 in the Kaidu River Basin from 1961 to 2018 were similar. SRI values of the basin were more positive around 1995, but the duration of droughts shown by SRI-3 was shorter and there were more frequent alterations between dry and wet periods than shown by SRI-12. From the comparison of the two basins, the general changes in the trend in SRI values between the Yarkant and Kaidu rivers were similar. Both basins have become wetter since the late 1990s, but the exact timing of the droughts and wet periods differed. In light of the above results, the seasonal hydrological characteristics of the watersheds, as well as drought risk, were subsequently analyzed using The SRI values exhibit different time scales (Figure 7a,b). Changes observed for SRI-1 are characterized by short intervals and frequent alternations between positive and negative values, which may lead to a misestimation of the occurrence of drought events. In contrast, SRI-3 possessed three-month intervals, which more closely reflected seasonal changes in the river basin. SRI-12 reflects the regularity of drought events throughout the year. Figure 7a shows that the general trends in SRI-3 and SRI-12 within the Yarkant River Basin from 1961 to 2018 were similar. Both of the two basins exhibited an increase in the SRI around 1993; values after 1993 were more positive. However, SRI-3 defined more droughts, separated by wet intervals, than SRI-12, which was characterized by more frequent changes from positive to negative SRI values. Figure 7b shows that the trends illustrated by SRI-3 and SRI-12 in the Kaidu River Basin from 1961 to 2018 were similar. SRI values of the basin were more positive around 1995, but the duration of droughts shown by SRI-3 was shorter and there were more frequent alterations between dry and wet periods than shown by SRI-12. From the comparison of the two basins, the general changes in the trend in SRI values between the Yarkant and Kaidu rivers were similar. Both basins have become wetter since the late 1990s, but the exact timing of the droughts and wet periods differed. In light of the above results, the seasonal hydrological characteristics of the watersheds, as well as drought risk, were subsequently analyzed using the SRI-3 data. The run theory applied in this study divided the frequency of drought events in the two study basins into four categories of differing magnitude (Figure 8). Between 1961 to 2018, both the Yarkant and Kaidu river basins were dominated by mild and moderate hydrological droughts. Extreme droughts were rare. In addition, the frequency of extreme drought in the Kaidu River Basin was significantly lower than in the Yarkant River Basin. As the SRI time scale increased, the frequency of mild and moderate droughts decreased, while the frequency of severe and extreme droughts increased. These trends indicate that both severe and extreme droughts were associated with prolonged droughts. The above research results show that the drought characteristics in the two watersheds are different. The run theory applied in this study divided the frequency of drought events in the two study basins into four categories of differing magnitude (Figure 8). Between 1961 to 2018, both the Yarkant and Kaidu river basins were dominated by mild and moderate hydrological droughts. Extreme droughts were rare. In addition, the frequency of extreme drought in the Kaidu River Basin was significantly lower than in the Yarkant River Basin. As the SRI time scale increased, the frequency of mild and moderate droughts decreased, while the frequency of severe and extreme droughts increased. These trends indicate that both severe and extreme droughts were associated with prolonged droughts. The above research results show that the drought characteristics in the two watersheds are different.

Correlation Analysis
After the four hydrological drought variables were extracted for each river, the Kendall rank, Spearman's rank-order, and Pearson product-moment correlation coefficients were used to measure the correlation between hydrological drought variables in each river basin (Table  3). Coefficients closer to 1 indicate stronger correlation between the two variables. The results show that the correlation coefficients among the groups of variables in both basins were close to 1. Thus, these four drought variables were capable of establishing the joint distribution function using the copula function to analyze the hydrological drought characteristics in two study basins.

Selection of a Suitable Marginal Distribution
Selection of a suitable marginal distribution function for each variable is required to effectively construct the copula joint distribution function. In this study, the most suitable marginal distribution of the four variables based on a BIC goodness-of-fit analysis (Bayesian Information Criterion) was initially selected. Then 20 functions (e.g., extreme value, generalized Pareto, exponential, and normal) were used to fit the marginal probability

Correlation Analysis
After the four hydrological drought variables were extracted for each river, the Kendall rank, Spearman's rank-order, and Pearson product-moment correlation coefficients were used to measure the correlation between hydrological drought variables in each river basin (Table 3). Coefficients closer to 1 indicate stronger correlation between the two variables. The results show that the correlation coefficients among the groups of variables in both basins were close to 1. Thus, these four drought variables were capable of establishing the joint distribution function using the copula function to analyze the hydrological drought characteristics in two study basins.

Selection of a Suitable Marginal Distribution
Selection of a suitable marginal distribution function for each variable is required to effectively construct the copula joint distribution function. In this study, the most suitable marginal distribution of the four variables based on a BIC goodness-of-fit analysis (Bayesian Information Criterion) was initially selected. Then 20 functions (e.g., extreme value, generalized Pareto, exponential, and normal) were used to fit the marginal probability distribution of each variable. Finally, the parameters were adjusted by using the maximum likelihood method to obtain the best-fit function of each variable, as shown in Table 4. The Chi-square test was also applied for four variables at a significance level of 5%. The analysis showed that the marginal distribution function of Dd in the Yarkant River Basin was a T Location Scale, the marginal distribution function of Ds was a logistic function, and the marginal distribution functions of Di and Dp were a generalized extreme value function. The marginal distribution function of Dd in the Kaidu River Basin was exponential, whereas the marginal distribution function of Di was a T location scale function. The marginal distribution functions of Ds and Dp were generalized extreme value functions. Figure 9 shows the fitting map and quantile-quantile (Q-Q) plots of the drought variables' marginal distributions for both rivers. These plots demonstrate that the fitting results were acceptable.

Selection of the Suitable Copula
The Gumbel, Frank and Clayton copula functions were chosen to make up the pairwise 2D joint distribution of Dd, Ds, Di and Dp. The parameters were calculated using an inversion of the Kendall's tau method. The RMSE and AIC values for each group of joint distributions were also calculated. The RMSE and AIC values of the copulas are provided in Table 5. For a particular copula function, the smaller the AIC value of the objective function value, the closer the RMSE value is to 0, and the better the copula function simulates the data. According to this criterion, a suitable copula function for each group of variables for the two basins were selected and are highlighted in bold (Table 5). Table 5 shows that the Frank copula provides a good fit for Di_Dp, Dd_Di and Dd_Dp for the Yarkant River, and Dd_Di for the Kaidu River. The Clayton copula fit the data well for Di_Dp and Dd_Dp of the Kaidu River. The comparison between the empirical frequency and the theoretical frequency of the copula functions is shown   The curves fitted to the drought variables demonstrate that the cumulative probability of drought duration (Dd), drought severity (Ds), drought intensity (Di), and drought peak (Dp) changed systematically at the top, implying long duration drought events (10-20 M); average drought severity values of 15-25 and drought peak values between 2-2.5 were rare. The cumulative probability curve for drought severity (Ds) was the steepest, mainly concentrated between 0 and 10. The probability distribution of drought severity was relatively uniform; there was no obvious area of concentration.

Selection of the Suitable Copula
The Gumbel, Frank and Clayton copula functions were chosen to make up the pairwise 2D joint distribution of Dd, Ds, Di and Dp. The parameters were calculated using an inversion of the Kendall's tau method. The RMSE and AIC values for each group of joint distributions were also calculated. The RMSE and AIC values of the copulas are provided in Table 5. For a particular copula function, the smaller the AIC value of the objective function value, the closer the RMSE value is to 0, and the better the copula function simulates the data. According to this criterion, a suitable copula function for each group of variables for the two basins were selected and are highlighted in bold (Table 5). Table 5 shows that the Frank copula provides a good fit for Di_Dp, Dd_Di and Dd_Dp for the Yarkant River, and Dd_Di for the Kaidu River. The Clayton copula fit the data well for Di_Dp and Dd_Dp of the Kaidu River. The comparison between the empirical frequency and the theoretical frequency of the copula functions is shown in Figure 10. The Frank copula and Clayton copula were the most suitable copula functions for the 2D joint distribution in the study basins.
The 2D joint distribution probability can only represent the relationship and probability of two hydrological drought variables. However, in order to more comprehensively describe and analyze hydrological drought risk, it is essential to set up the multi-dimensional joint distribution of the drought variables. Therefore, a three-dimensional (3D) symmetric copula was selected to form the correlation structure of the 3D hydrological drought variables. Three of the four drought variables (Dd, Ds, Di and Dp) were chosen to form a group. Specifally, three drought variables were chosen from four drought variables to form four pairs of copula three-dimensional joint distributions for the two watersheds. The copula parameters were calculated using a inversion of Kendall's tau method (Table 6). Since the smaller the RMSE value, the better the fit of the model, the three-dimensional joint copula distribution established by Frank copula provided a better result for Dd_Di_Dp. Thus, the probalibity of hydrological drought risk was analyzed using the Dd_Di_Dp 3D joint distribution construct using the Frank copula function for the two study basins.

Drought Risk Probability Assessment
The joint probability of the copula function was used as the evaluation index of hydrological drought risk and was analyzed using the 2D and 3D joint distribution probabilities of the hydrological drought variables, as described below.
(1) Two-Dimensional Joint Distribution Functions Figure 11 shows the joint probability distribution of each two variable group used to describe droughts. The two-dimensional joint probability (drought risk) differs from the three variable groups. Specifically, the joint probability of Di_Dp for the Yarkant River Basin was evenly distributed between 0.01-0.9, whereas the joint probability distribution of Dd_Di and Dd_Dp was between 0.01-0.2 and 0.7-0.9. The joint probability of Dd_Di and Dd_Dp in the Kaidu River Basin was evenly distributed between 0.01 and 0.9, and the joint probability of Di_Dp was mainly distributed between 0.01-0.3 and 0.7-0.9. In addition, for the 2D joint distribution of hydrological drought variables, when the value of one hydrological drought variable was fixed, the value of another variable increased. In addition, the joint probability of both increased toward a specific value. For example, when Dd was fixed, the greater the value of Di or Dp, then the greater the drought joint probability of Dd_Di or Dd_Dp.
Water 2020, 12, x FOR PEER REVIEW 18 of 24

Drought Risk Probability Assessment
The joint probability of the copula function was used as the evaluation index of hydrological drought risk and was analyzed using the 2D and 3D joint distribution probabilities of the hydrological drought variables, as described below.
(1) Two-Dimensional Joint Distribution Functions Figure 11 shows the joint probability distribution of each two variable group used to describe droughts. The two-dimensional joint probability (drought risk) differs from the three variable groups. Specifically, the joint probability of Di_Dp for the Yarkant River Basin was evenly distributed between 0.01-0.9, whereas the joint probability distribution of Dd_Di and Dd_Dp was between 0.01-0.2 and 0.7-0.9. The joint probability of Dd_Di and Dd_Dp in the Kaidu River Basin was evenly distributed between 0.01 and 0.9, and the joint probability of Di_Dp was mainly distributed between 0.01-0.3 and 0.7-0.9. In addition, for the 2D joint distribution of hydrological drought variables, when the value of one hydrological drought variable was fixed, the value of another variable increased. In addition, the joint probability of both increased toward a specific value. For example, when Dd was fixed, the greater the value of Di or Dp, then the greater the drought joint probability of Dd_Di or Dd_Dp.  show that there are obvious differences in the two-dimensional joint probability between two basins. In terms of multi-year average joint probability, the 2D joint probability of hydrological drought variables in the Yarkant River Basin was larger than that in the Kaidu River Basin. The joint probability of Dd_Dp in the Yarkant River Basin was the largest (0.53), followed by Di_Dp  obvious differences in the two-dimensional joint probability between two basins. In terms of multi-year average joint probability, the 2D joint probability of hydrological drought variables in the Yarkant River Basin was larger than that in the Kaidu River Basin. The joint probability of Dd_Dp in the Yarkant River Basin was the largest (0.53), followed by Di_Dp and Dd_Di (0.50). The joint probability of Di_Dp in Kaidu River Basin exhibited the largest joint probability (0.45), followed by the Dd_Dp and Dd_Di (0.37). In terms of the trend in joint probability for multiple years, the joint probability of the Yarkant River Basin showed a falling-rising-falling pattern from 1961 to 2018; that is, a downward trend from 1961-1969 and from 1990-1999, and an upward trend from 2000-2009, after which it declined again. Data from the Kaidu River Basin increased from 1961-1969 and from 1970-1979, and generally decreased after 1970-1979. The 2D joint probability of hydrological drought variables in the Yarkant River Basin was larger than that in the Kaidu River Basin. Further, the fluctuations in hydrological drought and joint distribution probability values in the Yarkant River Basin was more obvious, indicating that the overall drought risk in the Yarkant River Basin was larger than that in the Kaidu River Basin.
Water 2020, 12, x FOR PEER REVIEW 19 of 24 1970REVIEW 19 of 24 -1979. The 2D joint probability of hydrological drought variables in the Yarkant River Basin was larger than that in the Kaidu River Basin. Further, the fluctuations in hydrological drought and joint distribution probability values in the Yarkant River Basin was more obvious, indicating that the overall drought risk in the Yarkant River Basin was larger than that in the Kaidu River Basin. (2) Three-Dimensional Joint Distribution Functions In order to obtain comprehensive hydrolgocial drought risk information, the 3D joint distribution probability was developed using three drought variables. Figure 12 displays the temporal variations of the joint distribution probability of Dd_Di_Dp in the two watersheds from 1961 to 2018. The hydrological drought risk in the 3D joint distribution is similar to the 2D analysis. The 3D joint distribution probability in the Yarkant River Basin was greater than that in the Kaidu River Basin; the value in the Yarkant River Basin is 0.52 as compared to 0.36 for the Kaidu River Basin. In terms of the changes through time, the 3D joint distribution probability, the hydrological drought risk in the Yarkant River Basin generally increased (0.007/10 yrs); however, it decreased in the Kaidu River Basin (0.035/10 yrs). The 3D joint distribution probability of two river basins also possessed a falling-rising-falling and risingfalling trend, as did the 2D joint probability (Figure 12), that is, a downward trend from 1961-1969 to 1990-1999, and an upward trend from 2000-2009, after which it declined again. Data from the Kaidu River Basin increased from 1961-1969 to 1970-1979, and generally decreased after 1970-1979, then it keep stable from 1990-1999.

Drought Event Return Period
The drought event return period is a significant component of drought risk analysis. Therefore, this paper further analyzed the return period of hydrological drought events identified by run theory in the two river basins. Specifically, the return period (T) of the single variable was defined as 2 yrs, 5 yrs, 10 yrs, 20 yrs, 50 yrs and 100 yrs, and the co-occurrence return period and joint return period of the two-variable and three-variable joint distributions were obtained by calculating the univariate drought eigenvalues.
The results showed that the univariate return period was between the joint return period and the co-occurrence return period, indicating that the joint return period and the cooccurrence return period can be used as the upper and lower thresholds of the univariate return (2) Three-Dimensional Joint Distribution Functions In order to obtain comprehensive hydrolgocial drought risk information, the 3D joint distribution probability was developed using three drought variables. Figure 12 displays the temporal variations of the joint distribution probability of Dd_Di_Dp in the two watersheds from 1961 to 2018. The hydrological drought risk in the 3D joint distribution is similar to the 2D analysis. The 3D joint distribution probability in the Yarkant River Basin was greater than that in the Kaidu River Basin; the value in the Yarkant River Basin is 0.52 as compared to 0.36 for the Kaidu River Basin. In terms of the changes through time, the 3D joint distribution probability, the hydrological drought risk in the Yarkant River Basin generally increased (0.007/10 yrs); however, it decreased in the Kaidu River Basin (0.035/10 yrs). The 3D joint distribution probability of two river basins also possessed a falling-rising-falling and rising-falling trend, as did the 2D joint probability (Figure 12), that is, a downward trend from 1961-1969 to 1990-1999, and an upward trend from 2000-2009, after which it declined again. Data from the Kaidu River Basin increased from 1961-1969 to 1970-1979, and generally decreased after 1970-1979, then it keep stable from 1990-1999.

Drought Event Return Period
The drought event return period is a significant component of drought risk analysis. Therefore, this paper further analyzed the return period of hydrological drought events identified by run theory in the two river basins. Specifically, the return period (T) of the single variable was defined as 2 yrs, 5 yrs, 10 yrs, 20 yrs, 50 yrs and 100 yrs, and the co-occurrence return period and joint return period of the two-variable and three-variable joint distributions were obtained by calculating the univariate drought eigenvalues.
The results showed that the univariate return period was between the joint return period and the co-occurrence return period, indicating that the joint return period and the co-occurrence return period can be used as the upper and lower thresholds of the univariate return period. Taking the univariate 10-yr return period as an example, the three-variables joint return period of the Yarkant River Basin was approximately 5.92 yrs, while the co-occurrence return period was at least 30.61 yrs.
The calculated return periods (Table 7), combined with the classification criteria of hydrological drought, show that the Dp of the Yarkant River was 0.91 and the Dp of Kaidu River was 0.98, which are in the range of mild drought. The magnitude of hydrologic drought in the two basins was relatively low, and the degree of drought in the Kaidu River Basin was higher than that in the Yarkant River Basin. When T = 2 yrs, the co-occurrence return periods of Di_Dp, Dd_Di, Dd_Dp, and Dd_Di_Dp in the Yarkant River Basin were calculated to be 2.09 yrs, 2.41 yrs, 2.27 yrs, and 2.51 yrs, respectively. The return period of a mild drought event lasting for 2.18 months was about 2 yrs, whereas the co-occurrence return periods of Di_Dp, Dd_Di, Dd_Dp, and Dd_Di_Dp in the Kaidu River Basin were estimated to be 2.12 yrs, 2.6 yrs, 2.39 yrs, and 2.5 yrs, respectively. When T = 5 yrs, the drought duration in the Yarkant River was about 3.85 yrs, and the Dp was about 1.43 yrs. At this time, the Dd_Dp co-occurrence return period was 7.32 yrs, indicating that the return period of a moderate drought event lasting for 3.85 months in the Yarkant River Basin was about 8 yrs. Similarly, the return period of a moderate drought event lasting 8.61 months in Kaidu River Basin was about 8 yrs. When T = 10 yrs, the return period of a severe drought event lasting 5.36 months in the Yarkant River Basin was about 20 yrs, and the return period of a severe drought event lasting 12.3 months in the Kaidu River Basin was about 23 yrs. When T = 20 yrs, the return period of an extreme drought event lasting 7.59 months in the Yarkant River Basin was about 60 yrs, and the return period of an extreme drought event lasting 16.02 months in the Kaidu River Basin was about 74 yrs. Note: Dd represents hydrological drought duration, Di represents hydrological drought intensity, Dp represents hydrological drought peak. (I_P)T and represents joint return period of Di and Dp, (I_P)T or represents co-occurrence return period of Di and Dp.

Application of the SRI
This paper calculated the standardized runoff index (SRI) based on monthly runoff data. The results show that the minimum SRI values in the Yarkant River Basin appeared in 1965, 1972, 1980, 1989, and 1993, indicating that drought events occurred during these years. According to the drought levels defined in Table 1, extreme drought events in the basin occurred in 1965 and 1993. These dates are consistent with the droughts recorded in the Yarkant River Basin Conservancy Journal, indicating that the results of the SRI analysis in the basin are reasonable, and SRI is an appropriate hydrological drought index to be used in these study basins.

The Application of Copula Functions
The hydrological drought risk determined using the 3D joint distribution is similar to the results produced using the 2D distributions. These similarities verify the rationality of selecting the 3D copula functions in this study. The joint probability of multidimensional variables can describe the correlation structure between variables, and can be well-applied to the study of hydrological drought in arid areas with complex hydrological processes.

Differences in Hydrological Drought Risk in Two Study Basins
The cause(s) of the differences in the hydrological drought risk (probability) are complex. Not only are they caused by reduced precipitation, but also because other factors affect the water balance in these watersheds [40]. In arid and semi-arid regions, the response of river runoff to climate and cryosphere changes are complicated. By comparing the hydrological drought risk of the two river basins, fluctuations in hydrological drought risk can be identified, which may reveal the impact of climate change in each basin. According to the calculated SRI values and the climatic factors analyzed above, the difference in hydrological drought risk between the basins was mainly caused by the different recharge sources of runoff in the two basins. The runoff in the Yarkant River Basin is derived predominantly from glaciers and snow meltwater [26], whereas runoff in the the Kaidu River Basin is derived from a mixture of ice and snow meltwater and precipitation. Therefore, temporal trends in the SRI and precipitation in the Kaidu River Basin are similar, while the correlation between the SRI and the climatic variables (temperature and precipitation) is low in the Yarkant River Basin. The changes in runoff in the Yarkant River Basin was mainly caused by glacier melt runoff at middle and high altitudes. Based on the First and Second Chinese Glacier Inventory (CGI, http://westdc.westgis.ac.cn), the Randolph Glacier Inventory (RGI 6.0, http://www.glims.org/RGI/), Landsat 8, and Landsat TM data, the fourth-stage data of glacier area in the Yarkant River Basin were interpreted and revised, and it was found that the area of glacial coverage in the Yarkant River Basin intially decreased and then increased, which is consistent with the trends in hydrological drought in the basin. Later, during a decrease in the area of glacial ice, a reduction in meltwater resulted in a reduction in runoff, which, in turn, increased the probability of hydrological drought.

Drought Risk and Severity in the Study Basins
Drought is a significant threat to oasis irrigated agriculture [41]. The lower region of the Yarkant River Basin is the largest oasis in Xinjiang, with an irrigated area of 4500 km 2 and a population of about 1.7 million people. The Yarkant River provides ecological water for downstream areas, and also provides a large amount of water for agricultural production and domestic use within the downstream oasis. The Kaidu River irrigation area is only 850 km 2 . The Kaidu River primarily provides ecological water to the Kongque River and Bosten Lake. This analysis demonstrated that the Yarkant River has a greater risk of hydrological drought. Although the degree of drought is mainly mild, increases in downstream agriculture and domestic water consumption has exacerbated the shortage of water resources. It is thus necessary for watershed managers to put forward coping strategies according to the hydrological drought characteristics and drought return periods within the analyzed watersheds.

Conclusions
In this paper, a hydrological drought index based on the SRI was determined for two river basins in an arid region using runoff data collected at hydrological stations. Four drought variables (Ds, Dd, Di and Dp) were defined based on SRI-3 and run theory. Their joint probability was constructed using a copula function, which was used to develop an index of hydrological drought risk in the study basins. The hydrological drought risk probablity and hydrological drought event return period in the basins during the past half century were analyzed. The analysis resulted in the following conclusions: (1) Temporal changes in climate were analyzed from 1961 to 2016 in the study basins. Termporal trends were not statistically significant, although a slight increase in temperature was apparent in data from the Yarkant River Basin (0.027 • C/yr), whereas precipitation increased in Kaidu River Basin (1.19 mm/yr). The SRI of the two basins before the mutation year were generally negative, while more positive values occurred after the mutation year. The general temperature and precipitation patterns through time within the two basins was similar.
(2) The drought risk determined using the three-dimensional joint distribution was similar to the results obtained using the two-dimensional joint distribution; the probability of the two-variable joint distribution was greater than the probability of the three-variable joint distribution. Specifically, joint probability of hydrological drought risk within the Yarkant and Kaidu river basins showed a falling-rising-falling and rising-falling trend from 1961 to 2018, respectively. The joint probability of hydrological drought risk in the Yarkant River Basin was greater than in the Kaidu River Basin, and the degree of variation in joint probability in the Yarkant River Basin was slightly greater than that in the Kaidu River Basin.
(3) The return period of the mild, moderate, severe and extreme hydrological drought events in the two study basins were calculated. The co-occurrence return period of various levels of drought in the Yarkant River Basin was shorter than that of the Kaidu River Basin, and the joint return period was larger than that of the Kaidu River Basin. The single-variable return period of droughts in the two river basins was between the joint return period and the co-occurrence return period. The co-occurrence return period and the joint return period can be used as upper and lower thresholds of the single variable return period under the same drought-level conditions. The results of this study provide a scientific basis for the analysis and prediction of drought risk in the Yarkant and Kaidu river basins.

Funding:
The research is supported by the National Youth Thousand Talents Project (Y771071001).

Conflicts of Interest:
The authors declare no conflicts of interest.