## Abstract

This study focuses on a quantitative multi-source uncertainty analysis of multi-model predictions. Three widely used hydrological models, i.e., Xinanjiang (XAJ), hybrid rainfall–runoff (HYB), and HYMOD (HYM), were calibrated by two parameter optimization algorithms, namely, shuffled complex evolution (SCE-UA) method and shuffled complex evolution metropolis (SCEM-UA) method on the Mishui basin, south China. The input uncertainty was quantified by utilizing a normally distributed error multiplier. The ensemble simulation sets calculated from the three models were combined using the Bayesian model averaging (BMA) method. Results indicate the following. (1) Both SCE-UA and SCEM-UA resulted in good and comparable streamflow simulations. Specifically, the SCEM-UA implied parameter uncertainty and provided the posterior distribution of the parameters. (2) In terms of the precipitation input uncertainty, precision of streamflow simulations did not improve remarkably. (3) The BMA combination not only improved the precision of streamflow prediction, but also quantified the uncertainty bounds of the simulation. (4) The prediction interval calculated using the SCEM-UA-based BMA combination approach appears superior to that calculated using the SCE-UA-based BMA combination for both high flows and low flows. Results suggest that the comprehensive uncertainty analysis by using the SCEM-UA algorithm and BMA method is superior for streamflow predictions and flood forecasting.

## INTRODUCTION

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.

## METHODOLOGY

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.

### Hydrological models

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.

Parameter | Physical meaning | Prior range |
---|---|---|

Kc | Ratio of potential evapotranspiration to pan evaporation | 0.5–1.5 |

WUM | Water capacity in the upper soil layer | 10–40 |

WLM | Water capacity in the lower soil layer | 50–90 |

WDM | Water capacity in the deeper soil layer | 10–70 |

B | Exponent of the tension water capacity curve | 0.1–0.5 |

C | Coefficient of deep evapotranspiration | 0.1–0.3 |

EX | Exponent of the free water capacity curve | 1–1.5 |

SM | The free water capacity of the surface soil layer | 10–60 |

KI0 | Outflow coefficients of the free water storage to interflow | KI + KG = 0.7 |

KG0 | Outflow coefficients of the free water storage to groundwater | 0.1–0.5 |

CI0 | Recession constant of the lower interflow storage | 0.1–0.9 |

CG0 | Daily recession constant of groundwater storage | 0.9–0.999 |

CS0 | Recession constant for channel routing | 0.1–0.5 |

KE | Slot storage coefficient | 20–24 |

XE | Flow proportion factor | 0.1–0.5 |

Parameter | Physical meaning | Prior range |
---|---|---|

Kc | Ratio of potential evapotranspiration to pan evaporation | 0.5–1.5 |

WUM | Water capacity in the upper soil layer | 10–40 |

WLM | Water capacity in the lower soil layer | 50–90 |

WDM | Water capacity in the deeper soil layer | 10–70 |

B | Exponent of the tension water capacity curve | 0.1–0.5 |

C | Coefficient of deep evapotranspiration | 0.1–0.3 |

EX | Exponent of the free water capacity curve | 1–1.5 |

SM | The free water capacity of the surface soil layer | 10–60 |

KI0 | Outflow coefficients of the free water storage to interflow | KI + KG = 0.7 |

KG0 | Outflow coefficients of the free water storage to groundwater | 0.1–0.5 |

CI0 | Recession constant of the lower interflow storage | 0.1–0.9 |

CG0 | Daily recession constant of groundwater storage | 0.9–0.999 |

CS0 | Recession constant for channel routing | 0.1–0.5 |

KE | Slot storage coefficient | 20–24 |

XE | Flow proportion factor | 0.1–0.5 |

Parameter | Physical meaning | Prior range |
---|---|---|

Kc | Ratio of potential evapotranspiration to pan evaporation | 0.5–1.5 |

WUM | Water capacity in the upper soil layer | 10–40 |

WLM | Water capacity in the lower soil layer | 50–90 |

WDM | Water capacity in the deeper soil layer | 10–70 |

B | Exponent of the tension water capacity curve | 0.1–0.5 |

bx | Infiltration capacity distribution curve index | 0.1–2 |

