Yun-bio Wu , Lin-qing Xue ,c,*, Yun-hong Liu

a College of Hy drology and Water Resources, Hohai University, Nanjing 210098, China

b Department of Fundamental Education, Wanjiang University of Technology, Maanshan 243031, China

c College of Water and Architectural Engineering, Shihezi University, Shihezi 832003, China

Abstract This study developed a hierarchical Bayesian (HB)model for local and regional flood frequency analysis in the Dongting Lake Basin, in China.The annual maximum daily flows from 15 streamflow-gauged sites in the study area were analyzed with the HB model.The generalized extreme value (GEV)distribution was selected as the extreme flood distribution, and the GEV distribution location and scale parameters were spatially modeled through a regression approach with the drainage area as a covariate.The Markov chain Monte Carlo (MCMC)method with Gibbs sampling was employed to calculate the posterior distribution in the HB model.The results showed that the proposed HB model provided satisfactory Bayesian credible intervals for flood quantiles,while the traditional delta method could not provide reliable uncertainty estimations for large flood quantiles,due to the fact that the lower confidence bounds tended to decrease as the return periods increased.Furthermore,the HB model for regional analysis allowed for a reduction in the value of some restrictive assumptions in the traditional index flood method,such as the homogeneity region assumption and the scale invariance assumption.The HB model can also provide an uncertainty band of flood quantile prediction at a poorly gauged or ungauged site, but the index flood method with L-moments does not demonstrate this uncertainty directly.Therefore,the HB model is an effective method of implementing the flexible local and regional frequency analysis scheme,and of quantifying the associated predictive uncertainty.

Keywords:Flood frequency analysis; Hierarchical Bayesian model; Index flood method; Generalized extreme value distribution; Dongting Lake Basin

1.Introduction

Flood risk assessment and design of hydrologic structures(e.g., dams, bridges, and spillways)require robust estimation of return levels and associated uncertainties of extreme flooding (Lima and Lall, 2010; Bracken et al., 2016; Halbert et al., 2016; Li et al., 2019; Dong and Xie, 2004).Statistical modeling is an important approach (Katz et al., 2002; Huo et al., 2018).In such cases, a specified probability distribution is selected to fit the observed flood data, and then the quantile of the fitted probability distribution associated with a specified non-exceedance probability (e.g., 1% for the 100-year flood quantile)is estimated.However, the observed data series are rarely long enough to provide reliable flood quantile estimates to meet the requirements of engineering design(Halbert et al., 2016; Wu et al., 2019), leading to a large amount of uncertainties in the estimates.

