Abstract

This study introduces an approach to copula model parameter estimation within a Bayesian framework. The model is applied to estimate the return periods of the 2013–2015 drought in the Han River watershed (South Korea) as observed through a network of 18 rainfall stations in the watershed. For the univariate return periods, some of the 2013–2015 drought durations (i.e., 9–24 months) have return periods of around 30 years while the drought severities (which range from about 400 to 1,400 mm) at certain stations might be regarded as extreme events with return periods in the hundreds of years. The recent drought is extraordinary at many stations in the Han River watershed although the associated uncertainty is quite large (e.g., for the Seoul station, a range of about 180 to 2,250 years). Moreover, it is clearly seen that the unprecedented joint return periods of the 2013–2015 drought at many stations are consistent with the very high drought severity. The joint return period of the drought duration and severity is above 100 years for many of the stations and upwards of 1,000 years for a few stations.

INTRODUCTION

A number of univariate drought frequency methods have been developed and applied to estimate return periods of droughts (Loganathan et al. 1985; Engeland et al. 2004; Arena et al. 2006; Cebrián & Abaurrea 2006; Lana et al. 2008; Sadri & Burn 2012). Extreme drought events are particularly of interest from a hydrologic design perspective. They are best represented considering the multivariate dependence across the duration, severity, and intensity of an event (Favre et al. 2004; Genest et al. 2007; Renard & Lang 2007; Klein et al. 2009).

Modeling the multivariate dependence has been recognized as necessary for hydrologic risk estimation (Shiau 2006; Hao & Singh 2012, 2013; De Michele et al. 2013; Madadgar & Moradkhani 2014; Liu et al. 2015), especially for assessing risks associated with climate change and variability (Lall et al. 2015). For drought risk analyses, hydrologists are interested in how joint return periods or joint likelihoods for drought duration and severity could be estimated to understand their mutual dependences on drought risk (Shiau 2006; Madadgar & Moradkhani 2011; Yoo et al. 2012; Gräler et al. 2013).

Several studies have been conducted to estimate the joint or conditional probability of extreme events with multivariate joint distributions (Raynal-Villasenor & Salas 1987; Yue et al. 2001; Shiau & Shen 2001; Favre et al. 2004; Salas et al. 2005; Shiau 2006; Renard & Lang 2007; Yoo et al. 2012; De Michele et al. 2013; Hao & Singh 2013). A major difficulty in estimating the multivariate joint density function for drought risk analyses is that the marginal distributions of the variables of interest (e.g., duration and severity) are often different. Copula-based multivariate frequency techniques have been successfully applied to drought frequency analysis as a solution of the issue of dissimilar marginal distributions (Shiau 2006; Madadgar & Moradkhani 2011; Ma et al. 2013; Zhang et al. 2013; Xu et al. 2015). Copula functions were initially introduced by Sklar (1959) and they provide the flexibility to identify and fit a dependence structure to variables that have different marginal distributions. Despite many successful applications of copula-based frequency analysis in the hydrology community, the issue of estimating the parameter uncertainty is not often addressed. There have been limited studies on the use of Bayesian inference for copula function in hydrology including for estimation of bivariate return periods (Parent et al. 2014; Zhang et al. 2015; Bracken et al. 2016; Kwon & Lall 2016; Kwon et al. 2016; Sarhadi et al. 2016).

Over the last decade, Korea has frequently experienced droughts. Furthermore, future droughts may be more severe due to climate change (Kim et al. 2011; Kwon et al. 2011; Lee et al. 2012; Yoo et al. 2012). Since 2013, South Korea's water supplies have been significantly stressed as a result of a peninsula-wide ongoing drought. The Han River watershed (i.e., Seoul, Gyeonggi, and Gangwon provinces) has been particularly affected by the recent drought. The region has received about half of its average annual rainfall (1,200 mm). In June 2015, the water level of the Soyang Dam (which is the largest dam in South Korea) was at the lowest level recorded (i.e., 152.25 m a.s.l.) since the dam was built in 1974. In these contexts, this study primarily uses a Bayesian parameter estimation approach to estimate the return periods and their uncertainties associated with the recent 2013–2015 drought events in South Korea.

The data and watershed are introduced immediately below. The theoretical framework for the modeling approach is presented in the section ‘Method and model formulation’. In the following two sections the model is applied to the synthetic data and the case study for the recent 2013–2015 drought in South Korea. The final section presents a summary and conclusions.

HAN RIVER WATERSHED AND DROUGHT DEFINITIONS

