Copula-based drought severity-area-frequency curve and its uncertainty, a case study of Heihe River basin, China

Drought is a stochastic natural phenomenon that arises from considerable deficiency in precipitation. With the global temperature increasing, it becomes more common and has great effects on terrestrial ecosystems (Mishra & Singh 2009; Fang et al. 2018a; Han et al. 2018). The frequency analysis of drought is important since it is capable of quantifying the drought characteristics (Rahmat et al. 2015). Bivariate or even multivariate frequency analysis is more informative than univariate analysis due to the complexity and multivariate features of drought (Yu et al. 2014; Jiang et al. 2015; Rahmat et al. 2015; Amirataee et al. 2018).

Previous literature has been most concerned with the severity and the duration of drought due to their ability to characterize drought from the temporal scale (e.g., Lee & Kim 2011; Beguería et al., 2014; Rahmat et al. 2015; Tan et al. 2015; Afshar et al. 2016; Homdee et al. 2016; Fang et al. 2018b; Hameed et al. 2018). However, the spatial extent of drought is also very important in drought assessments. Recent researches have shown that the drought area has been significantly increasing globally and regionally during the last several decades. The percentage of areas under drought had increased by about 1.74% per decade from 1950 to 2008 globally (Dai 2011). The drought-affected area increased nearly 12-fold from 1950 to 2009 in China (Wang et al. 2012). A significant increase in drought area has been detected in the northwest of China, with an increasing rate of 3.96% per decade since the late 1990s (Yu et al. 2014). Therefore, taking both severity and areal extent into account is very necessary for better understanding drought from both a temporal and spatial scale, and it would be valuable for drought-based decision-making.

Applying the concept of quantifying the frequency analyses of both severity and areal extent based on the copula theory, the drought severity-area-frequency (SAF) curves can be derived. These curves are very useful in drought planning and management since they can provide the probabilities and return periods of the drought areal extent with different severities over a region (Reddy & Ganguli 2013; Amirataee et al. 2018). However, uncertainties are usually involved in the estimation of drought SAF curves due to the effects of different copula functions, marginal distributions, input data, parameter estimations, etc.

The former uncertainty sources, copula functions and marginal distributions, are of major concern in previous literatures. For example, Halwatura et al. (2016) compared four marginal distribution functions and two types of copulas, and assessed the impacts of these uncertainties on the estimations of recurrence intervals (RIs) of drought events. They stated that no differences were found between the RIs derived from the different marginal distributions for mild drought events, and the effects of different copulas on RIs can nearly be negligible in their study. However, Zhao et al. (2017) drew a slightly different conclusion by comparing six marginal distributions and three Archimedean copulas in a case study of the Weihe River basin. They pointed out that the return period of drought varies depending on the selected marginal distributions and copula functions, and the differences become larger with the increase of return periods. Zhang et al. (2015) evaluated the influences of marginal distributions to the uncertainty of joint distribution by employing a Bayesian framework, and the results showed that the stronger the tail of marginal distribution, the greater the uncertainty of joint distribution, especially for extreme events. All the above studies made great contributions in enriching the uncertainty analysis in drought frequency curves. However, the uncertainties caused by copula parameter estimations and copula input data are rarely regarded. Sadegh et al. (2017) presented a multivariate copula analysis toolbox which has the ability to quantify the uncertainties in drought frequency curves associated with copula parameter estimations and copula functions based on Bayesian framework; whereas the fixed drop-down menu options in the toolbox and the non-editable display for the outputs restrict its wide application to some extent.

This paper aims to develop a drought SAF curve with copula theory and to evaluate the uncertainties in the curve caused by different copula functions, copula parameter estimations, and copula input data. The following issues will be addressed, namely: (i) to compare and select the copula function with best performance in modeling the dependence structure of drought severity and drought areal extent; (ii) to construct a drought SAF curve based on the conditional probability distribution in the study area; and (iii) to quantify the uncertainties in drought SAF curves from three uncertainty sources.

