Quantifying multi-source uncertainties in multi-model predictions using the Bayesian model averaging scheme

Hydrological models have been widely used in the watershed hydrological processes simulation and flood forecasting, and impact study of climate change and land-use change (Hailegeorgis & Alfredsen 2015; Emam et al. 2016; Jie et al. 2016). They play important roles in understanding the complex hydrologic cycle and solving practical hydrologic problems (Singh & Woolhiser 2002). Since the 1850s, hydrological models have experienced abundant development from empirical models through lumped conceptual models to physically based distributed models (Todini 2011). Nowadays, the precision of hydrological prediction has increased with the development of the model structure and improvement of the input data precision. However, in the hydrological processes simulation and flood forecasting, there still inevitably exist different modeling uncertainties, i.e., parameter uncertainty, input uncertainty, and model structural uncertainty (Beven 2001). Quantification and reduction of these uncertainties in hydrological modeling remain as challenges for hydrologists.

Numerous studies have recently focused on the itemized analysis of uncertainties of hydrological modeling (Krzysztofowicz 1999; Kavetski et al. 2006; Duan et al. 2007; McMillan et al. 2011; Dong et al. 2013; Liang et al. 2013; Yen et al. 2014a, 2015a, 2015b; Zhou et al. 2016). They highlighted that input error quantification, parameter optimization, and multi-model ensemble strategies are the three most popular methods used to demonstrate the impacts of hydrological prediction uncertainties. Rainfall is the most important input data for a hydrological model; thus, adequate characterization of rainfall is fundamental for the success of rainfall–runoff modeling. The true value of the amount of watershed rainfall in practice is often unknown because of its high spatial variability and insufficient rain gauge observations. Hence, an accurate statistical representation of watershed rainfall errors is critical for the estimation of uncertainty of rainfall inputs, which affect streamflow simulations. Kavetski et al. (2006) introduced a normally distributed error multiplier to reduce the precipitation input uncertainty. McMillan et al. (2011) evaluated the multiplicative error model of rainfall uncertainty and implied the dependence of rainfall error structure on the time-step data. Yen et al. (2015a) assessed the effects of the latent variables on the model simulations and implied the improvement of the model results is still limited. In hydrological modeling, model parameters often need to be calibrated based on observed hydrographs. Two main parameter calibration methods are currently used. In the first method, only one optimal parameter set can be obtained for a basin and model, and the typical algorithms are genetic algorithm (Wang 1991); shuffled complex evolution (SCE-UA, Duan et al. 1992) and dynamically dimensioned search (Tolson & Shoemaker 2007). In the other method, the model parameter involves one set of random variables that follow a certain joint probability distribution, and the typical algorithms are generalised likelihood uncertainty estimation (Beven & Binley 1992); shuffled complex evolution metropolis (SCEM-UA, Vrugt et al. 2003) and differential evolution adaptive metropolis (Vrugt et al. 2009). Different optimization algorithms demonstrate different convergence speed and behavioral statistics in model parameter calibration and uncertainty analysis (Xu et al. 2013; Yen et al. 2014a). Among the mentioned optimization algorithms, the SCE-UA and SCEM-UA approaches have been widely used in parameter calibration and uncertainty analysis in the literature, but the effects of the two algorithms on the deterministic simulation and probability prediction still need to be evaluated and compared further. This consideration has motivated our current study.

Different hydrological models have diverse foci in describing hydrological physical processes. No one model can sufficiently describe the principles of watershed rainfall–runoff in all conditions (Chen et al. 2013). Hence, an ensemble strategy based on multiple models has been considered as an effective method to reduce the uncertainty of model structures and improve the precision of hydrological predictions. Different model combination methods, such as neural network (Shamseldin et al. 1997), fuzzy system (Xiong et al. 2001), and Bayesian model averaging (BMA, Raftery et al. 2005), have emerged. Of these, BMA is the representative method that can consider the weighted average of the individual predictions from various models. It has been widely used in hydrological ensemble prediction studies. For example, Raftery et al. (2005) applied BMA to dynamic numerical weather predictions and attained valuable results. Duan et al. (2007), Liang et al. (2013), Dong et al. (2013), Yen et al. (2015b), Arsenault et al. (2015), and Zhou et al. (2016) successfully used BMA to combine multi-model/multi-method simulations to obtain more robust streamflow series and more reliable probability predictions. Jiang et al. (2012, 2014) also applied BMA to merge multi-satellite precipitation-based streamflow simulations to improve the hydrological utility of satellite precipitation products.

