Probabilistic prediction in ungauged basins (PUB) based on regional parameter estimation and Bayesian model averaging

Predictions in ungauged basins (PUB) are widely considered as an important and challenging research topic in the hydrological sciences (Sivapalan et al. 2003; Hrachowitz et al. 2013). One of the primary research objectives of the PUB initiative was to improve the ability of existing hydrological models to predict in ungauged basins with reducing uncertainties (Merz & Bloschl 2004; Young 2006; Heuvelmans et al. 2006; Wagener & Wood 2006; Jin et al. 2009; Dong et al. 2013a, 2013b; Parkes et al. 2013). Regionalization of parameters is diffusely used to simulate runoff in PUB, which is regarded as the process of transferring parameter values from a donor gauged basin to the target ungauged basin (Vandewiele & Elias 1995; Xu 1999; Xie et al. 2007; Parajka et al. 2007; Bastola et al. 2008; Kizza et al. 2013; Hailegeorgis & Alfredsen 2015). Four regionalization approaches have been typically used for choosing the donor gauged basin whose calibrated and optimized parameter values are transferred to simulate runoff for the target ungauged basin: regression analysis, spatial proximity, physical similarity, and a mixture of them.

Various studies have been performed to determine which is the best approach between the spatial proximity method and the physical similarity approach. Post & Jakeman (1999) investigated the relationships between the model parameters of a lumped conceptual rainfall-runoff model and the basin landscape attributes of similarly sized basins. Burn & Boorman (1993), Johansson (1994), Sefton & Howarth (1998) and Kokkonen et al. (2003) derived the relationships between model parameters and physical catchment descriptor indices using geographical information system. Croke et al. (2004) adopted a simple hydrologic approach to simulate runoff adaption to land-use changes in ungauged basins. Goswami et al. (2007) developed a pooling method of regional parameter estimation (RPE) coupled with soil data and SMAR model in order to simulate flow in ungauged basins in France. Zhang & Chiew (2009) evaluated the disadvantages and advantages of different regionalization methods using two rainfall runoff models, Xinanjiang and SIMHYD in 210 Australian basins. The study showed that the best approach between the spatial proximity method and the physical similarity approach is hard to identify. Most studies suggested that the use of more information (such as remotely sensed (RS) vegetation data, soil data, climate, and land cover) or a mixture of them can improve the accuracy of runoff simulation in ungauged basins (Ao et al. 2006; Kay et al. 2006; Duan et al. 2006; Götzinger & Bárdossy 2007; Oudin et al. 2008; Bulygina et al. 2009; Li et al. 2009; Kling & Gupta 2009; Samaniego et al. 2011; Li et al. 2014; Ghahraman & Davary 2014).

Previous studies of PUB mostly focused on the comparison of prediction of individual parameter transfer schemes and their weight averages. Our research aims to compare the probabilistic prediction generated by the Bayesian model averaging (BMA) with that of each individual parameter transfer scheme, in order to see if BMA can effectively reduce the uncertainty in runoff PUB and improve the prediction reliability.

The paper is organized as follows: a brief introduction to the study area; then a general description of the main steps and procedures including: (1) variable infiltration capacity-three layers (VIC-3L) model description, (2) regional parameter transfer schemes, (3) probabilistic prediction using BMA, and (4) performance evaluation; comparison between BMA and its three individual parameter transfer schemes is then discussed, finishing with the conclusions drawn from this study.

The Hanjiang River is the largest tributary of the Yangtze River and it passes through the provinces of Shannxi and Hubei in China, and merges into the Yangtze River at Wuhan city. The river is of length 1,570 km and area 159,000 km². The basin has a sub-tropical monsoon climate and has, as a result, a dramatic diversity in its water resources. Annual precipitation varies from 700 to 1,100 mm, with 70–80% of the total amount occurring in the wet season from May to October. The Hanjiang basin plays a critical role in flood control and water supply in central China. The Danjiangkou reservoir located in the middle reach of the Hanjiang River is the source of water for the middle route of the South-North Water Division Project (SNWDP), and the Jianghan plain in the down basin is one of the most important bases for commodity grain production.

