Evaluation of the HYMOD model for rainfall – runoff simulation using the GLUE method

The Yalong River is the third largest base of the 13 hydropower bases in China. Long-time series of river discharge records are essential for the design of hydropower stations and water resource management. The existing monitoring network is scarce and cannot provide sufficient hydrological information for the basin. Rainfall–runoff models are popular tools for extending hydrological data in both space and time. In this paper, the feasibility of applying a conceptual rainfall–runoff model, HYdrological MODel (HYMOD), to the upper Yalong River basin was evaluated. The generalized likelihood uncertainty estimation (GLUE) was employed for model calibration and uncertainty analysis. The results show that simulated discharge matches the observations satisfactorily, indicating the hydrological model performs well and the application of HYMOD to estimate long time-series of river discharge in the study area is feasible.


INTRODUCTION
Located in western Sichuan Province and southeastern Qinghai-Tibet Plateau, the Yalong River is one of the biggest tributaries of the Jinsha River. The Yalong has a total length of 1500 km and annual runoff volume of about 61 billion m³. Possessing abundant hydropower resources, the basin is among the 13 largest hydropower bases of China. Until 2025, 21 levels of hydroelectric power stations will be designed and developed. Long time-series of river discharge data are essential for the design of hydropower stations, water resource management and ecological effects assessment. At present, river discharge monitoring stations within the upstream region of Yalong River basin are scarce (Zhang et al., 2000). The hydrological data is insufficient for hydrological analysis, which is important for the design of hydropower stations and operation rules.
Rainfall-runoff models are widely applied to estimate river discharge. Researchers have widely applied the distributed hydrological model in the study of the river basin hydrological cycle and obtained good results in simulations of hydrological features and runoff estimation (Wu et al., 2007;Zanon et al., 2013;Zhao et al., 2013). At present, several studies have been carried out to analyse the hydrological features of Yalong River basin. Yu et al. (2007) applied a distributed hydrological model, the WEP model to the water cycle in the Yalong River basin. For the upper stream region of the basin, Chen et al. (2001) estimated characteristics of river discharge at the monthly scale based on the HBV conceptual model.
The research of estimating river discharge at the daily scale in this region is scarce. In this paper, a daily step conceptual rainfall-runoff model, HYMOD, is applied in the upstream region of Yalong River basin for the purpose of evaluating its feasibility for estimating river discharge of long time series.

STUDY AREA
In this paper, the upper reaches of the Yalong River basin above the Ganzi in situ gauging station is selected as the study area. This area has a drainage area of about 33 000 km 2 with a river length of about 690 km. The dominant climate in this basin is continental climate. Most of the rainfall events happen between May and October, which account for 90% to 95% of total annual rainfall, varying from 500 to 650 mm. Most part of the upstream region of Yalong River is inaccessible to humans due to the high elevation, which makes the construction of hydrological stations difficult.
There is only one national hydrological station in the entire upstream of Yalong River basin located in Ganzi prefecture.

METHODOLOGY
(a) HYMOD and data for the model HYMOD model is a parsimonious daily step model with typical conceptual hydrological components, based on the theory of runoff yield under excess infiltration. The runoff generation process is described by a simple rainfall excess model based on the probability-distributed principle (Moore, 1985).The model structure is illustrated in Fig. 1. To describe the spatial variability of water storage capacity within the basin, the following cumulative distribution function is used: where, F is the cumulative rate of storage capacity of some point of the basin; C is the water storage capacity; Cmax and Bexp are two parameters describing basin maximum water storage capacity (mm) and the degree of spatial variability within the basin, respectively. The portion of precipitation that exceeds the water storage capacity is treated as runoff. The evapotranspiration is equal to potential evapotranspiration (PET) if enough water is available. Otherwise, it equals the available water storage. Based on the parameter Alpha, the runoff is divided into quick flow and slow flow, which are routed through three identical quick flow tanks (Q1, Q2, Q3) and a parallel slow flow tank, respectively. The flow rates in the routing system are described by the resident time in the quick tanks Kq (day) and the slow tank Ks (day), respectively. The parameters related to runoff generation (Cmax, Bexp and Alpha) were assumed to be identical in every sub-basin. The two routing parameters of each sub-basin describe the process of water flow movement from the sub-basin to the outlet of the whole basin. The travel time to the basin outlet varies in different sub-basins. For each time step, the output of the revised HYMOD is the sum of the water flow from each sub-basin, which arrives at the outlet of the basin in that step.
The input forcing data include basin average rainfall and PET. In this paper, 8 years (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987) of rainfall data of each sub-basin is selected from four national weather stations (Ganzi, Seda, Shiqu, Qingshui). The PET data of each sub-basin were obtained from a global 30-min grid PET data set: Ahn and Tateishi (1994) PET. In order to account for the spatial heterogeneity of the meteorological condition in the target basin, the study area was divided into four subbasins. The observed runoff data (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987) from the Ganzi in situ gauging station were used for calibration and validation. The calibration period is from 1980 to 1983 and the validation period is 1984-1987.