There are also some studies on assessment of the effects of different uncertainty sources on hydrological modeling (Kavetski et al. 2006; Ajami et al. 2007; Yen et al. 2014b). However, the comprehensive assessment of the effects of different uncertainty sources on deterministic simulation and probability prediction is still limited. Thus, the current study focuses on uncertainty analysis of multi-source and multi-model hydrological prediction. The innovations of the study include the following: (1) it considers rainfall input uncertainty, parameter estimation uncertainty, and model structural uncertainty by using three models, i.e., Xinanjiang (XAJ), hybrid rainfall–runoff (HYB), and HYMOD (HYM) models; (2) it compares the effects of SCE-UA and SCEM-UA algorithms on the hydrological prediction results; and (3) it investigates the superiority of the BMA multi-model ensemble strategy over the individual modeling approach.

The flowchart for the multi-source uncertainty analysis of multi-model predictions is shown in Figure 1. We adopted three different simulation cases to systematically consider the three sources (i.e., parameter uncertainty, input uncertainty, and model structural uncertainty) of hydrological modeling uncertainties. In case I, the model parameter uncertainty (hereafter ‘Para’) using SCE-UA and SCEM-UA algorithms for three hydrological models, i.e., XAJ, HYB, and HYM, was determined. In case II, a normally distributed error multiplier and combined parameter optimization algorithms were introduced to consider the model input and model parameter uncertainties (hereafter ‘Para + input’). In case III, the simulations calculated from case II were combined using BMA to comprehensively determine the model input, model parameter, and model structure uncertainties (hereafter ‘Para + input + struc’). The detailed methodologies are as follows.

XAJ model, hereinafter referred to as XAJ, is a well-known conceptual hydrological model developed by Zhao in the 1970s in China (Zhao 1992). In the present study, a sub-basin-structured semi-distributed XAJ model for streamflow simulation was constructed. The simulation was performed by computing the runoff from each sub-basin, and the slope and river network convergence processes were then integrated to obtain the streamflow series of the hydrologic station. A HYB model, hereinafter referred to as HYB, is a modified version of the XAJ model (Hu et al. 2005). Numerous field studies have shown that runoff within a basin is mainly generated by infiltration excess (Horton) runoff and saturation excess (Dunne) runoff (Ren et al. 2008). The HYB model combines the two runoff generation mechanisms by introducing spatial distribution curves of soil tension water storage capacity and infiltration capacity. Detailed descriptions of the mechanisms and applications of the HYB model were reported by Hu et al. (2005). HYMOD, hereinafter referred to as HYM, is a simple conceptual lumped hydrological model developed by Moore in the 1980s (Moore 1985). HYM consists of a simple rainfall excess model, which is connected to two series of linear reservoirs to route surface and subsurface flow. In the present study, an evaporation reduction factor K and a river network routing Muskingum–Cunge model were added to the original HYM. These three hydrological models have different complex model structures and different runoff generation mechanisms. They have been successfully and widely used in different river basins for streamflow simulation and flood forecasting (Ajami et al. 2007; Ren et al. 2008; Najafi et al. 2011; Jie et al. 2016; Xu et al. 2016). Tables 1–3 show the parameters and the prior ranges of the three models.

The models were operated on a daily time step within the 15 sub-basins in Mishui basin. The calibration period (CP) was from January 2000 to December 2005, and the period from January 2006 to December 2008 was used as the validation period (VP). This period of data was considered to be more representative of the current climate and land use situation of the study region.

SCE-UA is an effective and efficient global optimization algorithm proposed by Duan et al. (1992). It has been widely used in hydrological model parameter optimization. SCE-UA combines the direction searching of deterministic, non-numerical methods and the robustness of stochastic, non-numerical methods. It adopts the competition evolution theory, concepts of controlled random search, complex shuffling method, and downhill simplex procedures to obtain a global optimal estimation. Detailed calculation steps of SCE-UA are found in Duan et al. (1992).

