Estimation of return period and its uncertainty for the recent 2013–2015 drought in the Han River watershed in South Korea

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.

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 km², and the annual average discharge of the basin is about 613 m³/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.

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.

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.

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.

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.

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).

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).

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 (p_D) 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.

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.

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.

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 τ).

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.

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).

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.

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.

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.

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 ([email protected]). The authors declare that there is no conflict of interest regarding the publication of this article.