Uncertainty of Hydrological Drought Characteristics with Copula Functions and Probability Distributions : A Case Study of Weihe River , China

This study investigates the sensitivity and uncertainty of hydrological droughts frequencies and severity in the Weihe Basin, China during 1960–2012, by using six commonly used univariate probability distributions and three Archimedean copulas to fit the marginal and joint distributions of drought characteristics. The Anderson-Darling method is used for testing the goodness-of-fit of the univariate model, and the Akaike information criterion (AIC) is applied to select the best distribution and copula functions. The results demonstrate that there is a very strong correlation between drought duration and drought severity in three stations. The drought return period varies depending on the selected marginal distributions and copula functions and, with an increase of the return period, the differences become larger. In addition, the estimated return periods (both co-occurrence and joint) from the best-fitted copulas are the closet to those from empirical distribution. Therefore, it is critical to select the appropriate marginal distribution and copula function to model the hydrological drought frequency and severity. The results of this study can not only help drought investigation to select a suitable probability distribution and copulas function, but are also useful for regional water resource management. However, a few limitations remain in this study, such as the assumption of stationary of runoff series.


Introduction
Hydrological drought refers to a lack of water in the hydrological system and gives rise to negative impacts on river ecosystems and human lives [1].Hydrological droughts are typically defined as periods when streamflow below a pre-defined threshold, called the threshold level method (TLM) [2].Advantages of the TLM are (i) no a priori knowledge of probability distributions is required, and (ii) it directly produces drought characteristics (e.g., duration, severity, frequency).When the variable of interest (x) (i.e., soil moisture, groundwater storage, or discharge) is below a predefined threshold (τ), a drought is assumed to have occurred.A constant or a threshold can be used, and because a variable threshold level takes seasonal patterns into account, it has been widely used [3,4].The threshold is usually assumed to be equal to a given percentile of the flow duration curve between the 30 percentile flow (Q 30 ) and the 5percentile flow (Q 5 ) [5], or the threshold can be obtained by fitting some kind of statistical function through the data (normal, gamma, beta, etc.) [6].The threshold is 25percentile Water 2017, 9, 334 2 of 15 (Q 25 ), meaning that a drought occurs when the streamflow is below the 25th percentile flow from the duration curve.The number of drought events obtained from the Q 30 , Q 25 , Q 20 , Q 15 , Q 10 , and Q 5 were compared with a previous study [7].Based on the yearbook of "the drought history in China ", which reported that the main drought years in the Weihe River are 1962, 1969, 1971-1972, 1977, 1982, 1986, 1994, 1995, and 1997, drought events defined by Q 25 are the closest to the historical drought conditions.Therefore, the threshold level in this study selected Q 25 to define drought events and examine hydrological drought characteristics in the Weihe River Basin.
Since hydrological droughts are complex events, one index cannot provide a comprehensive evaluation.A bivariate distribution is derived for describing the correlated hydrology drought characteristic variables.For instance, bivariate normal distribution [8], bivariate exponential distribution [9], and bivariate gamma distribution [10] are often applied to study flood problems.However, because the hydrology drought characteristic variables are highly correlated and may obey different marginal distribution functions, one of the drawbacks for these bivariate distributions is that the same family is needed for each marginal distribution [11].
Copulas are functions that join univariate distribution functions to form multivariate distribution functions [12].Sklar [13] introduced copulas, and they have been used in insurance and finance [14].Because of the flexibility of copulas, they became popular in hydrological analysis [15].Especially, Archimedean copulas, such as Clayton, Frank, and Gumbel-Hougard copulas, are the most popular family used in hydrology [16] and have been widely applied to model the dependence between hydrological variables [17].When using copula functions analysis of hydrological droughts, several marginal distribution functions have been commonly used to fit drought duration and drought severity, such as exponential (EXP), gamma (GAMA), lognormal (LOGN), general Pareto (GP), generalized extreme value (GEV), and Weibull (WBL) [18] distributions.The inversion of Kendall's method was usually used to estimate the copula parameter [19], and the goodness of fit was usually tested by the least root mean square error (RMES), Kolmogorov-Smirnov (KS) test, Anderson-Darling (AD) test, Akaike information criterion (AIC), ordinary least squares (OLS), and Bayesian information criterion (BIC) [20].
However, the selection of a threshold directly influences the duration and severity of drought events, and it is an uncertainty factor in modeling the bivariate distribution of hydrological droughts.Due to the changing environment and the impact from human activities, the variability of hydrological drought properties is also an important uncertainty factor.Since there are some goodness of fit test methods for choosing appropriate margin probability distributions and copula functions, the final selections may be different depending on the different candidate margins distributions and copulas.Therefore, it is necessary to investigate the uncertainty and sensitivity in modeling the joint distribution of hydrological droughts.In this paper, we choose six commonly-used univariate probability distributions as the candidate margins for drought duration and severity, and three Archimedean copulas are employed to match the joint distributions.
The remaining parts of this paper are organized as follows: Section 2 introduces the case study and the data from three hydrometric stations, the methodologies for defining the hydrological drought, six univariate models, and the theory of copulas and the formula of the return period; Section 3 shows the primary results, including the univariate and bivariate drought frequency analysis and the calculation of return periods; Section 4 discusses the sensitive and uncertainty of drought frequency analysis.;and, finally, the main conclusions are given in Section 5.