SCEM-UA was built upon the principles of SCE-UA. Vrugt et al. (2003) combined the strengths of the Monte Carlo Markov Chain sampler with the concept of complex shuffling from SCE-UA to form the SCEM-UA algorithm, which not only provides the most probable parameter set, but also estimates the uncertainty associated with estimated parameters. SCEM-UA can simultaneously identify the most likely parameter set and its associated posterior probability distribution in every model run (Ajami et al. 2007). SCEM-UA has been successfully used in hydrologic and climate applications, such as rainfall–runoff model parameter calibration and uncertainty analysis (Ajami et al. 2007; Jiang et al. 2014). Detailed calculation steps of SCEM-UA are found in the work of Vrugt et al. (2003). In the present study, initial samples were obtained and then computations using SCEM-UA were performed using datasets with 5,000 and 10,000 samples.

BMA is a scheme for model combination that derives consensus predictions from competing predictions using likelihood measures as model weights. BMA has been primarily used to generalize linear regression applications. Raftery et al. (2005) successfully applied BMA to dynamic numerical weather predictions. Duan et al. (2007) and Ajami et al. (2007) used the BMA scheme to combine multiple models for hydrologic ensemble prediction that can provide more skillful and reliable predictions. The advantage of BMA is that the weights are directly bound with individual model simulation, that is, a well performing model can receive a higher weight than a poorly performing one. A more stable result can be obtained when the BMA method is used to combine different simulations. In the present study, we use BMA to merge the streamflow simulations from the three different hydrological models. Detailed calculation steps of the BMA method are found in the studies of Duan et al. (2007) and Ajami et al. (2007). For the sake of completeness, a brief description of the essence of the BMA scheme is presented as follows.

For SCE-UA-based simulation, the BMA weights and the variances of each model in the combination process were calculated, and then Monte Carlo Markov Chain sampling method was used to calculate the prediction uncertainty interval (Duan et al. 2007). Based on the repeated sampling experiments, we set the sampling times as 1,000. For SCEM-UA-based simulation, 15,000 streamflow series in the BMA combination process were simulated, and then normal population interval estimation method was used to calculate the prediction uncertainty interval (Ajami et al. 2007).

Mishui basin, a tributary of the Xiangjiang River, with a drainage area of 9,972 km² above the Ganxi hydrologic station, was selected as the study area (Figure 2). The basin is located southeast of Hunan Province in southern China and extends from longitudes 112.85 °E to 114.20 °E and latitudes 26.00 °N to 27.20 °N. The basin has a complex topography, with elevations ranging from 49 m to 2,093 m above sea level. The climate is of humid subtropical monsoon type, with an annual average temperature of approximately 18.0 °C and mean annual precipitation of approximately 1,561.0 mm. Temporal and spatial distributions of precipitation in the study region are uneven because of atmospheric circulation and most of the annual precipitation occurs between April and September. During these months, particularly in June, basin-wide heavy rains continuously occur, thereby resulting in flash floods. This multi-model ensemble prediction method can reduce the streamflow prediction and flood forecasting uncertainties, thus it is important to decision support systems for such river basins to prevent flood disasters and reduce flood damage.

The daily precipitation data from 2000 to 2008 were obtained from 35 rain gauge stations in the Mishui basin. For the same period, daily streamflow and potential evapotranspiration data were collected from the Ganxi hydrologic station and Wulipai evaporation station, respectively. This period of data was considered to be more representative of the current climate and land use situation of the study region. The inverse distance weighting of the three nearest rain gauges was used to obtain the spatially distributed precipitation database of 15 sub-basins for the Mishui basin. The 30 arc-second global digital elevation model data were obtained from the US Geological Survey. The vegetation-type data obtained from the International Geosphere-Biosphere Program were calculated and showed the land use distribution in the basin as forest and shrubs (54.4%), grasslands (33.5%), cropland (11.8%), and urban and water (0.3%).