A procedure coupling a regional parameter transfer module with probabilistic prediction module is developed to obtain probabilistic prediction in ungauged basins. It consists of two modules, namely a regional parameter transfer module and probabilistic prediction. The regional parameter transfer module is aimed at obtaining model parameter estimates from a limited number of calibrated basins and then regionalizing them to uncalibrated basins based on the spatial proximity approach, physical similarity approach (similar characteristics of climate as well as physicality), and multiple regression analysis, which is described in detail in the next section. Probabilistic prediction is designed to infer a prediction by weight averaging over many different regional parameter transfer schemes based on the BMA method, which is described in detail in the section ‘Hydrological probability prediction’. The main steps and procedures include: (1) VIC-3L model description; (2) regional parameter transfer scheme; (3) probabilistic prediction using BMA; and (4) performance evaluation.

The VIC-3L model has one kind of bare soil and different vegetation types in each grid cell (Liang & Xie 2001; Xie et al. 2003). It includes both the saturation and infiltration excess runoff processes in a grid cell with a consider­ation of the sub-grid scale soil heterogeneity, and the frozen soil processes for cold climate conditions. The one-dimensional Richard equation is used to describe the vertical soil moisture movement and the moisture transfer between soil layers obeys the Darcy law. The ARNO method is used to describe base flow which takes place only in the lowest layer. The routing model represented by the unit hydrograph method for overland flow and the linear Saint-Venant method for channel flow, allows runoff to be predicted (Liang et al. 1994).

The VIC-3L model has 10 hydrological parameters that need to be calibrated, as shown in Table 1. A similar description for VIC-3L parameters was made by Xie et al. (2007).

For a detailed description of methods for weather forcing data, vegetation dataset, and soil dataset based on regionalization and grid in this paper, readers are referred to Xie et al. (2007).

BMA is a statistical technique designed to infer a prediction by weight averaging over many different regional parameter transfer schemes. This method is not only a pathway for scheme combination but also a coherent approach for accounting for between-scheme and within-scheme uncertainty (Ajami et al. 2007). Below is a brief description of the basic ideas of this method.

To estimate BMA weight and scheme prediction variance ⁠, the expectation-maximization (EM) algorithm, which has proved to be an efficient technique for BMA calculation based on the assumption that K-member predictions are normally distributed, is described in this section (Duan et al. 2007).

It is difficult to maximize the function (6) by analytical methods. The EM algorithm is an effective method for finding the maximum likelihood by alternating between two steps, the expectation step and maximization step. The two steps are iterated to convergence when there is no significant change between two consecutive iterative log-likelihood estimations. In the EM algorithm, a latent variable (unobserved quantity) is used as an assistant for estimating BMA weight ⁠. For a detailed description of the EM algorithm for a BMA scheme, readers are referred to Dong et al. (2013b).

After estimating BMA weight and prediction variance ⁠, we use the Monte Carlo method to generate BMA probabilistic prediction for any time t (Hammersleym & Handscomb 1975). The procedures are described as follows.

1. Generate an integer value of k from [1, 2, …, K] with probability ⁠. A specific procedure is described as follows.

   • 1(a). Set the cumulative weight and compute for k = 1, 2, …, K.

   • 1(b). Generate a random number u between 0 and 1.

   • 1(c). If ⁠, this indicates that we choose the kth member of the ensemble predictions.

2. Generate a value of from the PDF of ⁠. Here, represents the normal distribution with mean and variance ⁠.

3. Repeat the above steps (1) and (2) for M times. M is the probabilistic ensemble size. In this paper, we set M = 100.

After generating the BMA probabilistic ensemble prediction, results are sorted in ascending order. From this, the 90% uncertainty intervals can be derived within the range of the 5% and 95% quantities.

For each individual scheme in the BMA model, the prediction uncertainty interval can also be constructed, with the Monte Carlo sampling method still being used to approximate the assumed PDF of ⁠.

There are three indices for evaluating the mean prediction (Dong et al. 2013b) presented as follows.

1. The Nash-Sutcliffe coefficient of efficiency (NS) – NS is not only an objective function but also a widely used performance criterion. It ranges from minus infinity to 1.0, with higher values indicating better agreement. The definition of NS is expressed in following equation.

   

   7

   where and are observed and simulated data at time t, is the average of observed data, T is the length of the data series.
2. Daily root mean square error (DRMS):

   

   8

1. Relative error of total runoff (RE):

   

   9