Many efforts have been devoted to reducing the uncertainties of estimates in flood frequency analyses.The most frequently used method is integrating additional information into flood frequency analyses, including two types of information extension (Gaume et al., 2010; Halbert et al., 2016).One is temporal extension of the data set to incorporate historical or paleoflood data (Hosking and Wallis, 1986a, 1986b;O'Connell et al.,2002;Reis and Stedinger,2005;Gaume et al.,2010; Nguyen et al., 2014).The other is spatial extension of the data set, which aggregates statistically homogeneous data to build a large regional data sample (Hosking and Wallis,1997; Renard, 2011; Najafi and Moradkhani, 2013).

Regional frequency analysis (RFA)is a popular method in hydrologic frequency analysis.It improves the accuracy of estimates by borrowing information from the neighboring sites,and has been shown to be superior to the standard at-site frequency analysis (Dalrymple, 1960; Hosking and Wallis,1997; Renard, 2011; Najafi and Moradkhani, 2013).The index flood method proposed by Dalrymple (1960)is a classical method of implementing RFA schemes, which is still widely used in engineering practice.Hosking and Wallis(1997)further improved this method by using L-moments.However, because the implementation of the index flood method is based on some restrictive assumptions, its application may be limited (Renard, 2011; Lima et al., 2016).For example, the scale invariance assumption is too restrictive in some cases and the delineation of homogeneous regions is rare in practice.In addition, the quantification of the total predictive uncertainty is another challenging task in the index flood method(Renard,2011;Najafi and Moradkhani,2013;Yan and Moradkhani, 2014; Lima et al., 2016).

Recently,an alternative approach based on the hierarchical Bayesian(HB)model has been proposed(Cooley et al.,2007;Renard, 2011; Najafi and Moradkhani, 2013, 2014; Gelman et al., 2014; Steinschneider and Lall, 2015; Bracken et al.,2016; Lima et al., 2016).The HB model has emerged as an RFA approach that is superior in many aspects to the traditional index flood method with L-moments.For example, this approach relaxes the assumption of scale invariance and quantifies uncertainties in parameter and flood quantile estimates directly (Najafi and Moradkhani, 2013, 2014; Lima et al., 2016).Cooley et al.(2007)first introduced the HB spatial model in extreme precipitation analysis, and Bracken et al.(2016)extended this work to model spatial extremes in a large domain using a Gaussian elliptical copula.Renard(2011)proposed a general HB framework to implement RFA schemes that overcome some difficulties in the traditional RFA method.Najafi and Moradkhani(2013)developed a spatial HB method to model the extreme runoffs by incorporating the latitude, longitude, elevation, and drainage area into the generalized Pareto distribution (GPD)scale parameter estimation.Yan and Moradkhani (2014)proposed a regional HB model for flood frequency analysis by adding a newL--moments layerin the model to avoid the subjective selection of flood distribution.Lima et al.(2016)proposed a multilevel HB model for improving local and regional flood frequency analysis.

Based on the advantages of the HB model,in this study,we used an HB model to analyze local(at-site)and regional flood frequency, and to quantify the uncertainties in the parameter estimation and flood prediction at both gauged and ungauged sites.The proposed HB model was built by modeling the parameters of flood probability distribution using the regression approach with spatial information for gauged sites, in which the specific catchment characteristics (e.g., drainage area, elevation, and slope)were considered as covariates.The desired flood statistics at ungauged sites in the HB model were predicted by transferring information about the gauged sites to the ungauged sites through the regression approach.

2.Study area and data

2.1.Study area

The Dongting Lake Basin is located in the middle and lower reaches of the Yangtze River(approximately 24°36′N(xiāo) to 30°27′N(xiāo), 107°26′E to 114°20′E), with a total area of 2.63 × 105km2(Fig.1).Dongting Lake, the second largest freshwater lake in China,with an area about 2625 km2,lies in the northeastern part of the basin.The basin is surrounded by mountains on the east, south, and west sides and opens to the north, forming a unique horseshoe-shape topography.Dongting Lake is fed by four major tributaries (the Xiangjiang,Zishui, Yuanjiang, and Lishui rivers)and three outfalls (the Songzi, Hudu, and Ouchi outfalls)of the Yangtze River.Eventually, the water drains back into the Yangtze River directly (Yuan et al., 2016).It plays a very important role in the diversion and storage of floods, and the comprehensive utilization of water resources in the middle reaches of the Yangtze River.

Fig.1.Spatial distribution of selected 15 streamflow-gauged sites in Dongting Lake Basin.

The Dongting Lake Basin lies in the subtropical monsoon climate zone, with a wet season between July and September and a dry season between November and the following February (Yuan et al., 2016).This area belongs to the middle part of three zones of China that frequently face disasters,and has become one of the most severe disaster-prone regions,especially for floods and droughts(Xiong et al.,2009).Floods have become one of the greatest obstacles to the sustainable development of agriculture in the Dongting Lake Basin.However, few studies have investigated the local and regional flood frequency and quantified the uncertainties in the flood prediction in this area.

2.2.Data

Annual maximum daily flows from 15 gauged sites within the study area were collected from the Hydrology and Water Resources Survey Bureau of Hunan Province,China.All these data are quality-controlled before their release.Data availability for all the gauged sites is shown in Fig.2 and the details of each site are shown in Table 1.The streamflow data are available over the period of 1951-2014,but not all sites have complete records.Most sites have data lengths of less than 30 years, with only sites 1, 7, and 15 having complete records.Site 6 has only five years of records.The incompleteness of data makes the traditional model parameter estimation challenging and highlights the necessity of using other new methods.

3.Methodology

3.1.Generalized extreme value (GEV)distribution

In extreme value theory, the maximum data in annual,seasonal, or monthly blocks are assumed to follow the generalized extreme value (GEV)distribution based on the asymptotic theory (Jenkinson, 1955; Coles, 2001; Najafi and Moradkhani, 2013).In this study, the GEV distribution was used to model flood extremes,because of its ability to capture a wide range of tail behaviors and its wide use in previous hydrologic studies (Katz et al., 2002; Morrison and Smith,2002; Coles et al., 2003; Renard et al., 2006; Ouarda and El-Adlouni, 2011; Nguyen et al., 2014; Assis et al., 2017;Bracken et al., 2018).

Table 1 Location,drainage area,and data length of 15 gauged sites used in this study.

The GEV cumulative distribution function is as follows:

where μ, σ, and ξ are the location, scale, and shape parameters, respectively, with μ∈R, σ ＞0, ξ∈R, and 1+ ξ（x-μ）/σ ＞0.The GEV distribution is equivalent to a Weibull distribution if ξ ＜0, to a Fr′echet distribution if ξ ＞0, and to a Gumbel distribution if ξ = 0.

In flood frequency analysis, hydrologists are usually interested in calculation of the designed flood for aT-year return period,i.e.,the predicted value of the quantilexp(p=1/T)of the GEV distribution.The quantile at each site can be obtained by inversing Eq.(1):

3.2.Hierarchical Bayesian model

IfQi,jdenotes the annual maximum flow for siteiin yearj,thenQi,jfollows a GEV distribution at sitei:

whereF（μi,σi,ξi）is a GEV distribution with location μi,scale σi,and shape ξi;i=1,2,…,I,withIbeing the total number of sites in the study area; andj= 1,2,…,Ni, withNibeing the recorded length for sitei.

Traditionally, the GEV distribution parameters and quantiles are estimated with the maximum likelihood estimation method,while their interval estimations are estimated with the delta method (Coles, 2001; Wu et al., 2019).However, the outputs of these methods are unstable,with high levels of bias and variance for small samples(Martins and Stedinger,2000).

To include more information in estimating parameters of the GEV model,in this study,we considered an HB model that incorporates spatial information (i.e., catchment characteristics)into model parameters through regression models.The model parameters are related to covariates through the monotonic link functiongk（·）:

wheregk（·）represents the one-to-one monotonic link function(e.g., identity function, and logarithm function); θk,iis thekth model parameter at sitei;hk（·） is the regression function,which can be linear or nonlinear; γkis a vector of regression coefficients;εk,iis the residual of regression;andxirepresents the set of covariates, which is related to the catchment characteristics (e.g., drainage area, elevation, and slope)of sitei.However, only the drainage area was considered in this study,and the extension framework of other characteristics was possible but was left for testing in the future.

Many studies have shown that there is a log-log linear relationship between the location and scale parameters and the drainage area (Gupta and Dawdy, 1995; Morrison and Smith,2002; Gaume et al., 2010; Lima and Lall, 2010; Villarini and Smith, 2010; Yan and Moradkhani, 2014; Lima et al., 2016).Therefore, Eq.(4)can be used to establish linear regression models between ln μ, ln σ, and lnA, where ln μ and ln σ correspond to the logarithm link function.Based on this knowledge,we considered independent normal distributions to be the prior distributions of parameters μ and σ:

where μiand σiare location and scale parameters of GEV at sitei;Aiis the drainage area of sitei;and α1, α2, β1, β2,,andare hyper-parameters.

In contrast to the transformed location and scale parameters,the shape parameter ξihas not been found to be dependent on the drainage areaAiin previous studies(Morrison and Smith,2002;Villarini and Smith, 2010).In general, the shape parameter is more difficult to estimate than the location and scale parameters(Yan and Moradkhani,2014),and the regional shape parameter is usually taken as an average of the local shape parameter at each site in many studies(Lima et al.,2016).Therefore,in this study,we considered a normal prior distribution for at-site shape parameters with a common mean across all sites:

The hyper-parameters in Eqs.(5)through (7)usually have no prior knowledge,and Gelman et al.(2014)suggested using non-informative independent prior distributions.As a result,the non-informative uniform priors in the range of（-∞,+∞）were selected for the hyper-parameters α1, α2, β1, β2,and,and the two-parameter inverse-gamma priors with parameters(0.01,0.01)were assigned for the hyper-parameters,and.

The Bayes'theorem was then used to calculate the posterior distribution of all the model parameters as follows:

whereqis the observed streamflow data set,θ=（μ,σ,ξ,α1,α2,β1,β2,）is the parameter vector in the model, Θ is the parameter space of parameter vector θ,p（θ） represents the prior distributions, andp（q|θ） is the likelihood function and has the following expression:

The joint posterior distribution of parameter vector θ can be obtained by substituting Eqs.(5)through (7)and (9)into Eq.(8):

whereqi,jis the observed annual maximum streamflow data for siteiin yearj.

The calculation of the posterior distribution in Eq.(10)is a challenging task because the dimensions increase with the number of sites.In this study, the Markov chain Monte Carlo(MCMC)method was employed to deal with this problem.This algorithm was implemented by using Gibbs sampling(Geman and Geman, 1984)in the software OpenBUGS (Lunn et al.,2009).Three parallel chains of a length of 2000 iterations were run for each model, and the first 1000 iterations were discarded as burn-in.The scale reduction factor ?R(Gelman and Rubin,1992)was used to test for convergence and the chain was regarded as having reached convergence if ?Rwas less than 1.1.

3.3.Local flood frequency analysis

Local(at-site)flood frequency analysis at a gauged site can be easily implemented with the proposed HB model.First,we draw samples of μi, σi, and ξifrom the joint posterior distributionp（θ|q） at sitei.Then, by substituting the samples of μi, σi, and ξiinto Eq.(2), we can obtain the posterior distribution of the quantile.Finally,local flood quantiles and their associated Bayesian credible intervals for different return periods can be estimated from the posterior distributions.

The predictive distribution of a future observation is also easily obtained through the posterior distribution in the HB model.Ifzdenotes the future flood design value,andp（θ|q）is the posterior distribution of θ for the observed valuesq,then the predictive density ofzcan be defined as follows(Coles,2001):

Compared with other predictive approaches, the predictive density based on the HB model has the advantage that it reflects uncertainty in the model through the posterior distributionp（θ|q） and uncertainty caused by the variability in future observations throughp（z|q） (Coles, 2001; Yoon et al., 2009).

In order to evaluate the performance of the HB model, the estimates of flood quantiles for different return periods and their associated uncertainties based on ordinary GEV distribution with the maximum likelihood estimate (MLE)at each site were used for comparison.

3.4.Regional flood frequency analysis

Flood quantile estimation at poorly gauged or ungauged sites usually depends on regional frequency analysis (RFA),which is performed by transferring the information arising from neighboring gauged sites (Hosking and Wallis, 1997).The index flood method proposed by Dalrymple(1960)is one of the most commonly used methods in engineering practice.For a pre-identified homogeneous region, the index flood method can be separated into the following four steps(Hosking and Wallis, 1997):(1)calculating the index flood λiby taking the mean or median of the annual maximum flow series for each site,then using the index flood λito standardize the annual maximum flow seriesQi,jof each site, i.e.,qi,j=Qi,j/λi; (2)estimating parameters of the common distribution(i.e., the GEV distribution in this study)at each site using the dimensionless rescaled dataqi,jwith L-moments; (3)building a regression equation linking the index flood λiand the catchment characteristics (i.e., drainage area in this study)from gauged sites for transferring the index flood λito an ungauged site;and(4)obtaining the quantile functionQi（F）of the frequency distribution at siteithroughQi（F） = λir（F）,wherer（F） is the regional growth curve of non-exceedance probabilityF.The limitations of the index flood method using the L-moments method are well-documented in the literature (Katz et al.,2002; Renard,2011;Lima et al.,2016).The HB model introduced in this study can overcome some of these issues.

For an ungauged or poorly gauged site, the predictive posterior distribution of μi, σi, and ξifrom Eqs.(5)through(7)is obtained first by using the drainage areaAiof the ungauged or poorly gauged site and the joint posterior distribution of the hyper-parameters α1, α2, β1, β2,, and.Then,the desired flood quantiles and their associated uncertainties can be estimated by substituting the obtained samples of μi, σi, and ξiinto Eq.(2).Compared with the traditional index flood method,the main advantage of prediction with the HB model is that it can predict the posterior distribution of some quantity of interest in a direct way.As a result,an uncertainty band of estimates can be quantified.For comparison, estimates of flood quantiles for different return periods based on the traditional index flood method using L-moments were also included for each selected site.

4.Results and discussion

4.1.Local flood frequency analysis

Fig.3.Logarithm of drainage area versus logarithm of GEV location parameter,logarithm of GEV scale parameter,and GEV shape parameter for 13 gauged sites (the black line indicates the ordinary least square regression).

Sites 4 and 11 (red triangles in Fig.1)were randomly chosen as cross-validation sites for the HB model regional flood frequency analysis, and were excluded from the modeling process.Hence,only the data from the remaining 13 gauged sites were used in the Bayesian inference.First, we drew scatter plots between the logarithm of MLEs of the GEV parameters and the logarithm of the drainage area.Then we fitted them with linear regression equations.It can be seen from Fig.3(a)and (b)that both the location μ and scale σ parameters of the GEV distribution show a well defined loglog linear relationship with the drainage areaA, and the coefficient of determinationR2of the regression models are 0.8702 and 0.6995, respectively.However, as expected, the shape parameter ξ does not present strong dependence on the drainage areaA(Fig.3(c)),with a negative slope of -0.6407,andR2is only 0.0453.The results indicate that the prior distributions of the three parameters of GEV in the HB model are reasonable.It is worth noting that there is an outlier on the top of Fig.3(c)(?ξ = 7.4031), which represents the MLE of the GEV shape parameter of site 6.The reason for this problem is that there are only five years of data available at this site,resulting in a poor MLE of the GEV shape parameter.

Fig.4(a)through (c)show the GEV parameter posterior density plots.Two gauged sites(sites 1 and 3)were randomly selected as examples.All the HB estimations (the median of the posterior distributions)and MLEs of GEV parameters for 13 gauged sites are shown in Table 2.

In order to compare the MLEs and HB model estimates of GEV parameters,Fig.5(a)through(c)show the boxplots of the parameters simulated by the HB model and the MLEs of GEV parameters (red stars).The Bayesian point estimates (the median of the posterior distributions)are represented by blue bars in the boxes.Notice that the distribution of parameters simulated by the HB model contains almost all the MLEs except for the scale parameter of site 8 (Fig.5(b))and the shape parameters of sites 6 and 8(Fig.5(c)),and that most of the MLEs are located in the box (between the first quantile and the third quantile).The reason for the MLE outliers at sites 6 and 8 is that data available from these two sites are limited (there are only five and 21 years,respectively)and cannot provide reliable maximum likelihood estimates.However,the HB estimates for the two sites that lie close to the average estimate from neighboring sites seem more plausible.In particular, the HB estimate=-0.0736 of site 6, is more reasonable than the MLE=7.4031.This shows the advantages of the HB model over the traditional MLE model for small samples.

Fig.6(a)through (d)show the estimation of flood quantiles and their associated uncertainties for four randomly chosen sites(sites 1,8,12,and 15).The black line is the flood quantile estimated by the HB model(i.e.,the median of posterior distribution)and the gray region is the associated 95%credible interval.The red solid line in the middle is the flood quantile estimated by the MLE method and the top and the bottom red dashed lines represent the 95%confidence interval based on the delta method.The empirical estimates are also shown by black dots for comparison.It can be seen from Fig.6(a)through(d)that the 95%confidence interval lower bound estimated by the delta method tends to decrease as the return period increases for each selected site, which indicates that the confidence interval for extreme events with long return periods is difficult to estimate with this method.However,the 95%credible interval estimated by the HB model is narrower than the corresponding confidence interval and still contains most of empirical estimates, highlighting the advantage of the HB model in uncertainty estimation.

The 100-year floods and their uncertainties estimated by the HB model and the traditional MLE method for all the 13 gauged sites are provided in Table 3.The MLEs of the GEV distribution parameters at site 6 are unreliable, leading to an unreasonably large quantile and large uncertainty,which is not shown in Table 3.The same problem also occurred at site 8,where the lower limit of the 100-year floods showed a highly negative value.Compared with the MLE method, the HB model shows that seven out of 13 sites(54%)show a reduction of uncertainty.The advantage of the HB model is that it can provide a more accurate description of parameter uncertainties and flood risk, and more realistic intervals of flood quantiles.This result corresponds to the previous study by Reis and Stedinger (2005).

Table 2 HB estimations (median of posterior distributions)and MLEs of GEV parameters for 13 gauged sites.

Fig.4.Posterior density plots of GEV parameters from HB model.

4.2.Regional flood frequency analysis

To test the validation of the HB model for regional flood frequency analysis, two sites (sites 4 and 11)were randomly selected as cross-validation sites, while the remaining 13 gauged sites were used to construct the HB model.The model prediction presented here is related to the problem of prediction in an ungauged basin(PUB)(Sivapalan,2003).The flood quantile predictions of the two cross-validation sites (sites 4 and 11)and the other two poorly gauged sites (sites 6 and 8)were predicted with the regional HB model.

Fig.6.Flood quantile estimates for four randomly chosen sites (sites 1, 8, 12, and 15)based on local HB model and MLE method.

Table 3 100-year floods and their uncertainties estimated by local HB model and MLE method for 13 gauged sites.

Fig.7(a)and (b)show the flood quantile prediction results for the two cross-validation sites(sites 4 and 11)based on the regional HB model.The black line represents the flood quantile estimated by the median of posterior distribution and the gray region represents the 95%credible interval for the HB model.For comparison, flood quantiles based on the index flood method with L-moments (the red dashed line)and the empirical estimates(black dots)are also shown in the figures.As we can see from Fig.7(a)and (b), all the empirical estimates (black dots)are distributed near the HB estimates and the index flood estimates, which indicates that the HB model can provide the same prediction accuracy as the index flood method with L-moments.However,the HB model can predict posterior distribution of flood quantiles, and as a result an uncertainty band of estimates can be quantified by the HB model,while the index flood method with L-moments does not demonstrate this uncertainty in a direct way.

The 100-year floods and their uncertainties estimated by the regional HB model and the index flood method for the two cross-validation sites (sites 4 and 11)and another two poorly gauged sites (sites 6 and 8)are listed in Table 4.As shown inthe table, the 95% credible interval of the flood predicted by the regional HB model contains the estimates of the index flood method at sites 6 and 8.Compared with the results of local analysis, the uncertainty estimation of the regional HB model is more reasonable than that of the delta method,which indicates that the prediction by the HB model at an ungauged or poorly gauged site is reliable.

Table 4 100-year floods and their uncertainties estimated by regional HB model and index flood method for two cross-validation sites(sites 4 and 11)and another two poorly gauged sites (sites 6 and 8).

5.Conclusions

(1)This study developed an HB model for local (at-site)and regional flood frequency analysis.The advantages of the HB model in uncertainty estimation of flood quantiles and prediction at poorly gauged or ungauged sites were shown in the Dongting Lake Basin.Integrating spatial information into GEV location and scale parameters with the drainage area was an effective way to improve the accuracy of parameter and flood quantile estimation, especially for the sites with short records and missing values.

(2)The HB model for local analysis can provide robust estimation of uncertainties for flood quantiles,while the traditional delta method cannot provide reliable uncertainty estimations for large flood quantiles, due to the fact that the lower confidence bound tends to decrease as the return periods increase.The HB model also showed satisfactory prediction results at ungauged sites.It predicted the posterior distribution of the desired flood quantile,and quantified uncertainty of estimates directly with an uncertainty band without any approximation, which the index flood method with L-moments could not do.

Fig.7.Flood quantile estimates for two cross-validation sites (sites 4 and 11)based on regional HB model and index flood method.

(3)The HB model for regional analysis, as an alternative approach of the traditional index flood method,improved upon many aspects,for example by relaxing the homogeneity region assumption and providing the total predictive uncertainty at ungauged sites.

Further work can be done in the future.Climate variables(e.g., El Ni~no Southern Oscillation (ENSO), and Pacific Decadal Oscillation)or other catchment characteristics (e.g.,elevation and slope)could be considered, or climate variables and catchment characteristics could be considered simultaneously in the GEV parameters in the model.Then,the model could be used to investigate the flooding risk of specific climate regimes (e.g., ENSO)as well as human interference.