The model parameters' prior ranges are defined in Tables 1–3, according to the physical meanings of the parameters and the actual hydro-climatic conditions of the Mishui basin. The SCE-UA algorithm gives a set of optimal solutions of the model parameters, while the SCEM-UA algorithm estimates the posteriori PDFs of the model parameters, which can reflect the effect of the model parameters' uncertainty on simulation result. The parameter frequency histograms were plotted by using the 10,000 group model parameters after the SCEM-UA algorithm convergence, in which the peak value of the posterior PDFs is the optimal parameter value for all samples. The marginal posterior probability distribution of the XAJ parameters estimated by SCEM-UA in case I is shown in Figure 3 and the statistical indices of the posterior probability distribution of the parameters estimated by SCEM-UA and the optimal parameters estimated by SCE-UA in case I are shown in Table 4. The histograms of XAJ parameters suggested that 12 parameters such as Kc, WDM, and so on (including all the sensitive parameters) approximately follow the normal distribution or the log-normal distribution. The rest of the parameters such as WLM and EX have two or more modal values, and this will increase the uncertainty of parameter optimization. Table 4 shows that the parameters WDM, EX, and CS0 have large CV values, implying that the mean value of the three parameters has poor representative power and large uncertainty. Some optimal parameters estimated by SCE-UA and SCEM-UA have some differences, and the possible reason may be due to the correlation between parameters and the ‘equifinality concept’ that different parameter sets may produce similar hydrologic behaviors (Beven & Binley 1992). Similar to the XAJ model results, most parameters of the HYB model and all parameters of the HYM model approximately follow the normal distribution or the log-normal distribution, which explains the effectiveness of the SCEM-UA optimization algorithm. Generally, the HYM model has a lesser number of parameters, which easily obey normal distribution. The XAJ and HYB models have more parameters than HYM model. For the influence of the correlation between parameters, the parameter uncertainty of XAJ model and HYB model is larger than that of the HYM model.

In order to consider the parameter and input uncertainty together, two rain input error modeling parameters m and are added to model parameter sets and further estimate the posterior PDFs simultaneously in case II. Figure 4 shows the marginal posterior probability distribution of the XAJ parameters estimated by SCEM-UA in case II. Table 5 demonstrates the statistical indices of the posterior probability distribution of the parameters estimated by SCEM-UA and the optimal parameters estimated by SCE-UA in case II. Comparing the parameter posterior PDFs of case II with that in case I, it can be concluded that the boundary of the models' parameters posterior distribution moves to a much more reasonable direction, and their posterior distributions are much closer to normal distribution. The rain input parameter is hard to concentrate to a single value, and it is difficult to optimize its value. This proved that there were rain input errors in the modeling, and the rain input error multiplier can describe the input errors to a certain extent. The two rain input parameters may introduce some new parameter estimating uncertainty and increase the difficulty of parameter optimization.

For comprehensive consideration of the model input, model parameter, and model structure uncertainties, we used the BMA to combine the three models' simulations for case II. Figure 5 displays the weight estimates of different models calculated using the BMA method. For the SCE-UA-based simulations, the weights of the XAJ, HYB, and HYM models are 0.36, 0.31, and 0.33, respectively. For the SCEM-UA-based simulations, the mean values of the weights of the XAJ, HYB, and HYM models are 0.35, 0.32, and 0.33, respectively. The weight of the BMA method is directly bound to individual model simulation, that is, a well performing model can receive a higher weight than a poorly performing one in theory. In this study, the XAJ model had the highest weight value, followed by the HYM model and the HYB model. The HYM model had a higher weight value than that of the HYB model, which may be due to the similar model structure of the XAJ model and the HYB model (Ren et al. 2008). By using the BMA combination, we can obtain deterministic streamflow series and probability predictions, which are comprehensively considered the multi-source uncertainties.