f0 | The average maximum infiltration capacity | 5–30 |

fc | The average stability infiltration capacity | 0.1–10 |

k | Infiltration capacity attenuation coefficient | 0.1–0.9 |

CS | Recession constant for channel routing | 0.1–0.5 |

CG | Daily recession constant of groundwater storage | 0.9–0.999 |

C | Coefficient of deep evapotranspiration | 0.1–0.3 |

KE | Slot storage coefficient | 20–24 |

XE | Flow proportion factor | 0.1–0.5 |

Parameter | Physical meaning | Prior range |
---|---|---|

Kc | Ratio of potential evapotranspiration to pan evaporation | 0.5–1.5 |

WUM | Water capacity in the upper soil layer | 10–40 |

WLM | Water capacity in the lower soil layer | 50–90 |

WDM | Water capacity in the deeper soil layer | 10–70 |

B | Exponent of the tension water capacity curve | 0.1–0.5 |

bx | Infiltration capacity distribution curve index | 0.1–2 |

f0 | The average maximum infiltration capacity | 5–30 |

fc | The average stability infiltration capacity | 0.1–10 |

k | Infiltration capacity attenuation coefficient | 0.1–0.9 |

CS | Recession constant for channel routing | 0.1–0.5 |

CG | Daily recession constant of groundwater storage | 0.9–0.999 |

C | Coefficient of deep evapotranspiration | 0.1–0.3 |

KE | Slot storage coefficient | 20–24 |

XE | Flow proportion factor | 0.1–0.5 |

Parameter | Physical meaning | Prior range |
---|---|---|

Kc | Ratio of potential evapotranspiration to pan evaporation | 0.5–1.5 |

Cmax | Max height of soil moisture accounting tank | 1–1,000 |

bexp | Distribution function shape parameter | 0.1–2 |

Alpha | Quick-slow split parameter | 0.1–0.99 |

Nq | Number of quick-flow routing tanks | 1–8 |

Rs | Slowflow routing tanks rate parameter | 0.001–0.1 |

Rq | Quick-flow routing tanks rate parameter | 0.1–0.99 |

KE | Slot storage coefficient | 20–24 |

XE | Flow proportion factor | 0.1–0.5 |

Parameter | Physical meaning | Prior range |
---|---|---|

Kc | Ratio of potential evapotranspiration to pan evaporation | 0.5–1.5 |

Cmax | Max height of soil moisture accounting tank | 1–1,000 |

bexp | Distribution function shape parameter | 0.1–2 |

Alpha | Quick-slow split parameter | 0.1–0.99 |

Nq | Number of quick-flow routing tanks | 1–8 |

Rs | Slowflow routing tanks rate parameter | 0.001–0.1 |

Rq | Quick-flow routing tanks rate parameter | 0.1–0.99 |

KE | Slot storage coefficient | 20–24 |

XE | Flow proportion factor | 0.1–0.5 |

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.

### Input error modeling

*et al.*2007). In this study, we adopted an error multiplier to determine the precipitation input uncertainty: where and are the measured and modified precipitation at time step t, respectively, is a normal error multiplier with a mean value of

*m*and a variance of at time step t. Based on the research of Ajami

*et al.*(2007), we assume that and .

### Parameter optimization

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

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.

*y*is BMA prediction, are observed data sets (in which

*X*denotes input forcing data and

*Y*is observed streamflow data) and is the ensemble of the

*K*-member predictions. The posterior distribution of the BMA prediction

*y*is given as: where is the posterior probability of the prediction given the input data

*D*, and it reflects how well model fits

*Y*. Actually is the BMA weight , and better performing predictions receive higher weights than the worse performing ones, and all weights are positive and should add up to 1. is the conditional probability density function (PDF) of the prediction

*y*conditional on and

*D*. Thus, the posterior mean and variance of the BMA prediction could be expressed as: where is the variance associated with model prediction with respect to observation

*D*. Compared with the deterministic multi-model combination method, BMA can better describe the uncertainty of analog variable. In this study, we use the expectation-maximization algorithm to estimate the BMA weight and model prediction variance (Ajami

*et al.*2007).

### Prediction uncertainty interval

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

### Evaluation statistics