We consider the Han River watershed as a case study due to the importance of the watershed from a water security perspective (the watershed is the primary water supply for the Seoul metropolitan area). The total watershed area of the Han River watershed is approximately 35,770 km2, and the annual average discharge of the basin is about 613 m3/s. The monthly rainfall data of 18 weather stations for the Han River watershed were obtained from the Korea Meteorological Administration (KMA) (http://web.kma.go.kr/eng), spanning the period 1974–2015. The locations of the gauging stations used in this study are shown in Figure 1.

Figure 1

A map showing Han River watershed along with weather stations (filled circle).

Figure 1

A map showing Han River watershed along with weather stations (filled circle).

Meteorological drought can be defined as an extended period of deficient precipitation. In this study, the 6-month rainfall is considered to evaluate seasonal droughts that occur due to rainfall deficit in wet (or monsoon) season and dry season (or non-monsoon). The accumulated 6-month monthly rainfall was transformed to anomaly time series determined by subtracting the long-term mean of the monthly rainfall time series as written in Equation (1). The monthly long-term mean of the full record was calculated as follows: 
formula
(1)
where is the accumulated 6-month rainfall at the end of a month i and year n.

We considered continuous negative precipitation anomalies to mark a drought, and the first instance of positive precipitation anomalies to mark the end of a drought. In general, the consecutive drought episodes are assumed to be independent realizations. We tested the independence of the extracted consecutive drought episodes (or variables) using autocorrelation function. It was found that the drought variables are linearly independent. As a result of this simple definition, we can define a time series of drought conditions and extract four historic characteristics: drought duration, drought severity, drought intensity (i.e., drought severity/drought duration), and drought interval time. The drought severity is determined as the accumulated rainfall deficit over the corresponding duration. We extracted drought characteristics from the time series of precipitation anomalies for all the stations in Han River watershed (not all are shown here due to space limitations; please see Figure S1, available with the online version of this paper). Information on the gauging stations used in this study and the basic statistics for the drought characteristics are summarized in Table 1.

Table 1

Information on the gauging stations used in this study and the basic statistics for the drought characteristics

Station no. Station name Location (Lat., Lon.) Annual rainfall
 
Drought
 
Mean (mm) Standard deviation (mm) Max duration (month) Max severity (mm) 
90 Sokcho 38.15, 128.34 1,358.36 252.87 22 (19) 670.62 (655.61) 
100 Daeqwalyeong 37.41, 128.43 1,682.77 486.59 35 (20) 1,122.27 (1,068.64) 
101 Chuncheon 37.54, 127.44 1,304.07 312.62 22(22) 1,163.71(1,163.71) 
105 Gangneung 37.45, 128.53 1,420.05 308.60 34 (18) 687.47 (565.16) 
108 Seoul 37.34, 126.57 1,380.63 373.14 24(24) 1,205.18(1,205.18) 
112 Incheon 37.29, 126.37 1,179.16 292.91 22(22) 929.48(929.48) 
114 Wonju 37.20, 127.57 1,322.65 291.60 17 (9) 570.01 (406.57) 
119 Suwon 37.16, 126.59 1,283.14 272.62 22(22) 785.20(785.20) 
131 Chongju 36.38, 127.26 1,228.57 261.84 24 (8) 568.42 (470.30) 
201 Ganghwa 37.42, 126.27 1,308.09 334.11 24(24) 1,407.72(1,407.72) 
202 Yangpyeong 37.29, 127.30 1,364.01 343.72 22(22) 1,133.29(1,133.29) 
203 Icheon 37.16, 127.29 1,329.46 305.61 18 (9) 570.65 (492.96) 
211 Inje 38.04, 128.10 1,160.97 300.43 24(24) 899.67(899.67) 
212 Hongcheon 37.41, 127.53 1,327.09 338.32 24(24) 1,260.22(1,260.22) 
221 Jecheon 37.09, 128.11 1,346.76 348.62 16 (9) 578.26 (536.95) 
226 Boeun 36.29, 127.44 1,274.54 329.78 23 (9) 785.86 (476.89) 
272 Youngju 36.52, 128.31 1,266.32 294.94 18 (9) 634.07 (483.17) 
273 Mungyeong 36.37, 128.09 1,239.07 289.11 26 (9) 843.19 (499.42) 
Station no. Station name Location (Lat., Lon.) Annual rainfall
 
Drought
 
Mean (mm) Standard deviation (mm) Max duration (month) Max severity (mm) 
90 Sokcho 38.15, 128.34 1,358.36 252.87 22 (19) 670.62 (655.61) 
100 Daeqwalyeong 37.41, 128.43 1,682.77 486.59 35 (20) 1,122.27 (1,068.64) 
101 Chuncheon 37.54, 127.44 1,304.07 312.62 22(22) 1,163.71(1,163.71) 
105 Gangneung 37.45, 128.53 1,420.05 308.60 34 (18) 687.47 (565.16) 
108 Seoul 37.34, 126.57 1,380.63 373.14 24(24) 1,205.18(1,205.18) 
112 Incheon 37.29, 126.37 1,179.16 292.91 22(22) 929.48(929.48) 
114 Wonju 37.20, 127.57 1,322.65 291.60 17 (9) 570.01 (406.57) 
119 Suwon 37.16, 126.59 1,283.14 272.62 22(22) 785.20(785.20) 
131 Chongju 36.38, 127.26 1,228.57 261.84 24 (8) 568.42 (470.30) 
201 Ganghwa 37.42, 126.27 1,308.09 334.11 24(24) 1,407.72(1,407.72) 
202 Yangpyeong 37.29, 127.30 1,364.01 343.72 22(22) 1,133.29(1,133.29) 
203 Icheon 37.16, 127.29 1,329.46 305.61 18 (9) 570.65 (492.96) 
211 Inje 38.04, 128.10 1,160.97 300.43 24(24) 899.67(899.67) 
212 Hongcheon 37.41, 127.53 1,327.09 338.32 24(24) 1,260.22(1,260.22) 
221 Jecheon 37.09, 128.11 1,346.76 348.62 16 (9) 578.26 (536.95) 
226 Boeun 36.29, 127.44 1,274.54 329.78 23 (9) 785.86 (476.89) 
272 Youngju 36.52, 128.31 1,266.32 294.94 18 (9) 634.07 (483.17) 
273 Mungyeong 36.37, 128.09 1,239.07 289.11 26 (9) 843.19 (499.42) 

The values in parentheses represent the maximum duration and severity calculated from the 2013–2015 drought. The bold numbers refer to the case where the maximum duration and severity occurred in 2013–2015.

We also estimated an areal rainfall by using the Thiessen polygon method and computed watershed-wide historic drought durations and severities (Figures 2 and 3). The drought duration and severity variables are positively skewed (Figure 3), with the recent drought state being the most severe and longest duration drought in the post-1974 record.

Figure 2

Drought variables used in this study. Left panel indicates drought duration (month) and right panel represents drought severity (-mm) during the duration.

Figure 2

Drought variables used in this study. Left panel indicates drought duration (month) and right panel represents drought severity (-mm) during the duration.

Figure 3

Scatter plot of duration and severity based on observed precipitation anomalies from 1974 to 2015. The symbol ‘filled star’ indicates current drought state. Gray dots represent randomly generated synthetic pairs using a Frank copula with Lognormal (i.e. duration) and Gamma (i.e. severity) as the two marginal distributions.

Figure 3

Scatter plot of duration and severity based on observed precipitation anomalies from 1974 to 2015. The symbol ‘filled star’ indicates current drought state. Gray dots represent randomly generated synthetic pairs using a Frank copula with Lognormal (i.e. duration) and Gamma (i.e. severity) as the two marginal distributions.

METHODS AND MODEL FORMULATION

Copula function overview

A copula is a multivariate probability function to describe interdependence among multivariate random variables (Sklar 1959). The copula has been widely used for multivariate analysis of hydrologic variables because of its flexibility and power to establish a joint probability function with the different marginals. Favre et al. (2004) showed the efficacy of using copula functions in multivariate hydrologic frequency analysis. Next, we briefly discuss the definition of copulas, and in the remaining part of this section, we present the candidate parametric copula models used in this study and discuss Bayesian parameter estimation.

Suppose we have m-dimensional random variables and their marginal cumulative distribution functions (CDFs) (i.e., ) are continuous. The random variables can be transformed into uniformly distributed marginal distributions by applying the marginal CDFs to the random variables as in Equation (2): 
formula
(2)
The copula of random variables is then defined as a joint CDF (Equation (3)): 
formula
(3)
The marginal CDF describes underlying marginal distribution of variable i, while the copula function C describes the dependence structure between the random variables . A main advantage of the copula approach is that the reverse of the above steps can be applied to efficiently simulate multivariate random samples. Specifically, the required multivariate random variables can be sampled from a uniformly distributed random vector derived from the defined copula function (Equation (4)): 
formula
(4)
where is the quantile function (or inverse CDF) of the marginal distribution, and Equation (3) can be expressed as Equation (5): 
formula
(5)

Marginal distribution

Lognormal (Equation (6)) and Gamma (Equation (7)) distributions were selected as marginal distributions for drought durations (x) and drought severity (y), respectively, based on the Bayesian information criterion (BIC) (see Supplementary information Table S1, available with the online version of this paper). The distribution with the lowest BIC is preferred: 
formula
(6)
 
formula
(7)
where is the mean, is the standard deviation, k is the shape parameter, is the rate parameter, and is the gamma function evaluated at k.

The scatter plot of the synthetic pairs using a Frank copula with Lognormal and Gamma as the two marginal distributions is shown in Figure 3, for visual representation of the dependence structure. A further discussion on the copula selection is provided in the following section. The empirical CDFs along with the theoretical marginal CDFs for four stations (i.e., 105, 108, 201, and 203) are illustrated in Figure 4.

Figure 4

A comparison of CDFs between empirical and theoretical CDFs of the drought duration (i.e. Log-normal distribution) and severity (i.e. Gamma distribution), for four stations (i.e. 105, 108, 201 and 203) as representative.

Figure 4

A comparison of CDFs between empirical and theoretical CDFs of the drought duration (i.e. Log-normal distribution) and severity (i.e. Gamma distribution), for four stations (i.e. 105, 108, 201 and 203) as representative.

Copula types and parameter estimation

This study employed three types of Archimedean copulas to model dependency structure of drought variables, namely, the Clayton, Frank, and Gumbel. All three of the copula functions are defined by one parameter (θ). The Archimedean copula is very flexible and convenient in terms of simulation because a closed functional form is available. The copula parameter, θ, can also be estimated by using general bivariate dependency measures such as Kendall's rank correlation coefficient and Spearman's rank correlation coefficient. The functional forms of the Clayton , Frank , and Gumbel copulas are given in Equations (8)–(12) as follows: 
formula
(8)
 
formula
(9)
 
formula
(10)
and 
formula
(11)
 
formula
(12)
where and indicate marginal CDFs for the drought durations (x) and drought severity (y), respectively.
A Bayesian approach to copula and marginal distribution parameter estimation is proposed to better estimate parameter uncertainty. The posterior distribution of the parameter vector, Θ, ) is given by Bayes' theorem (Equation (13)): 
formula
(13)
where is the parameter vector of the distribution to be fitted, is the likelihood function, D is the observation vector, and is the prior distribution.
A main advantage of the copula approach is that one can estimate joint distributions in two stages. The marginal distributions can first be fit by assuming that the variables are independent identically distributed vectors, then the dependence between the variables can be incorporated through the copula functions. As shown in the previous section, this study assumed that the marginal distributions for drought duration and severity characteristics are Lognormal and Gamma distribution, respectively. Hence, the joint log-likelihood function for the drought characteristics based on the Clayton , Frank , and Gumbel copula functions along with the Lognormal and the Gamma distributions can be written as Equations (14)–(16): 
formula
(14)
 