Table 6 shows the statistical performances of the streamflow simulations based on the SCE-UA and SCEM-UA algorithms of the three simulation cases (in which the value set in bold font refers to the optimum performance in the column). Figures 6–8 show the BMA combined streamflow series from the SCE-UA-based simulations and the SCEM-UA-based simulations of the three simulation cases, respectively. From Table 6 and Figures 6 and 7, we can see that the three models showed a good hydrologic prediction applicability in the Mishui basin, in which the XAJ model performed best, followed by the HYB model, and lastly, the HYM model. Especially for the high flow simulations, the XAJ model and the HYB model performed much better than the HYM model simulation. Generally, both parameter optimization algorithms generated good and comparative streamflow simulations. The SCEM-UA implied parameter uncertainty and provided the posterior distribution of the parameters. Using the 15,000 simulation sets, SCEM-UA showed a certain advantage over the SCE-UA algorithm in the calculation of the prediction uncertainty bounds. Given the precipitation input uncertainty in case II, the precisions of the simulated streamflows using the three models were not remarkably enhanced. This phenomenon may have been caused by the relatively small precipitation input uncertainty due to the dense rain gauge observations in the Mishui basin. Moreover, in the model parameters, an evaporation reduction factor parameter K was set, and this parameter could imply some precipitation input uncertainty. Our results are quite consistent with those of Yen et al. (2015a), who reported that the use of error multiplier to incorporate input uncertainty might not be the proper alternative choice in terms of generating better results. In case III, for both the SCE-UA and SCEM-UA algorithms, BMA combinations of the simulation sets improved the precision of streamflow predictions, especially during the VP. This condition was indicated by the high NSE and the small BIAS and RMSE values from BMA combinations compared with those from each single model (see Table 6). The daily NSE, BIAS, and RMSE values of the SCE-UA-based BMA combination in case III for the CP were 0.91, 0.04%, and 35.99 m³/s, respectively; and the corresponding values for the VP were 0.88, 3.85%, and 56.32 m³/s. The daily NSE, BIAS, and RMSE values of the SCEM-UA-based BMA combination in case III for the CP were 0.92, 0.16%, and 34.66 m³/s, respectively; and the corresponding values for the VP were 0.87, 3.49%, and 59.93 m³/s. Using BMA in combining multiple models to conduct ensemble streamflow simulation can effectively improve the precision of streamflow simulations, especially for the VP.

Table 7 shows the reliability performance of the calculated 95% confidence interval of the three simulation cases. Figures 6–8 show the 95% confidence interval from the SCE-UA-based simulations (sampling done 1,000 times) and from the SCEM-UA-based simulations of the three simulation cases, respectively. Both parameter optimization algorithms generated a certain precision of prediction uncertainty interval. However, the 95% confidence interval of the SCEM-UA-based simulation was much better than that of the SCE-UA-based simulation. With higher CR and lower D values, the SCEM-UA algorithm had an advantage in the estimation of prediction uncertainty bounds compared with the SCE-UA algorithm. Given the precipitation input uncertainty in case II, the performance of the calculated 95% confidence intervals of the three models showed minimal improvement in terms of higher CR values, especially for the VP. In case III, for both the SCE-UA and SCEM-UA algorithms, the reliability performance of the 95% confidence interval calculated from the BMA combined streamflows was much better than the performance of the interval from each signal model (see Table 7). The daily CR, B, and D values of the SCE-UA-based BMA combination for the CP were 90.19%, 315.60 m³/s, and 56.70 m³/s, respectively; and the corresponding values for the VP were 90.97%, 348.56 m³/s, and 69.74 m³/s. The daily CR, B, and D values of the SCEM-UA-based BMA combination for the CP were 95.62%, 271.15 m³/s, and 55.03 m³/s, respectively; and the corresponding values for the VP were 95.17%, 303.04 m³/s, and 66.06 m³/s. The calculated 95% confidence interval from the BMA combination had higher CR and better D values than those of each single model. In addition, it also had a higher B value. The increase in the uncertainty interval CR value was accompanied by an increase in B value, which has already been discussed by Xiong et al. (2009) and Dong et al. (2013). Thus, using BMA in combining multiple models to perform the ensemble hydrologic simulations can effectively calculate more reliable uncertainty bounds.