*n*is the number of simulation days.

*et al.*2009). CR, expressed as percentage, denotes the ratio of the number of observed streamflows enveloped by prediction bounds to the total number of observed hydrographs. B represents the average bandwidth of the whole prediction bounds. With a certain confidence level, a lower B value indicates a better prediction bound. D denotes the actual discrepancy between the trajectories consisting of the middle points of the prediction bounds and the observed hydrograph. It also shows the symmetry with respect to the observed discharges and the middle point of the prediction bounds. The formulas for CR, B, and D are given as: where is the number of observed streamflows enveloped by prediction bounds,

*n*is the total number of observed hydrographs, and and are the upper and low boundaries of the prediction bounds at time step i, respectively.

## STUDY AREA AND DATA

### Study area

Mishui basin, a tributary of the Xiangjiang River, with a drainage area of 9,972 km^{2} 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.

### Data

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

## RESULTS AND DISCUSSION

### Parameter uncertainty analysis

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.

Parameter | Kc | WUM | WLM | WDM | B | C | EX |
---|---|---|---|---|---|---|---|

Mean | 1.32 | 39.43 | 80.76 | 31.17 | 0.49 | 0.28 | 1.31 |

SD | 0.03 | 0.67 | 6.03 | 8.22 | 0.01 | 0.02 | 0.19 |

CV | 0.02 | 0.02 | 0.07 | 0.26 | 0.01 | 0.05 | 0.15 |

SCE-UA | 1.49 | 39.99 | 50.01 | 10.06 | 0.47 | 0.20 | 1.42 |

SCEM-UA | 1.34 | 39.21 | 86.06 | 40.20 | 0.49 | 0.29 | 1.46 |

Parameter | SM | KG0 | CI0 | CG0 | CS0 | KE | XE |

Mean | 24.08 | 0.38 | 0.84 | 0.99 | 0.13 | 20.17 | 0.50 |

SD | 1.53 | 0.02 | 0.02 | 0.01 | 0.03 | 0.27 | 0.01 |

CV | 0.06 | 0.05 | 0.02 | 0.00 | 0.23 | 0.01 | 0.00 |

SCE-UA | 36.64 | 0.47 | 0.81 | 0.99 | 0.16 | 20.08 | 0.50 |

SCEM-UA | 25.22 | 0.35 | 0.85 | 0.99 | 0.10 | 20.01 | 0.50 |

Parameter | Kc | WUM | WLM | WDM | B | C | EX |
---|---|---|---|---|---|---|---|

Mean | 1.32 | 39.43 | 80.76 | 31.17 | 0.49 | 0.28 | 1.31 |

SD | 0.03 | 0.67 | 6.03 | 8.22 | 0.01 | 0.02 | 0.19 |

CV | 0.02 | 0.02 | 0.07 | 0.26 | 0.01 | 0.05 | 0.15 |

SCE-UA | 1.49 | 39.99 | 50.01 | 10.06 | 0.47 | 0.20 | 1.42 |

SCEM-UA | 1.34 | 39.21 | 86.06 | 40.20 | 0.49 | 0.29 | 1.46 |

Parameter | SM | KG0 | CI0 | CG0 | CS0 | KE | XE |

Mean | 24.08 | 0.38 | 0.84 | 0.99 | 0.13 | 20.17 | 0.50 |

SD | 1.53 | 0.02 | 0.02 | 0.01 | 0.03 | 0.27 | 0.01 |

CV | 0.06 | 0.05 | 0.02 | 0.00 | 0.23 | 0.01 | 0.00 |

SCE-UA | 36.64 | 0.47 | 0.81 | 0.99 | 0.16 | 20.08 | 0.50 |

SCEM-UA | 25.22 | 0.35 | 0.85 | 0.99 | 0.10 | 20.01 | 0.50 |

*Notes:* In the table, SD indicates standard deviation, CV means variable coefficient, SCE-UA and SCEM-UA mean the optimal parameter values of the two algorithms, respectively.

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.

Parameter | Kc | WUM | WLM | WDM | B | C | EX | SM |
---|---|---|---|---|---|---|---|---|

Mean | 1.17 | 39.38 | 75.38 | 40.84 | 0.48 | 0.29 | 1.30 | 22.70 |