Xiong et al. (2009) and Dong et al. (2013b) presented a set of indices for assessing the probabilistic prediction generated by the uncertainty analysis methods. Three main indices are selected here to assess the probabilistic prediction produced by the BMA model as well as from each individual parameter transfer scheme.

1. Containing ratio (CR) – The containing ratio is used for assessing the goodness of the uncertainty interval. It is defined as the percentage of observed data points that are covered in the prediction bounds.

   

   10

   where and denote as upper and lower prediction bounds at time t, and is the number of observed data points that are covered in the prediction bounds.
2. Average band-width (B). Consider

   

   11

   where the average band-width B is also an index for measuring the performance of estimated uncertainty interval.
3. Average deviation amplitude (D) – The average deviation amplitude is an index to quantify the average deflection of the curve of the middle points of the prediction bounds from the observed streamflow hydrograph. It is defined as

   

   12

   with notations as defined previously.

Daily streamflow and weather data from 1980 to 1986 are used for calibration. The gauged sub-basins are selected as the primary basins to implement VIC-3L model calibration, which is achieved by matching the total annual streamflow volume and the shape of the mean daily hydrograph to the corresponding observations in the Ankang River basin. The two criteria, i.e., NS and RE, are used for model calibration. In the calibration study, the parameters of individual sub-basins with similar climate characteristics and underlying surface are assumed to have the same values. Ten hydrological parameters in the VIC-3L model have been calibrated for the 10 primary sub-basins. Table 3 shows the calibrated parameter values in 10 primary sub-basins in Ankang basin. Their typical ranges and the effect of each parameter on results of simulated streamflow are described below. (1) b typically ranges from 0 to 0.50. It describes the total of available infiltration capacity as a function of the relative saturated grid cell area and controls the quantity of runoff generation directly and the water balance. A lower value of b gives lower infiltration and yields higher surface runoff (the value of b in this paper is the inverse value of b in Xie et al. 2007). The highest value of b in sub-basins is only 0.450 for Lanhe sub-basin and the lowest value is 0.157 for Zhehe sub-basin. The rest of the values of b in sub-basins are very close to 0.257, because Ankang basin is in a humid region. (2) D_m typically ranges from 0 to 6 mm day⁻¹. (3) D_s typically ranges from 0 to 1. With a higher value of D_s, the base flow will be higher at lower water content in the lowest soil layer. (4) W_s typically ranges from 0 to 1. (5) Dep₁, Dep₂, and Dep₃ range from 0 to 0.40 m. In general, thicker soil depths slow down seasonal peak flows and increase the loss due to evapotranspiration. (6) x ranges from 0.05 to 0.30. (7) k ranges typically from 0.50 to 2.0. (8) ckg ranges typically from 0.5 to 1.0.

Table 4 lists the statistical results (NS, DRMS, and RE) for the 10 primary sub-basins using calibrated hydrological parameters. In terms of NS, DRMS and RE, the model in the calibration provides good simulation results for all sub-basins. In the next section, we will demonstrate that these calibrated hydrological parameters can be transferred to the ungauged sub-basins with reasonably good results.

In this paper, the BMA method is used to predict a new measurement value associated with a combination of probabilistic prediction in ungauged basins based on three individual parameter transfer schemes. The comparison between BMA (3) and its three individual parameter transfer schemes is made in terms of both mean prediction and probabilistic prediction in this study. The main conclusions are summarized as follows: (1) the mean prediction of BMA (3) is much closer to the observed data as compared with its best individual parameter transfer scheme (physical similarity approach) for two ungauged sub-basins; and (2) the probabilistic predictions of BMA (3) have larger containing ratio, larger average band-width, and smaller average deviation amplitude than any of its individual parameter transfer schemes for the two ungauged sub-basins. It is worth mentioning that further works will focus on transferring parameter approaches based on data mining and machine learning techniques, such as artificial neural networks and support vector machine, as well as choosing other hydrological models, then it is anticipated that the advantages of BMA can be generalized.

This study is financially supported by the International Cooperation in Science and Technology Special Project of China (2014DFA71910), National Natural Science Foundation of China (51509008, 51509141 and 51379223), Natural Science Foundation of Hubei Province (2015CFB217) and Open Foundation of State Key Laboratory of Water Resources and Hydropower Engineering Science in Wuhan University (2014SWG02).