formula
(15)
 
formula
(16)
 
formula
(16a)
 
formula
(16b)
 
formula
(16c)
where, x and y indicate drought durations and drought severity, respectively.
The prior distributions for the five parameters (e.g., ) are also assigned. Conjugate distributions are probability distributions whose prior and posterior distributions are in the same family. In particular, conjugate distributions are favored for computational reasons. In these contexts, the mean is assumed to be Gaussian with zero-mean and precision 10−3, the standard deviation, scale, and rate parameters are assumed to be gamma, and the copula parameter is assumed to be uniform. Here, the prior distributions for the parameters were given weak informative gamma priors (i.e., gamma (0.1, 0.1), which has mean 1 and variance 10. The Kendal rank correlation coefficients 0.9 is equal to 38.28, 18.00, and 10.00 for Frank, Clayton, and Gumbel copulas, respectively. In this regard, the copula parameter θ is given informative uniform prior distributions to reflect the expected positive dependence structure: 
formula
(17)
 
formula
(18)
 
formula
(19)
 
formula
(20)
 
formula
(21)
Combining the prior parameter distributions (as defined in Equations (17)–(21)) and the log-likelihood function (Equations (14)–(16)) into Equation (13) yields the joint log-posterior distribution of the parameters (Equation (22)): 
formula
(22)
The five parameters defined above are fit using the Gibbs sampler, a Markov Chain Monte Carlo (MCMC) method for estimating the posterior probability density function of model parameters (Gelman & Hill 2006). The Gibbs sampler aims to search the multivariate parameter space by: (1) holding constant all but one of the parameters, (2) simulating from the conditional distribution of the unfixed parameter, and (3) repeating for the rest of the parameters (Kwon et al. 2008). For further information on the use of the Gibbs sampler for Bayesian MCMC, please see Gilks et al. (1995).

Return periods

The univariate return period of a certain drought severity (Shiau & Shen 2001) and drought duration (Shiau 2006) can be estimated using the formula as follows: 
formula
(23a)
 
formula
(23b)
where E(L) is the expected drought interval time, s and d denote drought severity and drought duration, and and represent cumulative distributions of the drought severity and duration, respectively.
Shiau (2006) proposed a methodology that categorizes the return periods of bivariate distributed hydrologic events as joint and conditional return periods. The joint return period for drought duration and severity can be defined in two cases: return period for and and return period for or . The joint return period , the so-called ‘AND’ return period, is that the certain threshold drought duration (d) and severity (s) are exceeded by the respective random variables D and S (Requena et al. 2013). On the other hand, the joint return period , the so-called ‘OR’ return period is that the certain threshold drought duration (d) or severity (s) are exceeded by the respective random variables D and S (Requena et al. 2013): 
formula
 
formula
(24b)
where represents the joint probability between the drought duration and severity. There is an inequality between the univariate return periods and joint return period (Requena et al. 2013) as follows: 
formula
(25)
Detailed discussions on the relationships between the univariate and bivariate return periods can be found in Shiau (2006).

MODEL VALIDATION WITH SYNTHETIC DATA

An experimental study was performed by using synthetic data generated from the three copula functions with a set of predetermined parameters. We used and as the two marginal distributions and used a Spearman's rank correlation coefficient of 0.85 to define the dependence between the marginal distributions. Hence, the true copula parameters corresponding to the Spearman's rank correlation coefficient of 0.85 for Frank, Clayton, and Gumbel copulas are 4.066, 9.563, and 3.015, respectively, as summarized in Table 2. Five hundred correlated random variables were sampled for each of three copula functions (Figure 5).

Table 2

The experimental study for the three copula functions. The 95% credible intervals for the copula and marginal distribution are estimated from the posterior distributions

Copula Parameter True value 2.50% 50% 97.50% 
Clayton copula  3.000 2.784 2.992 3.215 
 1.000 0.912 0.984 1.061 
 1.000 0.970 1.003 1.039 
 0.500 0.516 0.497 0.479 
 4.066 3.474 3.910 4.394 
Frank copula  3.000 2.853 3.097 3.361 
 1.000 0.941 1.030 1.125 
 1.000 0.975 1.006 1.041 
 0.500 0.517 0.496 0.476 
 9.563 8.443 9.371 10.290 
Gumbel copula  3.000 2.613 2.817 3.054 
 1.000 0.867 0.940 1.020 
 1.000 0.977 1.013 1.046 
 0.500 0.532 0.512 0.492 
 3.015 2.777 3.020 3.273 
Copula Parameter True value 2.50% 50% 97.50% 
Clayton copula  3.000 2.784 2.992 3.215 
 1.000 0.912 0.984 1.061 
 1.000 0.970 1.003 1.039 
 0.500 0.516 0.497 0.479 
 4.066 3.474 3.910 4.394 
Frank copula  3.000 2.853 3.097 3.361 
 1.000 0.941 1.030 1.125 
 1.000 0.975 1.006 1.041 
 0.500 0.517 0.496 0.476 
 9.563 8.443 9.371 10.290 
Gumbel copula  3.000 2.613 2.817 3.054 
 1.000 0.867 0.940 1.020 
 1.000 0.977 1.013 1.046 
 0.500 0.532 0.512 0.492 
 3.015 2.777 3.020 3.273 
Figure 5

Synthetic bivariate random variables using (a) Clayton, (b) Frank and (c) Gumbel copula functions.

Figure 5

Synthetic bivariate random variables using (a) Clayton, (b) Frank and (c) Gumbel copula functions.

We generated three MCMC chains of length 10,000 and discarded the first 8,000 as burn-in. The estimated posterior densities of the parameters based on 2,000 MCMC samples from three chains are displayed in Table 2 and Figure 6. We then estimated the Bayesian credible region (defined as the interval that contains 95% of the posterior distribution).

Figure 6

The results of experimental study using synthetic dataset. Each panel shows the posterior density function of the parameters for Clayton, Frank and Gumbel copula model.

Figure 6

The results of experimental study using synthetic dataset. Each panel shows the posterior density function of the parameters for Clayton, Frank and Gumbel copula model.

The true values of the parameters all fall within the credible region in all experimental cases (Table 2). The biases in the posterior medians compared to the true values are relatively small and may be partially explained by sampling error in the generation of the bivariate samples. Table 2 and Figure 6 suggest that the proposed Bayesian approach is suitable for copula-based multivariate analysis.

Another central question is whether the posterior likelihood can be used to identify the correct generating copula model. To answer this, we again generated three bivariate random samples from the three known copulas. For each of these samples, we fit three models (Clayton, Frank, and Gumbel copulas) and compare the deviance information criterion (DIC) from the posterior likelihood functions of each of the models. The DIC is especially useful for hierarchical Bayesian model selection when the posterior distributions have been derived by MCMC simulation. The DIC is defined as the sum of the posterior mean of the deviance and the effective number of parameters (pD) estimated from the posterior distributions of the model parameters.

The estimated posterior distributions of the parameters along with the DIC value are summarized in Table 3. Lower DIC values indicate better model fit to the data. The results indicate that the proposed modeling scheme is capability of correctly identifying the data generating copula function.

Table 3

The experimental study for the selection of correct copula functions based on the DIC value. The presented values are median values estimated from the marginal posterior distributions of the parameters

Model Lognormal
 
Gamma
 
Copula Log-likelihood DIC 
Parameter  (3.0)  (1.0)  (1.0)  (0.5)  True/Sim. 
Case A (Clayton copula) 
 Clayton 2.725 1.155 0.989 0.741 4.066/3.900 36.996 329.450 
 Frank 2.723 1.114 0.984 0.732 9.563/10.663 31.810 340.180 
 Gumbel 2.721 1.161 0.999 0.752 3.015/2.353 20.739 361.840 
Case B (Frank copula) 
 Clayton 2.757 1.151 0.999 0.747 4.066/2.633 27.253 349.37 
 Frank 2.984 1.048 0.978 0.734 9.563/12.877 39.341 324.68 
 Gumbel 2.817 1.122 0.978 0.758 3.015/3.007 31.462 339.76 
Case C (Gumbel copula) 
 Clayton 3.472 0.795 0.878 0.662 4.066/2.445 24.585 310.72 
 Frank 3.259 0.865 0.891 0.646 9.563/9.680 30.173 301.19 
 Gumbel 3.457 0.794 0.868 0.612 3.015/3,004 35.907 284.17 
Model Lognormal
 
Gamma
 
Copula Log-likelihood DIC 
Parameter  (3.0)  (1.0)  (1.0)  (0.5)  True/Sim. 
Case A (Clayton copula) 
 Clayton 2.725 1.155 0.989 0.741 4.066/3.900 36.996 329.450 
 Frank 2.723 1.114 0.984 0.732 9.563/10.663 31.810 340.180 
 Gumbel 2.721 1.161 0.999 0.752 3.015/2.353 20.739 361.840 
Case B (Frank copula) 
 Clayton 2.757 1.151 0.999 0.747 4.066/2.633 27.253 349.37 
 Frank 2.984 1.048 0.978 0.734 9.563/12.877 39.341 324.68 
 Gumbel 2.817 1.122 0.978 0.758 3.015/3.007 31.462 339.76 
Case C (Gumbel copula) 
 Clayton 3.472 0.795 0.878 0.662 4.066/2.445 24.585 310.72 
 Frank 3.259 0.865 0.891 0.646 9.563/9.680 30.173 301.19 
 Gumbel 3.457 0.794 0.868 0.612 3.015/3,004 35.907 284.17 

The Log-likelihoods and DIC values in bold are the optimum copulas selected in the experimental study.

The observed drought data used in this work contain ‘Ties’ (i.e., the repeated measurements, especially concerning the duration), which may adversely affect the identification of the generating copula during the experimental study performed in this study. The identification problem when Ties are present has been investigated and the results were negative (De Michele et al. 2013). In this regard, we performed a number of simulations to test the sensitivity of the Ties' effects. As performed in the previous experimental study, three bivariate random samples from the three known copulas were generated, and the generated drought duration is then rounded to the nearest integer value to consider the Ties' effects, as displayed in Figure S2 (available with the online version of this paper). For each of these samples, we fit three models (Clayton, Frank, and Gumbel copulas) and compare the DIC from the posterior likelihood functions of each of the models. The estimated posterior distributions of the parameters along with the DIC value are summarized in Table S2 (available online). The results indicate that the proposed modeling scheme is still correctly identifying the data generating copula function. However, for the Clayton copula function, the DIC values are relatively similar to that of the other copulas, meaning that the Clayton copula function is likely sensitive to the Ties.

BIVARIATE DROUGHT FREQUENCY ANALYSIS

Model identification and parameter estimation

As with the synthetic dataset above, we use the Clayton, Frank, and Gumbel copulas as candidates for modeling dependency structure of drought variables in the Han River watershed. The marginal distributions for drought duration and severity were assumed to be the Lognormal and Gamma, respectively. The DIC value was calculated for the three copulas under the same marginal distributions and used to select the best model. This was done for a network of 18 weather stations (Table 4). The Frank copula is selected for 13 stations, and the Gumbel for the remaining stations. The lower DIC values for the Gumbel copula compared to the Frank copula indicate that upper tail dependence is an important feature of the drought duration and severity data. The selected copula functions are used for subsequent analysis.

Table 4

DIC values of the three different copulas for the weather stations

Station no. Station name DIC value
 
Clayton Frank Gumbel 
90 Sokcho 698.897 670.458 680.434 
100 Daeqwalyeong 664.055 650.905 650.357 
101 Chuncheon 845.887 819.694 828.498 
105 Gangneung 759.797 716.386 708.380 
108 Seoul 860.898 844.142 846.563 
112 Incheon 837.784 825.275 828.659 
114 Wonju 837.889 820.206 830.867 
119 Suwon 799.869 793.611 792.983 
131 Chongju 811.919 785.285 782.707 
201 Ganghwa 890.634 862.785 881.893 
202 Yangpyeong 859.953 844.063 845.178 
203 Icheon 844.625 832.102 842.248 
211 Inje 803.735 789.577 789.920 
212 Hongcheon 813.667 794.247 801.213 
221 Jecheon 825.814 821.015 813.988 
226 Boeun 745.422 721.203 724.536 
272 Youngju 820.933 801.918 805.088 
273 Mungyeong 714.144 696.728 698.273 
Station no. Station name DIC value
 
Clayton Frank Gumbel 
90 Sokcho 698.897 670.458 680.434 
100 Daeqwalyeong 664.055 650.905 650.357 
101 Chuncheon 845.887 819.694 828.498 
105 Gangneung 759.797 716.386 708.380 
108 Seoul 860.898 844.142 846.563 
112 Incheon 837.784 825.275 828.659 
114 Wonju 837.889 820.206 830.867 
119 Suwon 799.869 793.611 792.983 
131 Chongju 811.919 785.285 782.707 
201 Ganghwa 890.634 862.785 881.893 
202 Yangpyeong 859.953 844.063 845.178 
203 Icheon 844.625 832.102 842.248 
211 Inje 803.735 789.577 789.920 
212 Hongcheon 813.667 794.247 801.213 
221 Jecheon 825.814 821.015 813.988 
226 Boeun 745.422 721.203 724.536 
272 Youngju 820.933 801.918 805.088 
273 Mungyeong 714.144 696.728 698.273 

The DIC values in bold are the optimum copulas selected in this study.

The marginal posterior distributions for the parameters for each station are shown in Figure 7. The credible intervals of the marginal posterior distributions are obtained from three chains based on 10,000 iterations and with the first 8,000 discarded as burn-in. The parameters of the marginal distributions associated with drought duration (i.e., μ and ) appear to be more consistent across stations compared to the parameters of the drought severity (i.e., κ and τ).

Figure 7

Boxplots of the marginal posterior distributions for the parameters obtained from three chains based on 10,000 iterations. The symbol ‘-’ indicates a median value of the posteriors. The stations, labeled as 100, 105, 119, 131 and 221, are based on the Gumbel copula while Frank copula for the remains of stations.

Figure 7

Boxplots of the marginal posterior distributions for the parameters obtained from three chains based on 10,000 iterations. The symbol ‘-’ indicates a median value of the posteriors. The stations, labeled as 100, 105, 119, 131 and 221, are based on the Gumbel copula while Frank copula for the remains of stations.

Return period estimation for the 2013–2015 drought

The univariate drought duration and drought severity associated with the selected return levels (2, 5, 10, 30, 50, and 100-year) are shown in Table S3 (available with the online version of this paper). The expected recurrence interval between successive drought events, E(L), is from about 8 to 12 months across the weather stations. A recommended design criterion for water storage structures including reservoirs is the 30-year return period with either rainfall amounts or discharge in South Korea. The drought durations and severities corresponding to a 30-year return period for weather stations in the Han River watershed have a median (range) of about 23 (20–30) months and 508 (404–965) mm. Considering that many water storage structures are designed for the 30-year drought, a prolonged rainfall deficit of over 500 mm over 20 months could lead to a higher likelihood of sustained water shortage throughout the watershed. Given this, it is startling to note that the rainfall deficits of the 2013–2015 drought range from about 400 to 1,400 mm. The estimated univariate return period for the 2013–2015 drought's duration and severity among the stations analyzed here ranges from about 4 to 31 years and 17 to 6,105 years, respectively. The univariate and joint return periods of the recent 2013–2015 drought are displayed in Table 5. Some of the 2013–2015 drought durations (i.e., 9–24 months) have median return periods of around 30 years while the drought severities at certain stations may be regarded as extreme events with return periods in the hundreds of years. It is clear that the recent drought is extraordinary at many stations in the Han River watershed although the associated uncertainty is quite large (the 95% credible interval at the Seoul station, for example, has a range of about 180 to 2,250 years). Given these results and the inequality described in Equation (25), we see that the joint return periods are much higher than those of the univariate.

Table 5

The univariate and joint median return periods and their 95% credible intervals for the 2013–2015 drought event

Station no. Return periods
 
    
90 12.35 (8.58, 18.73) 62.10 (36.40, 142.14) 85.07 (36.90, 228.67) 11.70 (8.31, 17.48) 
100 12.74 (8.69, 19.69) 25.92 (15.59, 47.27) 26.17 (15.76, 47.69) 12.69 (8.65, 19.57) 
101 25.22 (15.24, 43.64) 1,484.67 (392.77, 7,366.58) 4,601.68 (859.00, 30,572.91) 24.90 (15.07, 43.16) 
105 14.25 (9.59, 21.70) 17.19 (11.11, 27.94) 18.92 (12.24, 30.49) 13.20 (8.97, 19.97) 
108 31.07 (18.27, 57.05) 554.38 (180.16, 2,253.62) 2,048.99 (472.31, 11,989.83) 29.72 (17.62, 54.57) 
112 24.65 (15.47, 41.56) 398.42 (131.63, 1,252.64) 1,125.54 (275.09, 4,812.24) 23.60 (14.83, 39.83) 
114 4.86 (3.85, 6.27) 19.38 (12.17, 33.50) 23.31 (13.84, 42.44) 4.67 (3.75, 5.90) 
119 30.03 (17.74, 54.81) 83.38 (39.65, 196.47) 84.35 (40.37, 197.69) 29.88 (17.67, 54.44) 
131 4.50 (3.61, 5.80) 22.43 (13.89, 38.48) 22.53 (13.98, 38.67) 4.49 (3.60, 5.79) 
201 29.68 (17.16, 55.15) 6,104.60 (1,165.77, 43,064.00) 22,264.77 (2,977.98, 227,804.71) 29.54 (17.10, 54.79) 
202 23.75 (14.71, 41.42) 393.50 (133.86, 1,524.36) 1,216.76 (309.56, 6,712.96) 22.65 (14.11, 39.29) 
203 4.88 (3.88, 6.38) 33.37 (19.29, 63.64) 39.88 (22.21, 80.03) 4.77 (3.82, 6.16) 
211 22.26 (13.82, 39.16) 176.96 (68.42, 552.63) 459.47 (129.50, 1,989.27) 20.62 (12.97, 35.62) 
212 22.01 (13.82, 38.03) 1,059.75 (272.35, 4,667.83) 2,765.07 (538.07, 16,453.46) 27.71 (13.63, 37.43) 
221 5.25 (4.04, 6.95) 17.75 (11.21, 29.17) 18.18 (11.46, 29.92) 5.21 (4.01, 6.88) 
226 4.22 (3.46, 5.29) 17.32 (10.84, 30.17) 18.33 (11.30, 32.56) 4.17 (3.44, 5.19) 
272 4.56 (3.68, 5.85) 24.12 (14.22, 44.90) 27.28 (15.47, 52.86) 4.47 (3.62, 5.66) 
273 4.38 (3.56, 5.56) 20.24 (12.26, 36.62) 21.94 (12.95, 41.23) 4.30 (3.53, 5.41) 
Station no. Return periods
 
    
90 12.35 (8.58, 18.73) 62.10 (36.40, 142.14) 85.07 (36.90, 228.67) 11.70 (8.31, 17.48) 
100 12.74 (8.69, 19.69) 25.92 (15.59, 47.27) 26.17 (15.76, 47.69) 12.69 (8.65, 19.57) 
101 25.22 (15.24, 43.64) 1,484.67 (392.77, 7,366.58) 4,601.68 (859.00, 30,572.91) 24.90 (15.07, 43.16) 
105 14.25 (9.59, 21.70) 17.19 (11.11, 27.94) 18.92 (12.24, 30.49) 13.20 (8.97, 19.97) 
108 31.07 (18.27, 57.05) 554.38 (180.16, 2,253.62) 2,048.99 (472.31, 11,989.83) 29.72 (17.62, 54.57) 
112 24.65 (15.47, 41.56) 398.42 (131.63, 1,252.64) 1,125.54 (275.09, 4,812.24) 23.60 (14.83, 39.83) 
114 4.86 (3.85, 6.27) 19.38 (12.17, 33.50) 23.31 (13.84, 42.44) 4.67 (3.75, 5.90) 
119 30.03 (17.74, 54.81) 83.38 (39.65, 196.47) 84.35 (40.37, 197.69) 29.88 (17.67, 54.44) 
131 4.50 (3.61, 5.80) 22.43 (13.89, 38.48) 22.53 (13.98, 38.67) 4.49 (3.60, 5.79) 
201 29.68 (17.16, 55.15) 6,104.60 (1,165.77, 43,064.00) 22,264.77 (2,977.98, 227,804.71) 29.54 (17.10, 54.79) 
202 23.75 (14.71, 41.42) 393.50 (133.86, 1,524.36) 1,216.76 (309.56, 6,712.96) 22.65 (14.11, 39.29) 
203 4.88 (3.88, 6.38) 33.37 (19.29, 63.64) 39.88 (22.21, 80.03) 4.77 (3.82, 6.16) 
211 22.26 (13.82, 39.16) 176.96 (68.42, 552.63) 459.47 (129.50, 1,989.27) 20.62 (12.97, 35.62) 
212 22.01 (13.82, 38.03) 1,059.75 (272.35, 4,667.83) 2,765.07 (538.07, 16,453.46) 27.71 (13.63, 37.43) 
221 5.25 (4.04, 6.95) 17.75 (11.21, 29.17) 18.18 (11.46, 29.92) 5.21 (4.01, 6.88) 
226 4.22 (3.46, 5.29) 17.32 (10.84, 30.17) 18.33 (11.30, 32.56) 4.17 (3.44, 5.19) 
272 4.56 (3.68, 5.85) 24.12 (14.22, 44.90) 27.28 (15.47, 52.86) 4.47 (3.62, 5.66) 
273 4.38 (3.56, 5.56) 20.24 (12.26, 36.62) 21.94 (12.95, 41.23) 4.30 (3.53, 5.41) 

The values in parentheses show the 95% credible intervals.

The joint distribution of drought duration and severity based on the copula functions and the Lognormal and Gamma marginal distributions was used to estimate the joint return periods for (duration and severity, D ≥ d and S ≥ s) and (duration or severity, D ≥ d or S ≥ s) using Equations 24(a) and 24(b) (Table 5). There was a range of 40 to 59 drought events (as per our precipitation anomaly-based definition in the section ‘Han River watershed and drought definitions') over the 42 years across the weather stations. The joint return period for the 2013–2015 drought is, in general, estimated to be greater than 100 years; estimated return periods are over 1,000 years at some stations. It is clear that the unprecedented joint return periods of the 2013–2015 drought at many stations are consistent with the very high drought severity of the current event. Contour plots of the joint return periods along with the 95% Bayesian credible intervals for selected stations are shown in Figure 8. The figure shows that the uncertainty ranges increase as the return periods increase. The contours also illustrate the variety of combinations of drought durations and severities associated with a given joint return period. The most severe state for the 2013–2015 drought (as measured by the joint return period) is associated with a drought severity of about 1,400 mm and a drought duration of 24 months at station 201 (Ganghwa).

Figure 8

The joint return periods of droughts based on Bayesian Copula model for the four stations (i.e. 105, 108, 201 and 203). The symbol ‘filled star’ indicates the 2013–2015 drought and the shaded area indicates the uncertainty associated with the return periods.

Figure 8

The joint return periods of droughts based on Bayesian Copula model for the four stations (i.e. 105, 108, 201 and 203). The symbol ‘filled star’ indicates the 2013–2015 drought and the shaded area indicates the uncertainty associated with the return periods.

A spatial distribution of the estimated median joint return period of the 2013–2015 drought for the Han River watershed (along with the Bayesian 95% credible interval (Figure S3, available online)) is illustrated in Figure 9. Here, a biharmonic spline interpolation technique is used to construct the map of return period. The western part of the Han River watershed has, in general, experienced the most severe droughts in recorded history over the last four decades.

Figure 9

Spatial distribution of the posterior median of the joint return periods. The symbol ‘+’ indicates the stations used in this study. The thick black line encloses sub-catchment region.

Figure 9

Spatial distribution of the posterior median of the joint return periods. The symbol ‘+’ indicates the stations used in this study. The thick black line encloses sub-catchment region.

Finally, the conditional return periods of the drought duration and severity, given drought severity and duration exceeding a certain threshold, were estimated for the four stations, respectively, as shown in Figure 10. Please refer to Shiau (2006) for more detailed formulation of the conditional return periods. Figure 10(a)10(d) demonstrate the conditional return periods of the drought duration given the drought severity while Figure 10(e)10(h) represent the conditional return periods of the drought severity given the drought duration. For example, the estimated conditional return period of drought severity given a duration of 24 months for the station 108 (Seoul) is over 100 years.

Figure 10

The left panels show the conditional return periods of drought duration given the drought severity while the right panels are the condition return periods of drought severity given the drought duration for the stations 105, 108, 201 and 203. The symbol ‘filled star’ represents the 2013–2015 drought state.

Figure 10

The left panels show the conditional return periods of drought duration given the drought severity while the right panels are the condition return periods of drought severity given the drought duration for the stations 105, 108, 201 and 203. The symbol ‘filled star’ represents the 2013–2015 drought state.

SUMMARY AND CONCLUSIONS

A copula-based multivariate frequency analysis was employed to model the dependences in drought duration and severity. Copula-based multivariate models have been widely applied to hydrologic frequency analyses, however, reliable estimation of the parameters is often a challenging task in the multivariate setting. Moreover, uncertainties associated with the model parameters in copula-based frequency models are not often properly addressed. We have used a Bayesian approach for parameter estimation in a copula-based bivariate drought frequency analysis. The proposed model was applied to estimate joint return periods on a network of 18 rainfall stations in the Han River watershed in South Korea. The key results and conclusions are summarized as follows:

  1. A validation of the proposed Bayesian parameter estimation approach was conducted based on synthetic data that were generated by a known set of parameters along with the different copulas (Clayton, Frank, and Gumbel). The validation exercise showed that the true underlying parameters could be reasonably estimated.

  2. For the univariate return periods, some of the 2013–2015 drought durations (i.e., 9–24 months) had return periods of around 30 years while the drought severities (which range from about 400 to 1,400 mm) at certain stations may be regarded as extreme events with return periods in the hundreds of years. It was clear that the recent drought was extraordinary at many stations in the Han River watershed although the associated uncertainty is quite large (e.g., for the Seoul station, a range of about 180 to 2,250 years).

  3. The main point of the paper was to show that the 2013–2015 drought is very rare. We confirmed this finding by doing the uncertainty analysis in the Bayesian framework. It was clear that the unprecedented univariate and joint return periods of the 2013–2015 drought at many stations were consistent with the very high drought severity of the current event. The estimated joint return period of the drought duration and severity were greater than 100 years for many stations and upwards of 1,000 years for a few stations. The most severe state for the 2013–2015 drought was associated with a drought severity of about 1,400 mm and a drought duration of 24 months at the station 201 (Ganghwa).

However, it should be noted that the estimated return periods for some stations are rather unrealistic and their uncertainties are quite large. This is mainly associated with a relatively small sample size, which introduces sampling uncertainty and variation about the parameters, in addition to the usual level of uncertainty. We have not intended to explore the potential reduction of uncertainties in estimating the return periods at this moment and future work will focus on the use of a hierarchical Bayesian model within a regional drought frequency analysis.

ACKNOWLEDGEMENTS

We thank Upmanu Lall and David Farnham from the Water Center, Columbia University, for providing feedback on the manuscript. We would like to thank the two anonymous reviewers for providing us with suggestions and insightful comments that improved the manuscript. This research was supported by a Grant (15AWMP-B079625-02) from Water Management Research Program funded by Ministry of Land, Infrastructure and Transport of the Korean Government. The data used in this study are available upon request from the corresponding author via email (hkwon@jbnu.ac.kr). The authors declare that there is no conflict of interest regarding the publication of this article.

REFERENCES

REFERENCES
Arena
C.
,
Cannarozzo
M.
&
Mazzola
M. R.
2006
Multi-year drought frequency analysis at multiple sites by operational hydrology – A comparison of methods
.
Physics and Chemistry of the Earth, Parts A/B/C
31
(
18
),
1146
1163
.
Bracken
C.
,
Rajagopalan
B.
,
Cheng
L.
,
Kleiber
W.
&
Gangopadhyay
S.
2016
Spatial Bayesian hierarchical modeling of precipitation extremes over a large domain
.
Water Resources Research
52
(
8
),
6643
6655
.
Cebrián
A. C.
&
Abaurrea
J.
2006
Drought analysis based on a marked cluster Poisson model
.
Journal of Hydrometeorology
7
(
4
),
713
723
.
De Michele
C.
,
Salvadori
G.
,
Vezzoli
R.
&
Pecora
S.
2013
Multivariate assessment of droughts: frequency analysis and dynamic return period
.
Water Resources Research
49
(
10
),
6985
6994
.
Favre
A. C.
,
El Adlouni
S.
,
Perreault
L.
,
Thiémonge
N.
&
Bobée
B.
2004
Multivariate hydrological frequency analysis using copulas
.
Water Resources Research
40
(
1
),
W01101
.
Gelman
A.
&
Hill
J.
2006
Data Analysis Using Regression and Multilevel/Hierarchical Models
.
Cambridge University Press
,
New York
.
Genest
C.
,
Favre
A. C.
,
Béliveau
J.
&
Jacques
C.
2007
Metaelliptical copulas and their use in frequency analysis of multivariate hydrological data
.
Water Resources Research
43
(
9
),
W09401
.
Gilks
W. R.
,
Best
N.
&
Tan
K.
1995
Adaptive rejection Metropolis sampling within Gibbs sampling
.
Applied Statistics
44
(
4
),
455
472
.
Gräler
B.
,
van den Berg
M.
,
Vandenberghe
S.
,
Petroselli
A.
,
Grimaldi
S.
,
De Baets
B.
&
Verhoest
N.
2013
Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation
.
Hydrology and Earth System Sciences
17
(
4
),
1281
1296
.
Hao
Z.
&
Singh
V.
2012
Entropy-copula method for single-site monthly streamflow simulation
.
Water Resources Research
48
(
6
),
W06604
,
doi:10.1029/2011WR011419
.
Hao
Z.
&
Singh
V.
2013
Modeling multisite streamflow dependence with maximum entropy copula
.
Water Resources Research
49
(
10
),
7139
7143
.
Klein
B.
,
Pahlow
M.
,
Hundecha
Y.
&
Schumann
A.
2009
Probability analysis of hydrological loads for the design of flood control systems using copulas
.
Journal of Hydrologic Engineering
15
(
5
),
360
369
.
Kwon
H. H.
,
Sivakumar
B.
,
Moon
Y. I.
&
Kim
B. S.
2011
Assessment of change in design flood frequency under climate change using a multivariate downscaling model and a precipitation-runoff model
.
Stochastic and Environmental Research and Risk Assessment
25
(
4
),
567
581
.
Liu
D.
,
Wang
D.
,
Wang
L.
,
Chen
Y.
,
Chen
X.
&
Gu
S.
2015
POME-copula for hydrological dependence analysis
.
Proceedings of the International Association of Hydrological Sciences
368
,
251
256
.
Loganathan
G.
,
Kuo
C.
&
McCormick
T.
1985
Frequency analysis of low flows
.
Hydrology Research
16
(
2
),
105
128
.
Ma
M.
,
Song
S.
,
Ren
L.
,
Jiang
S.
&
Song
J.
2013
Multivariate drought characteristics using trivariate Gaussian and student t copulas
.
Hydrological Processes
27
(
8
),
1175
1190
.
Madadgar
S.
&
Moradkhani
H.
2011
Drought analysis under climate change using copula
.
Journal of Hydrologic Engineering
18
(
7
),
746
759
.
Madadgar
S.
&
Moradkhani
H.
2014
Improved Bayesian multimodeling: integration of copulas and Bayesian model averaging
.
Water Resources Research
50
(
12
),
9586
9603
.
Parent
E.
,
Favre
A.-C.
,
Bernier
J.
&
Perreault
L.
2014
Copula models for frequency analysis what can be learned from a Bayesian perspective
?
Advances in Water Resources
63
,
91
103
.
Raynal-Villasenor
J.
,
Salas
J.
1987
Multivariate extreme value distributions in hydrological analyses
. In:
Water for the Future: Hydrology in Perspective
(
Rodda
J .C.
&
Matalas
N. C.
, eds).
IAHS Publication
,
Wallingford
,
UK
(164)
.
Requena
A.
,
Mediero Orduña
L.
&
Garrote de Marcos
L.
2013
A bivariate return period based on copulas for hydrologic dam design: accounting for reservoir routing in risk estimation
.
Hydrology and Earth System Sciences
17
(
8
),
3023
3038
.
Sadri
S.
&
Burn
D.
2012
Nonparametric methods for drought severity estimation at ungauged sites
.
Water Resources Research
48
(
12
),
W12505
.
Salas
J. D.
,
Fu
C. J.
,
Cancelliere
A.
,
Dustin
D.
,
Bode
D.
,
Pineda
A.
&
Vincent
E.
2005
Characterizing the severity and risk of drought in the Poudre River, Colorado
.
Journal of Water Resources Planning and Management
131
(
5
),
383
393
.
Sarhadi
A.
,
Burn
D. H.
,
Concepción Ausín
M.
&
Wiper
M. P.
2016
Time varying nonstationary multivariate risk analysis using a dynamic Bayesian copula
.
Water Resources Research
52
(
3
),
2327
2349
.
Shiau
J. T.
2006
Fitting drought duration and severity with two-dimensional copulas
.
Water Resources Management
20
(
5
),
795
815
.
Shiau
J. T.
&
Shen
H. W.
2001
Recurrence analysis of hydrologic droughts of differing severity
.
Journal of Water Resources Planning and Management
127
(
1
),
30
40
.
Sklar
M.
1959
Fonctions de répartition à n dimensions et leurs marges
.
Université Paris 8
,
Paris
,
France
.
Yoo
J.
,
Kwon
H. H.
,
Kim
T. W.
&
Ahn
J. H.
2012
Drought frequency analysis using cluster analysis and bivariate probability distribution
.
Journal of Hydrology
420
,
102
111
.
Yue
S.
,
Ouarda
T.
&
Bobee
B.
2001
A review of bivariate gamma distributions for hydrological application
.
Journal of Hydrology
246
(
1
),
1
18
.
Zhang
Q.
,
Xiao
M.
,
Singh
V. P.
&
Chen
X.
2013
Copula-based risk evaluation of hydrological droughts in the East River basin, China
.
Stochastic Environmental Research and Risk Assessment
27
(
6
),
1397
1406
.

Supplementary data