SD | 0.03 | 0.55 | 9.61 | 9.22 | 0.02 | 0.01 | 0.19 | 0.83 |

CV | 0.03 | 0.01 | 0.13 | 0.23 | 0.05 | 0.03 | 0.14 | 0.04 |

SCE-UA | 1.35 | 36.59 | 67.72 | 54.05 | 0.50 | 0.17 | 1.13 | 21.15 |

SCEM-UA | 1.16 | 39.14 | 85.56 | 28.32 | 0.50 | 0.29 | 1.43 | 22.59 |

Parameter | KG0 | CI0 | CG0 | CS0 | KE | XE | a | v |

Mean | 0.38 | 0.82 | 0.99 | 0.11 | 20.07 | 0.50 | 0.954 | 0.0004 |

SD | 0.04 | 0.05 | 0.00 | 0.01 | 0.07 | 0.00 | 0.004 | 0.0003 |

CV | 0.10 | 0.06 | 0.00 | 0.09 | 0.00 | 0.00 | 0.004 | 0.5800 |

SCE-UA | 0.34 | 0.86 | 0.99 | 0.12 | 20.00 | 0.50 | 0.950 | 0.0002 |

SCEM-UA | 0.35 | 0.86 | 0.99 | 0.11 | 20.02 | 0.50 | 0.951 | 0.0003 |

Parameter | Kc | WUM | WLM | WDM | B | C | EX | SM |
---|---|---|---|---|---|---|---|---|

Mean | 1.17 | 39.38 | 75.38 | 40.84 | 0.48 | 0.29 | 1.30 | 22.70 |

SD | 0.03 | 0.55 | 9.61 | 9.22 | 0.02 | 0.01 | 0.19 | 0.83 |

CV | 0.03 | 0.01 | 0.13 | 0.23 | 0.05 | 0.03 | 0.14 | 0.04 |

SCE-UA | 1.35 | 36.59 | 67.72 | 54.05 | 0.50 | 0.17 | 1.13 | 21.15 |

SCEM-UA | 1.16 | 39.14 | 85.56 | 28.32 | 0.50 | 0.29 | 1.43 | 22.59 |

Parameter | KG0 | CI0 | CG0 | CS0 | KE | XE | a | v |

Mean | 0.38 | 0.82 | 0.99 | 0.11 | 20.07 | 0.50 | 0.954 | 0.0004 |

SD | 0.04 | 0.05 | 0.00 | 0.01 | 0.07 | 0.00 | 0.004 | 0.0003 |

CV | 0.10 | 0.06 | 0.00 | 0.09 | 0.00 | 0.00 | 0.004 | 0.5800 |

SCE-UA | 0.34 | 0.86 | 0.99 | 0.12 | 20.00 | 0.50 | 0.950 | 0.0002 |

SCEM-UA | 0.35 | 0.86 | 0.99 | 0.11 | 20.02 | 0.50 | 0.951 | 0.0003 |

*Notes:* In the table, SD indicates standard deviation, CV means variable coefficient, SCE-UA and SCEM-UA mean the optimal parameter values of the two algorithms, respectively.

### Streamflow comparison between BMA ensemble and single model

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^{3}/s, respectively; and the corresponding values for the VP were 0.88, 3.85%, and 56.32 m^{3}/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^{3}/s, respectively; and the corresponding values for the VP were 0.87, 3.49%, and 59.93 m^{3}/s. Using BMA in combining multiple models to conduct ensemble streamflow simulation can effectively improve the precision of streamflow simulations, especially for the VP.

Cases | SCE-UA | SCEM-UA | ||||
---|---|---|---|---|---|---|

NSE | BIAS (%) | RMSE (m^{3}/s) | NSE | BIAS (%) | RMSE (m^{3}/s) | |

CP | ||||||

XAJ (Para) | 0.91 | −2.36 | 37.05 | 0.92 | 3.13 | 34.68 |

XAJ (Para + input) | 0.90 | 4.37 | 37.58 | 0.92 | 2.23 | 34.05 |

HYB (Para) | 0.88 | 2.50 | 42.49 | 0.89 | −1.08 | 39.41 |