Weihe River Basin and Data
The Weihe River Basin, with a drainage area of 134,800 km 2 , is the largest tributary of the Yellow River in North China (Figure 1) (between 104 • -110 • E and 34 • -38 • N).This basin originates from the north of Niaoshu Mountain with an altitude of 3485 m above sea level.The most important Water 2017, 9, 334 3 of 15 topographic feature of the Weihe River Basin is the Loess Plateau in the north, which is the main source of sediments in the river [21].The annual average temperature ranges between 9.3 • C and 14.4 • C, the annual mean precipitation amounts are in the range of 558-750 mm with a general increasing trend from north to south, and the annual mean runoff amounts are 10.37 billion m 3 .The runoff from July to September accounts for about 60-70% of the annual discharge [22].Agricultural losses due to local drought disasters occupy over 50% of the total losses [23].
Water 2017, 9, 334 3 of 14 topographic feature of the Weihe River Basin is the Loess Plateau in the north, which is the main source of sediments in the river [21].The annual average temperature ranges between 9.3 °C and 14.4 °C, the annual mean precipitation amounts are in the range of 558-750 mm with a general increasing trend from north to south, and the annual mean runoff amounts are 10.37 billion m 3 .The runoff from July to September accounts for about 60-70% of the annual discharge [22].Agricultural losses due to local drought disasters occupy over 50% of the total losses [23].The hydrological stations of Linjiacun (LJC), Xianyang (XY) and Huaxian (HX) are located in the upper, middle, and lower streams for the main river channel (Figure 1).
The data of three stations are provided by the Institute of Soil and Water Conservation (http://loess.geodata.cn/),Chinese Academy of Sciences and Ministry of Water Resources.The monthly streamflow covers the period between 1960 and 2012 years.Linjiacun and Xianyang station lack data in 2004.

Defining Drought Duration and Severity
A drought event starts when the variable ( x ) is below the threshold level (onset; 1 t = ) and the event continues until the threshold is exceeded again (recovery; t T = ).Each drought event ( i ) can be characterized by its duration and by some measure of the severity of the event (Figure 2).The hydrological stations of Linjiacun (LJC), Xianyang (XY) and Huaxian (HX) are located in the upper, middle, and lower streams for the main river channel (Figure 1).
The data of three stations are provided by the Institute of Soil and Water Conservation (http: //loess.geodata.cn/),Chinese Academy of Sciences and Ministry of Water Resources.The monthly streamflow covers the period between 1960 and 2012 years.Linjiacun and Xianyang station lack data in 2004.

Defining Drought Duration and Severity
A drought event starts when the variable (x) is below the threshold level (onset; t = 1) and the event continues until the threshold is exceeded again (recovery; t = T).Each drought event (i) can be characterized by its duration and by some measure of the severity of the event (Figure 2).
Water 2017, 9, 334 3 of 14 topographic feature of the Weihe River Basin is the Loess Plateau in the north, which is the main source of sediments in the river [21].The annual average temperature ranges between 9.3 °C and 14.4 °C, the annual mean precipitation amounts are in the range of 558-750 mm with a general increasing trend from north to south, and the annual mean runoff amounts are 10.37 billion m 3 .The runoff from July to September accounts for about 60-70% of the annual discharge [22].Agricultural losses due to local drought disasters occupy over 50% of the total losses [23].The hydrological stations of Linjiacun (LJC), Xianyang (XY) and Huaxian (HX) are located in the upper, middle, and lower streams for the main river channel (Figure 1).
The data of three stations are provided by the Institute of Soil and Water Conservation (http://loess.geodata.cn/),Chinese Academy of Sciences and Ministry of Water Resources.The monthly streamflow covers the period between 1960 and 2012 years.Linjiacun and Xianyang station lack data in 2004.

Defining Drought Duration and Severity
A drought event starts when the variable ( x ) is below the threshold level (onset; 1 t = ) and the event continues until the threshold is exceeded again (recovery; t T = ).Each drought event ( i ) can be characterized by its duration and by some measure of the severity of the event (Figure 2).The duration of a drought event is calculated by [24]: where δ(t) is a binary variable indicating a drought situation with respect to time t, x(t) is the hydrological variable on time t, τ(t) is the threshold level of that hydrological variable with respect to time t, D i is the duration of drought event i, and ∆t is the time step of t.
For fluxes the most commonly-used severity measure is the deficit volume, calculated; by summing up the differences between the actual flux and the threshold level over the drought period [25].The equation is: where s(t) is the deviation with respect to time t, and S i is the deficit of drought event i.
The streamflow (m 3 /s) is converted into runoff depth (mm) by: where R is the runoff depth in mm; Q is the streamflow in m 3 /s, t is the monthly time in seconds, and A is the drainage area in km 2 .

Marginal Distribution Model and Copula-Based Models
Six commonly-used univariate probability distributions are selected as the candidate margins for drought duration and drought severity.They are the exponential (EXP), gamma (GAM), log-normal (LOGN), generalized Pareto (GP), generalized extreme value (GEV), and Weibull (WBL) distributions.The cumulative distribution functions (CDF) of six univariate distributions are given in Table 1.The maximum likelihood method is used to estimate the parameters.Akaike information criterion (AIC) [26] and the Anderson Darling (AD) [27] test are used to select the best univariate distribution.
Developed by Sklar [13], copulas are functions that link univariate distribution functions to form multivariate distribution functions, in which the domain is [0, 1].In terms of two random variables, Sklar's theorem, states that if F X,Y (x, y) is a two-dimensional joint cumulative distribution function, F X (x) and F Y (y) are marginal distribution functions of variables X, Y, then there exists a copula C such that: In which c is the density function of C, defined as: In this study, we choose widely-used three mono-parameter Archimedean copulas (Clayton, Frank, and Gumbel-Hougard) as the candidates for modeling the joint distribution of hydrological drought properties.Three copulas are defined as shown in Table 2.The inference function for margins Water 2017, 9, 334 5 of 15 (IFM), suggested by Joe [28], is employed to estimate the parameter in copulas.The ordinary least squares (OLS) and Akaike information criterion (AIC) are used for testing the goodness of fit.

Return Period of Droughts
The return period of a variable, defined as the average elapsed time between occurrences of an event with a certain magnitude, or greater [29].In this study, the return period of droughts can be defined by drought duration or drought severity as [30]: where T D is the return period of drought duration, T S is the return period of drought severity; N is the length of data series; n is the numbers of drought events, and F D (d) and F S (s) are cumulative distribution function of drought duration and severity, respectively There are two cases of bivariate return periods; one is the drought duration exceeding a specific value and the drought severity exceeding another specific value (D ≥ d and S ≥ s), called co-occurrence return period (T 0 ).Another case is the drought duration exceeding a specific value or the drought severity exceeding another specific value (D ≥ d or S ≥ s), called the joint return period (T α ).Both return periods, in terms of the copula-based bivariate drought distribution, are described below [31]: Water 2017, 9, 334 6 of 15

Drought Duration and Severity Characteristics
Based on the Q 25 threshold level of monthly runoff, there are a total of 46, 51, and 68 drought events, respectively, in three stations (Linjiacun, Xianyang, Huaxian) during 1960-2012 (Figure 3). Figure 4 shows the number of drought events of different drought duration and severity at three stations.The main statistical characteristic values of drought duration and drought severity are listed in Table 3.The most severe drought lasted for 17 months in Linjiacun and Xianyang stations.The average drought durations are 3.391, 3.059, and 2.382 months, respectively, at the three stations.This is why the number of drought events of LJC station (46) and XY station (51) are less than those of HX station (68).The skewness coefficient (which is a measure of the asymmetry of the probability distribution of a real-valued random variable about its mean) of drought severity is also larger than the drought duration.In view of the drought duration, the dispersion degree of Linjiacun station is the smallest.In view of the drought severity, the dispersion degree of Xianyang station is the smallest.The relationship between drought duration and drought severity is very good (Table 4); they reflect that there is very strong correlation between drought duration and drought severity in three stations.The Pearson linearly-dependent coefficients in the three stations are over 0.88.The Spearman rank correlation coefficient exceeded 0.77.For the Kendall rank correlation coefficient, its surpass 0.63.For example, the Pearson linearly-dependent coefficient ρ is 0.883, 0.916, and 0.885, respectively, in the three stations, which suggests the importance of using copulas based on the drought frequency analysis method.The relationship between drought duration and drought severity is very good (Table 4); they reflect that there is very strong correlation between drought duration and drought severity in three stations.The Pearson linearly-dependent coefficients in the three stations are over 0.88.The Spearman rank correlation coefficient exceeded 0.77.For the Kendall rank correlation coefficient, its surpass 0.63.For example, the Pearson linearly-dependent coefficient ρ is 0.883, 0.916, and 0.885, respectively, in the three stations, which suggests the importance of using copulas based on the drought frequency analysis method.

Selected Marginal Distributions
Table 5 lists the estimated parameters for six theoretical distribution models, and the results of AIC and AD test at the 99% (α = 0.01) significant level have shown in Table 6.According to AIC and the AD test, the best fitted marginal distributions are WBL and GEV, WBL and WBL, and WBL and GAM, at the three stations, for drought duration and severity, respectively.Even though the GP model has relatively low AIC, the existence of "outliers" resulting in the GP distribution is not acceptable by the AD test for drought duration and drought severity, the EXP model is also not acceptable by the AD test for the fit of drought severity at Huaxin station.Figure 5 compares the six fitted distributions with empirical distributions of the identified drought characteristic variables in the three stations.The estimated theoretical cumulative probabilities for the best-fitting distribution are quite close to the empirical ones, which denotes that these probability distributions perform fairly well.Due to the existence of a few very small value (drought severity are 0.000761 mm, 0.03164 mm, 0.04114 mm, and 0.082 mm, respectively) of drought severity at HX station, the cumulative distribution function of LOGN only goes up to 0.7.

Selected Copula Functions
The results of the parameters and goodness of fit of the candidate bivariate copula functions are shown in Table 7.According to OLS and AIC, the Frank copula is found to be the best-fitted one for Xianyang and Huaxian stations, and the Gumbel copula is the best-fitted for Linjiacun station.Figure 6 shows the probability-probability (PP) plot.It turns out that the estimated cumulative probabilities agree well with the empirical ones.

Selected Copula Functions
The results of the parameters and goodness of fit of the candidate bivariate copula functions are shown in Table 7.According to OLS and AIC, the Frank copula is found to be the best-fitted one for Xianyang and Huaxian stations, and the Gumbel copula is the best-fitted for Linjiacun station.Figure 6 shows the probability-probability (PP) plot.It turns out that the estimated cumulative probabilities agree well with the empirical ones.

Return Period of Droughts
According to Equations ( 8) and ( 9), return periods of 2, 5, 10, 20, 50, and 100 years, defined by separate drought duration and drought severity, are summarized in Table 8.For drought duration, the return period at Xianyag station is the largest, however, for drought severity, the return period of larger than 20 years for duration and severity at Linjiacun station is the largest.Considering two parameters, using joint probability distributions and copulas to calculate the co-occurrence return period and joint return period (in Figure 7), it is clear that the co-occurrence return period is greater than both of the return periods defined by the drought duration and drought severity separately, while the joint return period is less than both of the return periods.Under the same increasing range in univariate models, the increasing range of the co-occurrence return period is larger than the joint return period.This shows that these two kinds of combination return periods can be regarded as two extreme conditions of marginal distribution return periods.It is possible to estimate the interval of the actual return period according to the co-occurrence return period and the joint return period.

Return Period of Droughts
According to Equations ( 8) and ( 9), return periods of 2, 5, 10, 20, 50, and 100 years, defined by separate drought duration and drought severity, are summarized in Table 8.For drought duration, the return period at Xianyag station is the largest, however, for drought severity, the return period of larger than 20 years for duration and severity at Linjiacun station is the largest.Considering two parameters, using joint probability distributions and copulas to calculate the co-occurrence return period and joint return period (in Figure 7), it is clear that the co-occurrence return period is greater than both of the return periods defined by the drought duration and drought severity separately, while the joint return period is less than both of the return periods.Under the same increasing range in univariate models, the increasing range of the co-occurrence return period is larger than the joint return period.This shows that these two kinds of combination return periods can be regarded as two extreme conditions of marginal distribution return periods.It is possible to estimate the interval of the actual return period according to the co-occurrence return period and the joint return period.

Return Period of Droughts
According to Equations ( 8) and ( 9), return periods of 2, 5, 10, 20, 50, and 100 years, defined by separate drought duration and drought severity, are summarized in Table 8.For drought duration, the return period at Xianyag station is the largest, however, for drought severity, the return period of larger than 20 years for duration and severity at Linjiacun station is the largest.Considering two parameters, using joint probability distributions and copulas to calculate the co-occurrence return period and joint return period (in Figure 7), it is clear that the co-occurrence return period is greater than both of the return periods defined by the drought duration and drought severity separately, while the joint return period is less than both of the return periods.Under the same increasing range in univariate models, the increasing range of the co-occurrence return period is larger than the joint return period.This shows that these two kinds of combination return periods can be regarded as two extreme conditions of marginal distribution return periods.It is possible to estimate the interval of the actual return period according to the co-occurrence return period and the joint return period.The significance level of 99% (α = 0.01) is used to determine where the marginal distribution is to be rejected or accepted based on the AD test.The return periods of drought duration and severity are calculated based on six marginal distributions at three stations.Figure 8 shows their variability compared with the empirical return period.There are significant differences in the return period based on different marginal distributions.Considering drought duration, the return period of GP, GEV, and WBL models are best-fitted to the empirical value, the GEV and WBL models are fully accepted by the AD test, the value are minimal and, additionally, the GP model is rejected by the AD test.Considering drought severity, the return period of WBL at Xianyang station and GAM at Huaxian station are nearest to the empirical values, they are perfectly acceptable by AD test and, especially, the AIC values are minimal.The results indicate that the drought return period is sensitive to the selection of the margins.

Effects of the Selection of Margin Distributions to the Return Period
The significance level of 99% ( 01 .0 = α ) is used to determine where the marginal distribution is to be rejected or accepted based on the AD test.The return periods of drought duration and severity are calculated based on six marginal distributions at three stations.Figure 8 shows their variability compared with the empirical return period.There are significant differences in the return period based on different marginal distributions.Considering drought duration, the return period of GP, GEV, and WBL models are best-fitted to the empirical value, the GEV and WBL models are fully accepted by the AD test, the AIC value are minimal and, additionally, the GP model is rejected by the AD test.Considering drought severity, the return period of WBL at Xianyang station and GAM at Huaxian station are nearest to the empirical values, they are perfectly acceptable by AD test and, especially, the AIC values are minimal.The results indicate that the drought return period is sensitive to the selection of the margins.

Effects of the Selection of Copula Functions to the Return Period
In order to analyze the difference between three copula functions, there are three drought samples (with a duration of 6.55 months and a severity of 1.04 mm (denoted by drought event 1), a duration of 9.27 months and a severity of 4.36 mm (denoted by drought event 2), and a duration of 10.90 months and a severity of 7.71 mm (denoted by drought event 3)) at Huaxian station.The cooccurrence return period and the joint return period based on 36 possible joint distributions of duration and severity, and with different Copula functions, are calculated for the three selected

Effects of the Selection of Copula Functions to the Return Period
In order to analyze the difference between three copula functions, there are three drought samples (with a duration of 6.55 months and a severity of 1.04 mm (denoted by drought event 1), a duration of 9.27 months and a severity of 4.36 mm (denoted by drought event 2), and a duration of 10.90 months and a severity of 7.71 mm (denoted by drought event 3)) at Huaxian station.The co-occurrence return period and the joint return period based on 36 possible joint distributions of duration and severity, and with different Copula functions, are calculated for the three selected drought events at Huaxian station.Figure 9 are the boxplots and show their variability compared with the empirical values.For hydrological drought event 2, the co-occurrence return period ranges from 6.56 to 31.49years for the Clayton copula, 5.71 to 11.41 years for the Frank copula, and 5.45 to 32.07 years for the Gumbel copula.There are large differences for the three copula functions.At Huaxian station, the best-fitted copula selected by marginal distributions GEV and GAM is the Frank copula, the co-occurrence return period of the best marginal distribution and best copula is 7.15 years, and the co-occurrence return period of the empirical value is 7.67 years.This indicates that the value of Frank, which is selected as the best copula function at Huaxian station, is nearest to the empirical values.In addition, the mean value of the Frank copula for 36 possible joint distributions are also near the empirical values.The same result is in other two hydrological drought events, and the joint return period also shows the same situations.
Water 2017, 9, 334 12 of 14 drought events at Huaxian station.Figure 9 are the boxplots and show their variability compared with the empirical values.For hydrological drought event 2, the co-occurrence return period ranges from 6.56 to 31.49years for the Clayton copula, 5.71 to 11.41 years for the Frank copula, and 5.45 to 32.07 years for the Gumbel copula.There are large differences for the three copula functions.At Huaxian station, the best-fitted copula selected by marginal distributions GEV and GAM is the Frank copula, the co-occurrence return period of the best marginal distribution and best copula is 7.15 years, and the co-occurrence return period of the empirical value is 7.67 years.This indicates that the value of Frank, which is selected as the best copula function at Huaxian station, is nearest to the empirical values.In addition, the mean value of the Frank copula for 36 possible joint distributions are also near the empirical values.The same result is found in other two hydrological drought events, and the joint return period also shows the same situations.

Effects of Human Activities on Drought Frequency
The tendency of average annual runoff is decreasing at the three stations, especially, due to human activities, the observed runoff has significantly decreased in the period of 1991-2013 [7], which also shows that after 1991 human activities mainly resulted in the short-term drought fluctuation.In this study, because the main aim is to analyze the uncertainty of marginal distributions and copula functions on hydrological drought, we did not spilt the time series into two periods based on the turning point year.In addition, human activities affect the drought frequency in many ways.Irrigation area expansion means subtractions of water have greatly increased, which have resulted in the decrease of observed runoff [32].The changing of land use has increased water demand, which has aggravated hydrological drought.
In addition, there are other uncertain factors, such as the selected threshold level of streamflow, climate change, and the length of time series.The threshold usually varies from 30 Q to 5 Q of streamflow, and drought events of higher severity are "nested" inside drought events of lower severity, such as a 10 Q threshold drought implying the occurrence of 20 Q and 30 Q threshold drought events.Under the changing climate, the precipitation series and runoff series are nonstationary; the surface runoff is not only sensitive to precipitation, but also sensitive to temperature [33].The length of the times series determines the number of drought events.These uncertainty

Effects of Human Activities on Drought Frequency
The tendency of average annual runoff is decreasing at the three stations, especially, due to human activities, the observed runoff has significantly decreased in the period of 1991-2013 [7], which also shows that after 1991 human activities mainly resulted in the short-term drought fluctuation.In this study, because the main aim is to analyze the uncertainty of marginal distributions and copula functions on hydrological drought, we did not spilt the time series into two periods based on the turning point year.In addition, human activities affect the drought frequency in many ways.Irrigation area expansion means subtractions of water have greatly increased, which have resulted in the decrease of observed runoff [32].The changing of land use has increased water demand, which has aggravated hydrological drought.
In addition, there are other uncertain factors, such as the selected threshold level of streamflow, climate change, and the length of time series.The threshold usually varies from Q 30 to Q 5 of streamflow, and drought events of higher severity are "nested" inside drought events of lower severity, such as a Q 10 threshold drought implying the occurrence of Q 20 and Q 30 threshold drought events.Under the changing climate, the precipitation series and runoff series are non-stationary; the surface runoff is not only sensitive to precipitation, but also sensitive to temperature [33].The length of the times series determines the number of drought events.These uncertainty factors are not key points in this study, as we will discuss their effects on hydrological drought frequency in our next study.

Conclusions
This study investigated the regional drought frequency analysis in the Weihe River Basin considering the spatio-temporal structure of drought with copula functions.The primary conclusions are given as follows: (1) There are more drought events at Huaxian station (lower basin) than at Linjiacun station (upper basin), but there are longer drought durations and greater severity at Linjiacun station.(2) Based on the AD test, five models (EXP, GAM, LOGN, GEV, WBL) are acceptable for the fit of the drought duration and drought severity at LJC and XY stations.The GP model is not acceptable for the goodness-of-fit of drought duration and drought severity at three stations, the EXP model is rejected for the fit of drought severity at Huaxian station.(3) Based on ordinary least squares (OLS) and Akaike information criterion (AIC), the Frank copula is the best joint distribution function at Linjiacun and Huaxian stations, while the Clayton copula is the best-fitted model at Huaxian station.(4) The co-occurrence return period is greater than both the return periods defined by drought duration and drought severity separately, while the joint return period is shorter than both of the return periods.This shows that these two kinds of combination return periods can be regarded as two extreme conditions of the marginal distribution return period.It is possible to estimate the interval of the actual return period according to the co-occurrence return period and the joint return period.(5) The drought return period is sensitive to the selected marginal distribution and different copula functions.Therefore, it is important to select proper marginal distributions and copula functions, and the sensitivity and uncertainty of hydrological droughts should be paid more attention on the modeling and designing of drought models with consideration to the condition of water resources and the requirement of water management.

Figure 1 .
Figure 1.The study basin map.

Figure 2 .
Figure 2. Threshold level method with a variable threshold to define the drought duration and drought severity.

Figure 1 .
Figure 1.The study basin map.

Figure 1 .
Figure 1.The study basin map.

Figure 2 .
Figure 2. Threshold level method with a variable threshold to define the drought duration and drought severity.

Figure 2 .
Figure 2. Threshold level method with a variable threshold to define the drought duration and drought severity.

Figure 4 .
Figure 4.The number events of different drought (a) durations and (b) severities at three stations.

Table 4 .
Correlation between the hydrological drought characteristics' indices.indicates the Pearson linearly-dependent coefficient, τ indicates the Kendall rank correlation coefficient; s ρ indicates the Spearman rank correlation coefficient.

Figure 4 .
Figure 4.The number events of different drought (a) durations and (b) severities at three stations.

Figure 9 .
Figure 9. Sensitivity of the co-occurrence return period and joint return period of three drought events at Huaxian station and the selection of the univariate distribution and the copulas.(a): co-occurrence return period of drought event 1; (b): co-occurrence return period of drought event 2; (c): cooccurrence return period of drought event 3; (d): joint return period of drought event 1; (e): joint return period of drought event 2; (f): joint return period of drought event 3.

Figure 9 .
Figure 9. Sensitivity of the co-occurrence return period and joint return period of three drought events at Huaxian station and the selection of the univariate distribution and the copulas.(a): co-occurrence return period of drought event 1; (b): co-occurrence return period of drought event 2; (c): co-occurrence return period of drought event 3; (d): joint return period of drought event 1; (e): joint return period of drought event 2; (f): joint return period of drought event 3.

Table 3 .
Hydrological drought index of single-feature statistical characteristic values.

Table 4 .
Correlation between the hydrological drought characteristics' indices.
Note: ρ indicates the Pearson linearly-dependent coefficient, τ indicates the Kendall rank correlation coefficient; ρ s indicates the Spearman rank correlation coefficient.

Table 5 .
The estimated parameters by the maximum likelihood method.

Table 6 .
The goodness of fit of the marginal distribution.

Table 7 .
Parameters and the goodness of fit of the bivariate distributions based on copulas.

Table 7 .
Parameters and the goodness of fit of the bivariate distributions based on copulas.

Table 8 .
Return periods defined by drought duration and severity in three stations.Effects of the Selection of Margin Distributions to the Return Period

Table 8 .
Return periods defined by drought duration and severity in three stations.