A number of indices have been developed so far for drought quantification, monitoring, and analysis, such as PDSI (Palmer Drought Severity Index), SPI (Standardized Precipitation Index), and SPEI (Standardized Precipitation-Evapotranspiration Index), etc. Under the current climate warming conditions, many studies show that, SPEI, developed by Vicente-Serrano et al. (2010), could give a better performance in drought assessments due to the fact that it combines the multi-scalar character with the capacity of involvement of temperature effects on droughts (e.g., Tan et al. 2015; Homdee et al. 2016). Thus, SPEI was selected and 12-monthly scale of SPEI will be analyzed in this study. The temperature-based Thornthwaite method is used to calculate potential evapotranspiration (PET) as it only requires monthly mean temperature data and the latitudinal coordinate of the location, which are readily available at most meteorological stations. SPEI is computed with the assumption that P-PET (precipitation minus potential evapotranspiration) series follow the Log-logistic distributions. Detailed information about SPEI and its calculation can be found in the studies of Vicente-Serrano et al. (2010), Beguería et al. (2014), and Hameed et al. (2018).

Different copulas differ in describing the dependence structures. Comparisons of copulas could be necessary before bivariate drought frequency analysis. Ten copulas are selected and given in Table 1, including not only one-parameter Archimedean families such as Clayton, Gumbel, Frank, and Joe, but also two-parameter copulas, BB1 (Clayton–Gumbel), BB6 (Joe–Gumbel), BB7 (Joe–Clayton), and BB8 (Joe–Frank). The more flexible structures of two-parameter copulas allow for different non-zero lower and upper tail dependence coefficients (Brechmann & Schepsmeier 2013). Many of the above copulas are commonly used in hydrological and climatological applications (e.g., Xu et al. 2015; Zhang et al. 2015, 2017; Zhao et al. 2017; Ayantobo et al. 2018).

The copula parameter (or and ⁠) in Table 1 helps to measure the extent of relationship between variables (Ayantobo et al. 2018). Its estimation from sample data is not straightforward. One of the commonly used approaches is maximum likelihood estimation (MLE) (Brechmann & Schepsmeier 2013). Although the solution to MLE involves the roots of a cubic equation, and cannot be written in a simple form, the MLE can be easily found using numerical routines. The standard errors for MLE copula parameter estimations are based on inversion of the Hessian matrix (Brechmann & Schepsmeier 2013).

According to MLE, the optimal value of copula parameter can be derived, together with its standard errors which can be used to represent the parameter uncertainty ranges (without considering the uncertainty in marginal distributions in severity and areal extent). The effects of copula parameter estimations to drought SAF curves are quantified by random generating copula parameters based on the asymptotic normality first, then calculating the drought curves according to the new generating copula parameter values, and finally identifying the and th percentiles of drought curves which are defined as the lower and upper limits of confidence intervals with given probabilities. These confidence intervals are the uncertainty ranges of drought curves caused by copula parameter estimations.

To quantify the effects of copula input data uncertainty to drought SAF curves, the non-parametric bootstrap method, which has usually been used in quantifying the sample uncertainty in previous studies (e.g., Hu et al. 2015; Vergni et al. 2017), is selected. The idea behind this method is that since the distribution from which the sample has been taken is unknown, the values in that sample are the best guide to the true distribution. The theoretical details of the bootstrap approach can be found in Efron (1979). The copula input data in this case consists of n sets of drought severity (S) and drought areal extent (A), denoted as (S₁, A₁), (S₂, A₂), …, (S_i, A_i), …, (S_n, A_n), in which n means the data length. The S_i and A_i series are fitted by theoretical probability distributions and the best distribution is selected based on goodness of fit tests. The specifics of bootstrap procedure are as follows:

1. In the jth bootstrap, resampling the original n set values of drought variables with replacement, and generating a bootstrapped sample of the same size n as ((S₁^j*, A₁^j*), (S₂^j*, A₂^j*), ……, (S_i^j*, A_i^j*), ……, (S_n^j*, A_n^j*).

2. Using the probability distribution type which has been pre-specified according to the original severity and areal extent series, to fit the S^j* and A^j* for obtaining their optimal parameter estimates and the corresponding marginal distributions, CDF(S^j*) and CDF(A^j*).

3. Using CDF(S^j*), CDF(A^j*) and the best candidate copula pre-selected according to the original marginal distributions, to estimate the copula parameters (i = 1, or i = 1, 2).

4. Calculating drought severities at given drought areal extents according to the conditional probability functions, and then obtaining the drought curve SAF^j*.

5. The above procedure is repeated m times to generate m bootstrap samples, m estimations of marginal distributions for severity and areal extent, m sets of copula parameters ⁠, and m drought curves SAF^j* (j = 1, 2, …, m).

6. From drought curves SAF^j* (j = 1, …, m), identifying the and th percentiles of drought curves which is defined as the lower and upper limits of confidence intervals with given probabilities, and then analyzing the contributions of copula input data to the uncertainty of drought curves.

The main steps of this study are as follows. The values of 12-month SPEI at the stations are calculated first. Then, the drought severity and drought areal extent are determined according to IDW spatial interpolation technique. The two drought variables are fitted by a variety of theoretical probability distributions. After that, copulas are used to compare the performances of capturing the dependence structures of the two drought variables with goodness of fit tests. Drought SAF curves are then extracted with the best copula. Finally, the uncertainties in SAF curves caused by different copulas, copula parameter estimations, and copula input data, are investigated.

The Heihe River basin, which falls between 97.1°–102.0°E and 37.7°–42.7°N, is the second largest typical inland river basin in the northwest of China with a controlled catchment area of 143,000 km² (Figure 1). It is one of the most drought-prone river basins in China (Wang et al. 2016; Ayantobo et al. 2018; Niu et al. 2019) with an annual precipitation of approximately 150 mm across the basin, an annual average air temperature of 6–8 °C, and annual potential evaporation of 1,000–2,000 mm. The runoff of the basin mainly comes from the upper reach, and consumed by large areas of irrigation in the middle reach and oasis ecosystem in the lower reach. Severe drought events would lead to inadequate water supply and thus affect the local agriculture and vulnerable oasis ecosystem. In this study, monthly precipitation and temperature data spanning from 1960 to 2015 were collected from nine meteorological stations within the basin. The locations and the detailed information of stations are shown in Figure 1 and Table 2.

Drought severity and areal extent series over the study area determined by Equations (1) and (2) are plotted in Figure 2. It is clear that, before 1997, the drought was mainly characterized by light drought conditions. However, several moderate or even severe droughts were captured after 1997, and furthermore, drought events with stronger severity and larger areal extent occurred more frequently than before.

It is necessary to evaluate the dependency between variables prior to the application of copula when analyzing the bivariate joint probability distribution of such random variables (Amirataee et al. 2018). For this purpose, Pearson's, Spearman's, and Kendall's correlation coefficients are applied to explore the dependency between drought severity and drought areal extent. These correlation coefficients are often used in measuring the monotone association. The first belongs to the parametric method and sensitive to outliers, and the latter two belong to the non-parametric methods and usually suggested for non-normally distributed data (Shong 2010).

The correlation coefficients between the two drought variables are 0.678, 0.722, and 0.549, respectively, and all the corresponding p-values are less than the significance level of 0.05. It is concluded that there is a significant correlation between the variables of drought severity and drought areal extent at 95% confidence level.

The first step in copula fitting consists of the identification of appropriate marginal distributions for the two drought variables. Since they are typically unknown, we first fit them by using several theoretical distributions. According to previous studies (e.g., Loukas & Vasiliades 2004; Serinaldi et al. 2009; Reddy & Ganguli 2013; Afshar et al. 2016; Amirataee et al. 2018), ten univariate probability distributions (Exponential, Normal, Lognormal, Logistic, Loglogistic, Gamma, Weibull, Log-Pearson III, Burr, and GEV distribution) are tested for fitting drought severity (S). Since the drought areal extent has definite ranges from 0 to 1, the bounded distributions, Beta, Power function, Pert, and Uniform, are tested as its candidate margins. Amirataee et al. (2018) stated that Beta distribution performed satisfactorily in fitting the drought area extent in their study. Kolmogorov–Smirnov (K-S), Anderson–Darling (A-D), and Chi-squared (C-S) are used to test the goodness of fit.

Results are shown in Table 3. The critical values for the three tests are 0.061, 2.502, and 15.507 at 0.05 significance level, respectively. As shown, GEV is selected finally as the marginal distribution for drought severity, because it not only passes the K-S, A-D, and C-S tests at 0.05 significance levels (all the statistic values are lower than the critical values for the three tests), but shows the best performance with the smallest statistic values. This distribution was also proved to be suitable in fitting the drought severity in many other studies (e.g., Lee & Kim 2011; Reddy & Ganguli 2013). As to the variable of drought areal extent, none of the distributions passes the tests. Then, we perform the fitting process based on its pseudo observations (Ghoudi & Rémillard 2004). In light of the references of Ghoudi & Rémillard (2004) and Kojadinovic et al. (2010), it is acceptable to replace the unknown marginal by their empirical counterparts obtained from the pseudo observations.

The results of modeling the joint probability distribution of the two drought variables by using ten alternative copula functions are summarized in Table 4, including the statistic values of goodness of fit tests (AIC, BIC, and log-likelihood), the optimal parameter estimations, and their standard errors. Clearly, nearly all the copulas have close AIC, BIC, and log-likelihood values except for Frank which shows the greatest departure from the others, indicating that most of the selected copulas have similar performance in capturing the dependence structures of drought severity and drought areal extent in this study case. Comparatively, Gumbel, among the one-parameter copula functions, provides the best fit with the smallest AIC (−473.45) and BIC (−469.24) values. Among the two-parameter copulas, BB7 performs the best. Student-t and BB8 copulas even perform inferior to Gumbel according to the corresponding AIC and BIC values.

To better understand the performance of each copula in modeling the drought variables in this study, the generated samples derived from the estimated theoretical copula by Monte Carlo technique are compared with the observed data for checking whether they have similar features. A set of 5,000 samples are generated, and the scatter plot of generated samples, together with the observed datasets are shown in Figure 3. Different copula functions perform a little differently. The Normal, Student-t, and Frank are symmetric and can be used to model symmetric correlation patterns of variables. Gumbel and Clayton are asymmetric and can model the asymmetrical tail correlation patterns. Gumbel exhibits greater correlation in the upper tail, whereas Clayton exhibits greater correlation in the lower tail (Aas et al. 2009).

From the scatterplot, the observed data in the two opposite corners appear clustering, seemingly displaying both upper- and lower-tail dependence, while the dependence structures in the tails appear slightly asymmetric and a little skewed. Therefore, asymmetric copulas may be a good choice. Combined with the results of goodness of fit presented in Table 3, Gumbel is then selected as the benchmark copula in this study. The two-parameter copulas, BB1, BB6, and BB7, are not selected herein since, on one hand, they show similar performance with Gumbel and, on the other hand, the higher model complexity (with more parameters) would make a stimulation over conditioning of the model to some extent (Sadegh et al. 2017).

Actually, the optimality of copula functions could be dataset- and case-specific (Afshar et al. 2016). For example, Reddy & Ganguli (2013) evaluated three copula families in modeling the joint dependence between the drought severity and drought areal extent in western Rajasthan of India, and found that the Gumbel–Hougaard copula best represents the drought properties in their study. Xu et al. (2015) constructed a trivariate joint distribution based on the copula theory and concluded that the Joe and Gumbel copulas are more suitable to estimate the joint distribution of drought duration, affected area, and severity in southwest China. Amirataee et al. (2018) regarded the Frank copula as the most appropriate in modeling the joint probability of drought severity and percent of area in the Lake Urmia basin, Iran. In our case, with considering the results of goodness of fit and the greater dependence in upper tails of the data, the one-parameter Gumbel is finally selected as the benchmark copula.

The drought SAF curve, developed by using conditional probability function based on Gumbel copula, is plotted in Figure 4. The horizontal axis represents the percentage of drought areal extent and the vertical axis represents the drought severity. Only droughts with conditional probability of less than 0.2 (corresponds to return period of more than five years) are considered in the analysis. It can be seen that the drought event with larger areal extent is generally associated with higher severity. For example, with a conditional probability of 0.1 (at the return period of ten years), when the drought areal extent amounts to 20%, the average drought severity is 0.98, belonging to light drought; while when the drought areal extent increases to 60% and 80%, the severity increases to 1.17 and 1.30, resulting in moderate droughts.

The drought SAF curve based on the conditional probability can allow a water planner to assess the probability of a drought with a certain areal extent and severity occurring simultaneously. For example, for a drought areal extent reaching 40%, the occurrence probability is 0.05 when the severity is 1.12, and 0.02 when the severity increases to 1.21. This information cannot be obtained from the univariate frequency analysis of droughts, in which only the probability of univariate occurrence can be derived. The copula-based drought SAF curve can provide decision-makers with probability maps of drought severity and useful information on drought areal extent in the study area, yielding important information for water management purposes.

In this subsection, whether different copulas lead to significantly different assessments of SAF curves is investigated. Figure 5 shows the conditional probability curves derived from different copula functions. Conditional probabilities of 0.005, 0.01, 0.02, 0.05, 0.1, and 0.2 are compared. Gumbel copula is selected as the benchmark as stated above.

With a conditional probability of 0.2, most of the curves are close to each other, and only slight differences are found between the benchmark and the alternatives, except for the two curves derived from Clayton and Frank. As the conditional probability decreases (conditional return period increases), these differences increase in general, being consistent with the findings of Zhao et al. (2017), who pointed out that with an increase of the return period, the drought differences from different copula functions become larger. Such differences are not very large in this study as a whole, except for those from Clayton and Frank.

As to Clayton and Frank, large differences are detected. Interestingly, it is also noted that the AIC and BIC values for Clayton are −442.62 and −438.41 (Table 3), very close to those for Gumbel (−473.45 and −469.24). However, as shown in Figure 3, the drought SAF curves derived from these two copulas have a considerably larger difference, and Clayton gives larger possibilities of drought severity than Gumbel.

BB6 gives quite close drought severity estimations with Gumbel no matter how much the conditional probability is. BB1, BB7, and Student-t copulas give similar curve shapes with Gumbel, whereas they overestimate the drought severity overall. A similar conclusion can be drawn from Normal when the drought areal extent is higher than 10%, overestimating the drought severities comparatively. Joe and BB8 overestimate the drought severity when the drought area is lower than about 10%, but underestimate the drought severity when the drought areal extent becomes higher.

In summary, the drought SAF curves are influenced by copula choices, and the curves have different conditional probabilities of occurrence depending on which copula is used. BB6 gives quite close estimations in drought severity to Gumbel. Normal, Student-t, BB1, BB7 overestimate the drought severities, while Joe and BB8 underestimate the drought severities as a whole. Clayton and Frank copulas depart the most, and overestimate the drought severity much more. As the conditional probability decreases, the differences in drought severity estimations from different copulas increase, although most of the differences are not large. It implies that the choice of copula should be implemented with caution since different copulas show different performances in tail behaviors, especially when the focus is on extreme events (with small probabilities).

After estimating the copula parameter(s) with the MLE method, 1,000 copula parameters are randomly generated from the Normal distribution with the mean and standard error derived from MLE. The estimated optimal parameter value and standard error for Gumbel copula is 2.27 and 0.114 (see Table 3). Such range of standard error in parameter estimation is then translated to drought SAF curves, rendering the analysis of uncertainty. Combining with the conditional probability function and the generated copula parameters, a set of conditional probability isolines are then derived, which will be used to analyze the SAF curve uncertainties caused by copula parameter estimations.

Figure 6 presents the 95% CIs of conditional probability isolines associated with the uncertain parameter estimations of Gumbel. To avoid the curve overlapping, the six conditional probabilities are separated into two plots. Obviously, the uncertainty in copula parameter estimation has a certain effect on the conditional probability isolines. The uncertainty ranges of drought SAF curves become wider with the conditional probability decreasing. This is because the lower the probability of the drought event, the less the number of the observed data, and the less information available for the analysis.

Kendall's tau method is also used to estimate the parameter in Gumbel, yielding an optimal parameter value of 2.22 with a standard error of 0.130, which does not deviate much from the results of the MLE method. The corresponding uncertainty ranges of conditional probability isolines are also quite close to the results from the MLE method.

Using bootstrap technique to generate 1,000 set of drought samples ((S₁^j*, ₁^j*), (⁠₂^j*, ₂^j*), …, (⁠_i^j*, _i^j*), …, (⁠₅₀₁^j*, ₅₀₁^j*) (⁠ = 1, ⁠, 1,000)). To help visual comparisons of dependence structures among the original and the generated copula input data, the marginal distributions of the original samples and bootstrapped samples of drought severity and drought areal extent are plotted in Figure 7. The density of points in any area represents the density of the copula. It is observed that the bootstrapped samples well preserve the dependence structure between the drought severity and drought areal extent.

Uncertainty bands in Figure 8 reveal the effects of copula input data uncertainty to the estimation of conditional drought severity. Obviously, with the conditional probability decreasing, the uncertainty bands are increasing, similar findings as those from the effects of copula parameter uncertainty. Table 5 presents the lower (2.5% percentile) and upper (97.5% percentile) limits and the ranges of estimated drought severity at a certain drought areal extent owing to both input data uncertainty and parameter estimation uncertainty. The ranges resulting from both uncertainty sources grow remarkably as the conditional probability decreases. Compared with those caused by copula parameter, the uncertainty bands caused by copula input data present more widely. This reflects that the uncertainty of input data brings more uncertainty to drought SAF curves than that of copula parameter estimations. In other words, to improve the quality of copula input data can reduce the uncertainty of drought SAF curves more effectively.

In this study, 12-month SPEI is calculated over the Heihe River basin in northwest China, the drought SAF curve is developed, and then the uncertainties in the curve associated with different copula functions, copula parameter estimations, and copula input data are investigated.

Results show that GEV is finally selected as the marginal distribution for drought severity, and pseudo observations for drought areal extent. Gumbel shows the best performance in capturing the dependence structure of drought severity and drought areal extent and thus is selected as the benchmark copula. The estimations of drought SAF curves are influenced by copula choices. BB6 is capable of giving close drought severity as Gumbel; Normal, Student-t, BB1, BB7 overestimate the drought severities; Joe and BB8 underestimate the drought severities when the drought areal extents become higher; Clayton and Frank produce the most deviation in drought SAF curves although Clayton even presents quite close AIC and BIC values with Gumbel. The differences in drought SAF curves caused by different copulas are increasing with the conditional probability decreasing. As to the parameter estimations and input data, both of them result in the drought SAF curves uncertainty. The uncertainty ranges caused by copula input data are larger than those caused by copula parameter estimations, and consequently, to improve copula input data quality is recommended for developing the reliability of such frequency analysis. Results also show that the uncertainty bands in the drought SAF curves increase as the conditional probability decreases. It indicates that larger uncertainties exist in the frequency analysis of extreme values (with small probability). This highlights the importance of uncertainty analysis in frequency study, given that most studies in hydrology and climatology focus on extreme values. Couple-based frequency analysis of drought with considering the uncertainty issues is capable of providing more beneficial probability information to water managers, and aiding in the development of mitigation strategies and disaster preparedness plans for different levels of drought.

The first author is grateful for financial support from the China Scholarship Council and China University of Geosciences (Beijing), as well as the support provided by CSIRO during her visit in Australia. This study is supported by the Fundamental Research Funds for the Central Universities (No. 35832015028) and NSFC (No. 41101038). The authors also thank the anonymous referees for their comments and suggestions that have led to the quality of the paper being improved.