HYB (Para + input) | 0.87 | −6.00 | 42.53 | 0.88 | −3.41 | 41.27 |

HYM (Para) | 0.85 | 1.38 | 46.51 | 0.85 | 1.31 | 46.63 |

HYM (Para + input) | 0.85 | 1.17 | 46.79 | 0.85 | 1.67 | 46.69 |

BMA (Para + input + struc) | 0.91 | 0.04 | 35.99 | 0.92 | 0.16 | 34.66 |

VP | ||||||

XAJ (Para) | 0.83 | 1.90 | 69.23 | 0.81 | 6.14 | 71.95 |

XAJ (Para + input) | 0.85 | 7.62 | 64.03 | 0.82 | 5.12 | 70.23 |

HYB (Para) | 0.80 | 6.64 | 74.37 | 0.82 | 3.35 | 69.35 |

HYB (Para + input) | 0.86 | −2.50 | 62.26 | 0.83 | − 0.70 | 67.04 |

HYM (Para) | 0.83 | 6.25 | 69.10 | 0.83 | 6.17 | 69.19 |

HYM (Para + input) | 0.82 | 5.79 | 69.42 | 0.83 | 6.26 | 69.08 |

BMA (Para + input + struc) | 0.8 | 3.85 | 56.32 | 0.87 | 3.49 | 59.93 |

Cases | SCE-UA | SCEM-UA | ||||
---|---|---|---|---|---|---|

NSE | BIAS (%) | RMSE (m^{3}/s) | NSE | BIAS (%) | RMSE (m^{3}/s) | |

CP | ||||||

XAJ (Para) | 0.91 | −2.36 | 37.05 | 0.92 | 3.13 | 34.68 |

XAJ (Para + input) | 0.90 | 4.37 | 37.58 | 0.92 | 2.23 | 34.05 |

HYB (Para) | 0.88 | 2.50 | 42.49 | 0.89 | −1.08 | 39.41 |

HYB (Para + input) | 0.87 | −6.00 | 42.53 | 0.88 | −3.41 | 41.27 |

HYM (Para) | 0.85 | 1.38 | 46.51 | 0.85 | 1.31 | 46.63 |

HYM (Para + input) | 0.85 | 1.17 | 46.79 | 0.85 | 1.67 | 46.69 |

BMA (Para + input + struc) | 0.91 | 0.04 | 35.99 | 0.92 | 0.16 | 34.66 |

VP | ||||||

XAJ (Para) | 0.83 | 1.90 | 69.23 | 0.81 | 6.14 | 71.95 |

XAJ (Para + input) | 0.85 | 7.62 | 64.03 | 0.82 | 5.12 | 70.23 |

HYB (Para) | 0.80 | 6.64 | 74.37 | 0.82 | 3.35 | 69.35 |

HYB (Para + input) | 0.86 | −2.50 | 62.26 | 0.83 | − 0.70 | 67.04 |

HYM (Para) | 0.83 | 6.25 | 69.10 | 0.83 | 6.17 | 69.19 |

HYM (Para + input) | 0.82 | 5.79 | 69.42 | 0.83 | 6.26 | 69.08 |

BMA (Para + input + struc) | 0.8 | 3.85 | 56.32 | 0.87 | 3.49 | 59.93 |

*Notes:* In the table, Para indicates considering model parameter uncertainty in case I, Para + input means considering model input and model parameter uncertainties in case II, Para + input + struc denotes considering model input, model parameter, and model structure uncertainties in case III. The value set in bold font refers to the optimum performance in the column.

### Prediction interval comparison between BMA ensemble and single model

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^{3}/s, and 56.70 m^{3}/s, respectively; and the corresponding values for the VP were 90.97%, 348.56 m^{3}/s, and 69.74 m^{3}/s. The daily CR, B, and D values of the SCEM-UA-based BMA combination for the CP were 95.62%, 271.15 m^{3}/s, and 55.03 m^{3}/s, respectively; and the corresponding values for the VP were 95.17%, 303.04 m^{3}/s, and 66.06 m^{3}/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.

Cases | SCE-UA (sampling 1,000 times) | SCEM-UA | ||||
---|---|---|---|---|---|---|