Figure 8 compares the BMA-combined streamflow mean values and the calculated 95% confidence interval with the observed hydrograph at the daily time scales from the SCE-UA-based simulation and SCEM-UA-based simulation for case III. Both SCE-UA- and SCEM-UA-based BMA combinations generated good streamflow simulations and reliable 95% confidence intervals. The precisions of streamflow simulations of the SCE-UA- and SCEM-UA-based simulations were comparatively good, but the reliability of SCEM-UA-calculated 95% confidence interval was much better than that of SCE-UA in terms of higher CR and lower B and D values (Table 7). Figure 8 also demonstrates that the SCEM-UA-calculated 95% confidence interval can preferably contain the observed high flows and this is very important for flood control decision-making. For the low flow series, the SCEM-UA-based method can give much better confidence intervals than that of the SCE-UA-based method. Thus, the results suggest that determining the model parameter uncertainties using the SCEM-UA algorithm can generate more reliable uncertainty bounds than that of the simulation from SCE-UA.

For the hydrological simulation and forecast, the hydrological model must go through the model parameters' calibration and validation stages. The hydrological model can be applied to practical use only on the condition that the calibrated model can also perform well in the VP (Singh & Woolhiser 2002). Different hydrological models have different instabilities in the CP and VP for their variant climatic conditions, respectively (Yan et al. 2013; Li et al. 2015). Most models cannot have the same performance in the VP as that in the CP. Figure 9 compares the hydrological models' simulation performances in the CP and VP for the three different simulation cases. Figure 10 shows the distribution of RMSE for XAJ, HYB, and HYM considering different uncertainty sources. From Figures 9 and 10, we can see that both at case I and case II, the three hydrological models have better simulation precision in the CP than that in the VP. While, at case III, by using the BMA combination of three hydrological models, it can effectively improve the precision of streamflow predictions in terms of high NSE value, small BIAS and RMSE values in the VP. Normally, hydrological modeling has higher uncertainties in the VP than in the CP, while the BMA multi-model ensemble strategy can effectively improve this phenomenon and give a higher skill and reliability forecasting for the future (Velázquez et al. 2011; Broderick et al. 2016). Thus, choosing appropriate hydrological models, considering the parameter uncertainties, and using the multi-model ensemble strategy, can improve the accuracy of the hydrological forecasting results.

This study performed a multi-source uncertainty analysis of hydrological prediction by using input error quantification, parameter optimization, and multi-model ensemble methods in a typical humid watershed in southern China. The results show that both the SCE-UA and SCEM-UA parameter optimization algorithms can make the XAJ, HYB, and HYM models generate good streamflow simulations with NSE values higher than 0.80 and BIAS values smaller than 7.62%. Specifically, the SCEM-UA can imply parameter uncertainty and provide the posterior distribution of the parameters. Thus, the SCEM-UA algorithm has advantages in the estimation of model parameter uncertainty and predicting reliable hydrological forecasts. Considering precipitation input uncertainty does not improve the precision of streamflow simulation in the selected Mishui basin, which is probably due to the availability of good quality and dense rain gauge stations. The BMA combination of the simulation sets calculated from single models not only improves the precision of streamflow predictions in terms of NSE and BIAS values, but also quantifies the uncertainty bounds for the simulation sets in terms of CR values. The improvement of the prediction precision of the BMA combination is much more evident in the VP than in the CP. This finding demonstrates that the hydrological modeling has more uncertainties in the VP, and that the BMA multi-model ensemble can effectively reduce these uncertainties. Comparison of the prediction uncertainty interval from the two different parameter optimization algorithms shows that the calculated 95% prediction interval from SCEM-UA-based BMA simulations is much better than that calculated from SCE-UA-based BMA simulations. Hence, these results suggest that the comprehensive uncertainty analysis concerning model parameters uncertainties and multi-model ensembles by using the SCEM-UA algorithm and BMA method is advantageous and of practical importance for streamflow predictions and flood forecasting, which can collectively provide more robust streamflow series and more reliable uncertainty bounds.

The current study was jointly supported by the National Key Research and Development Program approved by Ministry of Science and Technology, China (2016YFA0601504); the Programme of Introducing Talents of Discipline to Universities by the Ministry of Education and the State Administration of Foreign Experts Affairs, China (B08048); the National Natural Science Foundation of China (41501017, 51579066) and Natural Science Foundation of Jiangsu Province (BK20150815).