(b) The GLUE uncertainty analysis method
The GLUE (Generalized Likelihood Uncertainty Estimation) method is the most widely used and efficient method applied in studying the hydrological model uncertainty. The method is a global parameters sensitivity analysis method developed based on the importance sampling and regionalized sensitivity analysis of the RSA (Beven et al., 1992). GLUE distinguishes different parameters using likelihood. This method takes into account the facts that the optimal is the best, and avoids the risks of using a single optimal value to predict as well. The GLUE method is similar to the RSA method in employing the sensitivity analysis of comparing the distribution of the behavioural parameter sets and the original parameter sets. It can also use the cumulative likelihood for global sensitivity analysis (Blasone et al., 2008). The key point of the GLUE method is that the quality of the model simulation results is not determined by a single parameter but a set of model parameters of the model. Therefore, the GLUE method uses the Monte Carlo method to randomly sample the parameter values from the prior distribution within a given initial parameter distribution space. Then we form multi-group sets of parameter combinations to run the model (Beven et al., 2001) In this paper, the GLUE method is applied to analyse the uncertainty of the HYMOD model with the sequential daily runoff data (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987) from the Ganzi in situ gauging station. For the case study, five parameters are included in the model calibration. Using a Latin-Hypercube sampling algorithm, 50 000 parameter sets are randomly selected based on the parameter ranges listed in Table 1 using the Monte Carlo method. The prior distribution of parameters is assumed to have a uniform distribution. To assess the performance of each parameter set, the Nash-Sutcliffe (NS) coefficient is chosen as the objective function. After calculating the NS coefficient values for each parameter set, the likelihood threshold for rejecting the non-behavioural parameter sets is subjectively decided as 0.6 (Sun et al., 2010). Parameter sets with NS coefficient values larger than 0.6 are chosen as behavioral parameter sets. The cumulative probability distribution of discharge is carried out with the simulated discharge data derived by the behavioural parameter sets. Finally, the 5% and 95% quantile of the cumulative distribution are used as the lower and upper limit of the uncertainty band.
To evaluate the simulated results and analyse the model uncertainty, the P-factor and R-factor are calculated. The P-factor is the percentage of the measured data bracketed by the 95% predicted uncertainty. If the P-factor is bigger and closer to 1, it indicates that the uncertainty interval brackets more measured runoff data and the simulation result is better. The R-factor is the ratio of the average width of uncertainty intervals and the standard deviation of measured runoff data. A larger P-factor and a smaller R-factor indicates that the uncertainty interval is narrower and brackets more measured runoff data and thus the simulation results are better.

RESULTS AND DISCUSSION (a) Sensitivity analysis of parameters
The Monte Carlo method was employed to randomly generate 50 000 parameter sets. Then each parameter set was applied to run the model. Among the 50 000 parameter sets, the NS coefficient of simulated river discharge of 1953 sets of parameters reaches 0.6. These sets of parameters were identified as behavioural ones and were used to analyse the parameter uncertainties and for the ensemble simulation.
Although the GLUE method treats the combination of parameter values as a whole parameter set, the sensitivity of each model parameter can be inferred from the scatterplot (Fig. 2) of identified behavioural values of every parameter. Among the five model parameter, only the Cmax shows a single peak distribution. The values of the other four parameters span the whole parameter ranges, which is a common representation of equifinality. The prior and posterior distributions of the model parameters were also analysed (Fig. 3). If the posterior distribution deviates from the prior distribution (i.e. the uniform distribution) significantly, the parameter is considered as a sensitive parameter (Peters et al., 2003;Beven et al., 2001;Sun et al., 2012). The result shows that posterior distributions of Cmax, Bexp and Kq, are quite different from prior distribution, indicating that these parameters are more sensitive than the other two (Ks and Alpha).

(b) Results of the runoff modelling and uncertainty analysis
The simulated river discharge and the value of efficiency are shown in Fig. 4 and Table 2, respectively. The NS coefficients of the parameter sets with best simulation of the calibration and validation periods are 0.75 and 0.74, respectively. Judging from the hydrograph, the simulated river discharge could reproduce the observations reasonably. The scatterplot (Fig. 5) of the NS coefficients of the behavioural parameters demonstrates that a high NS coefficient in the calibration period is accompanied by a high NS coefficient in the validation period, indicating that the parameter sets performing well in calibration period also perform well in the validation period.
Regarding the simulation uncertainty issue, the P-factors in the calibration and validation period are 0.40 and 0.58, i.e. 40% of the observed data fall into the uncertainty interval of the ensemble simulation in the calibration period, while in the validation period, it is 58%. The Rfactors in the two periods are about 1.3. Comparing the observed and simulated daily runoff hydrograph (Fig. 4) of Ganzi gauging station shows that for most time of a year, the simulation matches observations well. However, in the winter and snow melting seasons, there is some deviation between the observations and simulations. There are two possible reasons: first, for such a large simulation area (about 33 000 km 2 ), precipitation data from only four rainfall stations are available for simulation. Secondly, HYMOD is simple in its structure and the snow melting process is not considered. Generally, even using such a limited number of precipitation data, the simulation results are acceptable. It is concluded that HYMOD is applicable for runoff simulation in the upstream region of the Yalong River basin.

CONCLUSIONS
In this study, the HYMOD model was applied to the upper reaches of the Yalong River basin, above the Ganzi in situ gauging station, and the GLUE method was used for model calibration and uncertainty analysis. The results show that the parameters Cmax, Bexp and Kq of the HYMOD model are relatively sensitive. The river discharge simulation fit the in situ observation well. Utilization of limited precipitation data and the fact that snow melting processes are not considered in the model structure are two possible reasons for the relatively low performance of the model simulation in winter and snow melting periods. In general, the model is feasible for rainfall-runoff process simulation of the study area. With more rainfall data, HYMOD is expected to perform better and be useful to estimate long time-series of river discharge at the daily scale in the upstream region of Yalong River basin.