CR% | B (m^{3}/s) | D (m^{3}/s) | CR% | B (m^{3}/s) | D (m^{3}/s) | |

CP | ||||||

XAJ (Para) | 59.31 | 152.87 | 58.20 | 78.65 | 169.17 | 52.31 |

XAJ (Para + input) | 74.86 | 200.15 | 60.78 | 79.06 | 169.78 | 51.79 |

HYB (Para) | 75.05 | 258.66 | 74.08 | 80.34 | 222.25 | 64.41 |

HYB (Para + input) | 81.07 | 273.01 | 65.39 | 78.97 | 225.23 | 67.60 |

HYM (Para) | 71.40 | 225.70 | 63.08 | 85.26 | 237.97 | 62.91 |

HYM (Para + input) | 68.57 | 212.49 | 63.29 | 87.68 | 254.06 | 64.01 |

BMA (Para + input + struc) | 90.19 | 315.60 | 56.70 | 95.62 | 271.15 | 55.03 |

VP | ||||||

XAJ (Para) | 62.32 | 183.64 | 71.21 | 80.47 | 188.68 | 64.41 |

XAJ (Para + input) | 73.81 | 220.50 | 74.63 | 81.48 | 190.31 | 63.07 |

HYB (Para) | 71.99 | 289.95 | 82.84 | 80.66 | 244.40 | 74.13 |

HYB (Para + input) | 82.66 | 285.44 | 71.71 | 80.38 | 249.24 | 77.14 |

HYM (Para) | 68.61 | 270.26 | 77.88 | 86.77 | 261.23 | 76.62 |

HYM (Para + input) | 69.16 | 252.23 | 76.84 | 88.96 | 278.24 | 77.31 |

BMA (Para + input + struc) | 90.97 | 348.56 | 69.74 | 95.17 | 303.04 | 66.06 |

Cases | SCE-UA (sampling 1,000 times) | SCEM-UA | ||||
---|---|---|---|---|---|---|

CR% | B (m^{3}/s) | D (m^{3}/s) | CR% | B (m^{3}/s) | D (m^{3}/s) | |

CP | ||||||

XAJ (Para) | 59.31 | 152.87 | 58.20 | 78.65 | 169.17 | 52.31 |

XAJ (Para + input) | 74.86 | 200.15 | 60.78 | 79.06 | 169.78 | 51.79 |

HYB (Para) | 75.05 | 258.66 | 74.08 | 80.34 | 222.25 | 64.41 |

HYB (Para + input) | 81.07 | 273.01 | 65.39 | 78.97 | 225.23 | 67.60 |

HYM (Para) | 71.40 | 225.70 | 63.08 | 85.26 | 237.97 | 62.91 |

HYM (Para + input) | 68.57 | 212.49 | 63.29 | 87.68 | 254.06 | 64.01 |

BMA (Para + input + struc) | 90.19 | 315.60 | 56.70 | 95.62 | 271.15 | 55.03 |

VP | ||||||

XAJ (Para) | 62.32 | 183.64 | 71.21 | 80.47 | 188.68 | 64.41 |

XAJ (Para + input) | 73.81 | 220.50 | 74.63 | 81.48 | 190.31 | 63.07 |

HYB (Para) | 71.99 | 289.95 | 82.84 | 80.66 | 244.40 | 74.13 |

HYB (Para + input) | 82.66 | 285.44 | 71.71 | 80.38 | 249.24 | 77.14 |

HYM (Para) | 68.61 | 270.26 | 77.88 | 86.77 | 261.23 | 76.62 |

HYM (Para + input) | 69.16 | 252.23 | 76.84 | 88.96 | 278.24 | 77.31 |

BMA (Para + input + struc) | 90.97 | 348.56 | 69.74 | 95.17 | 303.04 | 66.06 |

*Notes:* In the table, Para indicates considering model parameter uncertainty in case I, Para + input means considering model input and model parameter uncertainties in case II, Para + input + struc denotes considering model input, model parameter, and model structure uncertainties in case III. The value set in bold font refers to the optimum performance in the column.

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.

### The different performance in CP and VP

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.

## CONCLUSIONS AND SUGGESTIONS

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.

## ACKNOWLEDGEMENTS

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