Analysis of parameter uncertainty in hydrological and sediment modeling using GLUE method : a case study of SWAT model applied to Three Gorges Reservoir Region , China 1

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26  Z.Y. Shen∗, L. Chen, T. Chen State Key Laboratory of Water Environment Simulation, School of Environment, Beijing Normal University, Beijing 100875, P.R. China E-mail address:zyshen@bnu.edu.cn Tel./fax: +86 10 58800398. Abstract The calibration of hydrologic models is a worldwide challenge due to the uncertainty involved in the large number of parameters. The difficulty even increases in a region with high seasonal variation of precipitation, where the results exhibit high heteroscedasticity and autocorrelation. In this study, the Generalized Likelihood Uncertainty Estimation (GLUE) method was combined with the Soil and Water Assessment Tool (SWAT) to quantify the parameter uncertainty of the stream flow and sediment simulation in the Daning River Watershed of the Three Gorges Reservoir Region (TGRA), China. Based on this study, only a few parameters affected the final simulation output significantly. The results showed that sediment simulation presented greater uncertainty than stream flow, and uncertainty even greater in high precipitation conditions than during the dry season. The main uncertainty sources of stream flow came from the catchment process while a channel process impacts the sediment simulation greatly. It should be noted that identifiable parameters such as CANMX, ALPHA_BNK, SOL_K could be obtained with an optimal parameter range using calibration method. However, equifinality was also observed in hydrologic modeling in TGRA. This study demonstrated that care must be taken when calibrating the SWAT model with non-identifiable parameters because these may lead  * Corresponding author. Tel./fax: +86 10 58800398.


Introduction
Watershed hydrology and river water quality models are important tools for watershed management for both operational and research programs (Quilbé and Rousseau, 2007; Van et al., 2008;Sudheer and Lakshmi, 2011).However, due to spatial variability in the processes, many of the physical models are highly complex and generally characterized by a multitude of parameters (Xuan et al., 2009).Technically, the modification of parameter values reveals a high degree of uncertainty.Overestimation of uncertainty may lead to expenditures in time and money and overdesign of watershed management.Conversely, underestimation of uncertainty may result in little impact on pollution abatement (Zhang et al., 2009).In order to apply hydrological models in the practical water resource investigations, careful calibration and uncertainty analysis are required (Beven and Binley, 1992;Vrugt et al., 2003;Yang et al., 2008).
Much attention has been paid to uncertainty issues in hydrological modeling due to their great effects on prediction and further on decision-making (Van et al., 2008;Sudheer and Lakshmi, 2011).Uncertainty estimates are routinely incorporated into Total Maximum Daily Load (TMDL) (Quilbé and Rousseau, 2007).Usually, the uncertainty in hydrological modeling is from model structures, input data and parameters (Lindenschmidt et al., 2007).In general, structural uncertainty could be improved by comparing and modifying the diverse model components (Hejberg and Refsguard, 2005).The uncertainty of model input occurs because of changes in natural conditions, limitations in measurement, Published by Copernicus Publications on behalf of the European Geosciences Union.and lack of data (Berk, 1987).One way to deal with this issue is to use random variables as the input data, rather than the conventional form of fixed values (Yulianti et al., 1999).Currently, parameter uncertainty is a hot topic in the uncertainty research field (Shen et al., 2008;Sudheer et al., 2011).
The model parameters can be divided into the conceptual group and the physical group (Gong et al., 2011).The conceptual parameters such as CN 2 in the SCS curve method are defined as the conceptualization of a non-quantifiable process, and determined by the process of model calibration.Conversely, physical parameters can be measured or estimated based on watershed characteristics when intensive data collection is possible (Vertessy et al., 1993;Nandakumar and Mein, 1997).Because of the unknown spatial heterogeneity of a studied area and the expensive experiments which may be involved, the physical parameters are usually determined by calibrating the model against the measured data (Raat et al., 2004).However, when the number of parameters is large either due to the large number of sub-processes being considered or due to the model structure itself, the calibration process becomes complex and uncertainty issues appear (Rosso, 1994;Sorooshian and Gupta, 1995).It has been shown that parameter uncertainty is inevitable in hydrological modeling and a corresponding assessment should be conducted before model prediction in the decision making process.Studies of parameter uncertainty have been conducted in the area of integrated watershed management (Zacharias et al., 2005), peak flow forecasting (Jorgeson and Julien, 2005), soil loss prediction (Cochrane and Flanagan, 2005), nutrient flux analysis (Murdoch et al., 2005;Miller et al., 2006), assessment of the effect of land use change (Eckhardt et al., 2003;Shen et al., 2010;Xu et al., 2011) and climate change impact assessment (Kingston and Taylor, 2010), among many others.Nevertheless, parameter identification is a complex, non-linear problem and numerous possible solutions might be obtained by optimization algorithms (Nandakumar and Mein, 1997).Thus, the parameters cannot be identified easily.Additionally, different parameter sets may result in similar prediction which is known as the phenomenon of equifinality (Beven and Binley, 1992).However, to the best of our knowledge, there are few studies about parameter identifiability based on uncertainty analysis in hydrological modeling.
Several calibration and uncertainty analysis techniques have been applied in previous research work, such as the first-order error analysis (FOEA) (Melching and Yoon, 1996), the Monte Carlo method (Kao and Hong, 1996) and the Generalized Likelihood Uncertainty Estimation method (GLUE) (Beven and Binley, 1992).The FOEA method is based on linear-relationships and fails to deal adequately with the complex models (Melching and Yoon, 1996).The Monte Carlo method requires repeating model simulation according to the parameter sampling, resulting in tremendous computational time and human effort (Gong et al., 2011).However, the GLUE methodology determines the performance of the model focus on the parameter set, not on the individual parameters (Beven and Binley, 1992).The GLUE method can also handle the parameter interactions and non-linearity implicitly through the likelihood measure (Vazquz et al., 2009).In addition, GLUE is a simple concept and is relatively easy to implement.Therefore, GLUE is used in this study for parameter uncertainty analysis.
The Three Gorges Project -the largest hydropower project in the world -is situated at Sandoupin in Yichang City, Hubei Province, China.It is composed mainly of the dam, the hydropower station, the two-lane, five-stage navigation locks, and the single-lane vertical ship lift.While the Three Gorges Project provides flood control, power generation, and navigation benefits, it also has a profound impact on the hydrology and environment, such as river flow interruption and ecosystem degradation.Hydrological models have been used in this region to study the impact of the project (Lu and Higgitt, 2001;Yang et al., 2002;Wang et al., 2007;Shen et al., 2010).However, research on the uncertainty of hydrological models in such an important watershed is lacking.Due to the varying geographical locations and water systems (Xu et al., 2011), it is of great importance to study the uncertainty of model parameters that affect the hydrological modeling process.Previously we had conducted a parameter uncertainty analysis for nonpoint source pollution modeling in this region.In the present investigation, a further study in hydrological modeling was developed.
Hence, the main objective of this study was to identify the degree of uncertainty and uncertainty parameters for prediction of stream flow and sediment in a typical watershed of the Three Gorges Reservoir Region, China.In this study, a semi-distributed hydrological model, Soil and Water Assessment tool (SWAT) was combined with the GLUE (Generalized likelihood uncertainty estimation) method to quantify the uncertainty of parameters and to provide a necessary reference for hydrological modeling in the entire Three Gorges Reservoir region.
The paper is organized as follows: (1) a description of the study area and a brief introduction of the hydrological model and GLUE method; (2) both the impact of parameter uncertainty on model results and parameter identifiability are analyzed in the result and discussion section; (3) a conclusion is provided.

Site description
The Daning River Watershed (108 • 44 -110 • 11 E, 31 • 04 -31 • 44 N), lies in the central part of the Three Gorges Reservoir Area (TGRA) (Fig. 1), is in Wushan and Wuxi Counties, in the municipality of Chongqing, China and covers an area of 4426 km 2 .Mountainous terrain makes up 95 % of the total area and low hills contribute the other 5 %.The average altitude is 1197 m.The landuse in the watershed is 22.2 % cropland, 11.4 % grassland, and 65.8 % forest.Zonal yellow soil is the dominant soil of the watershed.This area is characterized by the tropical monsoon and subtropical climates of Northern Asia.A humid subtropical monsoon climate covers this area, featuring distinct seasons with adequate sunshine (an annual mean temperature of 16.6 • C) and abundant precipitation (an annual mean precipitation of 1124.5 mm).A hydrological station is located in Wuxi County, and this study focused on the watershed controlled by the Wuxi hydrological station, which has an area of approximately 2027 km 2 (Fig. 1).

SWAT model
The SWAT model (Arnold et al., 1998) is a hydrologic/water quality tool developed by the United States Department of Agriculture-Agriculture Research Service (USDAARS).The SWAT model is also available within the BASINS (Better Assessment Science Integrating point & Non-point Sources) as one of the models that the USEPA supports and recommends for state and federal agencies to use to address point and nonpoint source pollution control.The hydrological processes are divided into two phases: the land phase and the channel/floodplain phase.The SWAT model uses the SCS curve number procedure when daily precipitation data is used while the Green-Ampt infiltration method is chosen when sub-daily data is used to estimate surface runoff.The SCS curve number equation is: where Q surf is the accumulated runoff or rainfall excess (mm H 2 O); R day is the rainfall depth for the day (mm H 2 O); I a is the initial abstractions, which includes surface storage, interception, and infiltration prior to runoff (mm H 2 O); and S is the retention parameter (mm H 2 O).The retention parameter varies spatially due to changes in soil, land use, management, and slope and temporally due to changes in soil water content.The retention parameter is defined as: where CN is the curve number for the day.
The SWAT model uses the Modified Universal Soil Loss Equation (MUSLE) to estimate sediment yield at HRU (Hydrological Response Units) level.The MUSLE is defined as: where Q sed is the sediment yield on a given day (metric tons); Q surf is the surface runoff volume (mm H 2 O ha −1 ); q peak is the peak runoff rate (m 3 s −1 ); A hru is the area of the HRU (ha); K usle is the USLE soil erodibility factor; C usle is the USLE cover and management factor; P usle is the USLE support practice factor; L usle is the USLE topographic factor; and F CFEG is the coarse fragment factor.In order to efficiently and effectively apply the SWAT model, different calibration and uncertainty analysis methods have been developed and applied to improve the prediction reliability and quantify prediction uncertainty of SWAT simulations (Arabi et al., 2007).In this study, a parameter sensitivity analysis was performed prior to calibrating the model.Based on the sensitivity ranking results provided by the Morris Qualitative Screening Method (Morris, 1991), the 20 highest ranked parameters affecting stream flow and sediment yield (shown in Table 1) were selected for the following uncertainty analysis using the GLUE method.For modeling accuracy, parameters were calibrated and validated using the highly efficient Sequential Uncertainty Fitting version-2 (SUFI-2) procedure (Abbaspour et al., 2007).The initial parameter range was recommended from the SWAT manual.This calibration method is an inverse optimization approach that uses the Latin Hypercube Sampling (LHS) procedure along with a global search algorithm to examine the behavior of objective functions.The procedure has been incorporated into the SWAT-CUP software, which can be downloaded for free from the EAWAG website (Abbaspour et al., 2009).For the runoff, the Nash-Sutcliffe coefficients during the calibration period and validation period were 0.94 and 0.78, respectively.For the sediment yield, the Nash-Sutcliffe coefficients in the calibration and validation periods were 0.80 and 0.70, respectively.More details can be found in the study of Shen et al. (2008) and Gong et al. (2011).

GLUE method
The GLUE method (Beven and Freer, 2001) is an uncertainty analysis technique inspired by importance sampling and regional sensitivity analysis (Hornberger and Spear, 1981).In GLUE, parameter uncertainty accounts for all sources of uncertainty; i.e. input uncertainty, structural uncertainty, parameter uncertainty and response uncertainty.Therefore, this method has been widely used in many areas as an effective and general strategy for model calibration and uncertainty estimation associated with complex models.In this study, the GLUE analysis process consists of the following three steps:

Step 1: definition of likelihood function
The likelihood function was used to evaluate SWAT outputs against observed values.In our study, the Nash-Sutcliffe coefficient (E NS ) was picked because it was the most frequently used likelihood measure for GLUE based on the literature (Beven and Freer, 2001;Freer et al., 1996;Arabi et al., 2007).
where Q mea,i and Q sim,i are the measured and simulated values for the i-th pair, Q mea is the mean value of the measured values, and n is the total number of paired values.The E NS value ranges from −∞ to 1, with 1 indicating a perfect fit.

Step 2: sampling parameter sets
Due to the lack of a prior distribution of a parameter, uniform distribution was chosen due to its simplicity (Muleta and Nicklow, 2005;Lenhart et al., 2007;Migliaccio and Chaubey, 2008).The range of each parameter was divided into n overlapping intervals based on equal probability (Table 1) and parameters were identically chosen from spanning the feasible parameter range.The drawback of a typical GLUE approach is its prohibitive computational burden imposed by its random sampling strategy.Therefore in this study, an improved sampling method was introduced by combing Latin Hypercube Sampling (LHS) with GLUE.Compared to random sampling, LHS can reduce sampling times and provide a 10-fold greater computing efficiency (Vachaud and Chen, 2002).Therefore, LHS was used for random parameter sampling to enhance the simulation efficiency of the GLUE simulation.Values then were randomly selected from each interval.
If the initial sampling of the parameter space was not dense enough, the GLUE sampling scheme probably could not ensure a sufficient precision of the statistics inferred from the retained solutions (Bates and Campbell, 2001).Hence, a large number of sampling sets (10 000 times) were conducted.Because the SWAT module and the SWAT-CUP software were in different interfaces, all of the 10 000 simulations were calculated manually.The whole simulation period lasted six months on a Centrino Duo at 2.8 GHz computer.

Uncertainty of outputs
For the purpose of determining the extent to which parameter uncertainty affects model simulation, the degree of uncertainty of outputs was expressed by 95CI, which was derived by ordering the 10 000 outputs and then identifying the 2.5 % and 97.5 % threshold values.The 95CI for both stream flow and sediment period were shown in Fig. 2. It was evident that the 95CI for stream flow and sediment was 1-53 m 3 s −1 and 2000-7 657 800 t, respectively.In addition, sediment simulation presented greater uncertainty than stream flow, which might be due to the fact that sediment was affected and dominated by both stream flow processes as well as other factors, such as land use variability (Shen et al., 2008;Migliaccio and Chaubey, 2008).From Fig. 2, the temporal variation of outputs was presented in which an evidently clear relationship was obtained between the amount of the rainfall and the width of confidence interval.This result highlighted an increased model uncertainty in the high precipitation condition.The variability in the uncertainty of sediment was the same as that for runoff, because runoff affects both factors.This could be explained by the fact that uncertainty was inherent in precipitation due to variability in the time of occurrence, location, intensity, and tempo-spatial distribution (Shen et al., 2008).In a hydrology model such as SWAT, although a rainfall event may affect only a small portion of the basin, the model assumes it affects the entire basin.This may cause a larger runoff event to be observed in simulation although little precipitation was recorded due to the limited local extent of a certain precipitation event.In the Three Gorges Reservoir area, the daily stream flow changes frequently and widely, thus the measured value might not represent the actual value of the daily flow and the discrepancy between the measured mean value and simulated mean value would be high.However, more precise simulated flow would depend on designing accurate rain-gauge networks and the existence of fewer measurement errors (Chang et al., 2007).
From Fig. 2, it is clear that most of the observed values were bracketed by the 95CI, 54 % for stream flow outputs and 95 % for sediment.However, several stream flow observations were found to be above the 97.5 % threshold values (such as March, April, November 2004; March, April, May, June, July, August and October 2005; February, March, April, May and July 2006; March, May, June, July and August 2007).Conversely, only one observation (October 2006) was observed below the 2.5 % threshold of sediment output.Measured value was not entirely in the range of 95CI, indicating that the SWAT model could not fully simulate the flow and sediment processes.However, it was acknowledged that for a parameter, model structure and data input can also cause uncertainty in model simulation (Bates and Campbell, 2001;Yang et al., 2007).Based on the results presented in this study, it was not possible to tell the extent to which the errors in the input and model structure contribute on the total simulation uncertainty.However, as parameter uncertainty was only able to account for a small part of whole uncertainty in hydrological modeling, this study suggests further studies are needed on model structure and input in TGRA.Another concern in hydrologic modeling was the equifinality of model parameters (Beven and Binley, 1992;Wagener and Kollat, 2007).Table 2 showed multiple combinations of parameter values yield the same E NS during hydrologic modeling in TGRA.The so-called equifinality showed there was no unique parameter estimation and hence uncertainty in the estimated parameters in TGRA was obvious.This result agreed well with many other studies (Beven and Binley, 1992).This may be due to the fact that parameters obtained from calibration were affected by several factors such as correlations amongst parameters, sensitivity or insensitivity in parameters, spatial and temporal scales and statistical features of model residuals (Wagener et al., 2003;Wagener and Kollat, 2007).It could be inferred that the identifiability of an optimal parameter obtained from calibration should also be evaluated.For an already gauged catchment, a virtual study can provide a point of reference for the minimum uncertainty associated with a model application.This study highlighted the importance of the monitoring task for several important physical parameters to determine more credible results for watershed management.

Uncertainty of parameters
Figures 3 and 5 illustrate the variation of E NS for the Daing River watershed as a function of variation of each of the 20 parameters considered in this study.By observing the dotty plot from Fig. 3, it was evident that the main sources of streamflow uncertainty were initial SCS CN II value (CN2), available water capacity of the layer (SOL AWC), maximum canopy storage (CANMX), base flow alpha factor for bank storage (ALPHA BNK), saturated hydraulic conductivity (SOL K), and soil evaporation compensation factor (ESCO).Among the above six parameters, SOL AWC and CANMX were the most identifiable parameters for the Daing River watershed.This could be explained by the fact that SOL AWC represented soil moisture characteristics or plant available water.This parameter plays an important role in evaporation, which is associated with runoff (Burba and Verma, 2005).It has also been suggested that the soil water capacity has an inverse relationship with various water balance components (Kannan et al., 2007).Therefore, an increase in the SOL AWC value would result in a decrease in the estimate of base flow, tile drainage, surface runoff, and hence, water yield.As shown in Fig. 3, the optimal range of SOL AWC was between [0, 0.2] and better results could be obtained in this interval.By using calibration methods, optimal parameter ranges could also be obtained without much difficulty for other identifiable parameters (CANMX [0, 30], ALPHA BNK [0.3, 1], SOL K [80,300]).However, presence of multiple peaks in the Nash-Sutcliffe model efficiency for CN2 and ESCO indicated that estimation of these parameters might not be feasible.However, it should be noted that non-identifiability of a parameter does not indicate that the model is not sensitive to these parameters.Generally, CN2 was considered as the primary source of uncertainty when dealing with stream flow simulation (Eckhardt and Arnold, 2001;Lenhart et al., 2007).This study showed that CN2 exhibited non-identifiability in the stream flow simulation.This is similar to the study proposed by Kannan et al. (2007).The potential cause would be that there was an explicit provision in the SWAT model to update the CN2 value for each day of simulation based on available water content in the soil profile.Therefore, a change in the initial CN2 value would not greatly affect water balance components.Estimation of non-identifiable parameters, such as CN2 and ESCO for the Daning River watershed, would be difficult as there may be many combinations of these parameters that would result in a similar model performance.Instead of the process calibration, a decision regarding modeling could deal with these non-identifiable parameters by setting a confidence interval on model output.
Figures 4 and 6 illustrate the cumulative parameter frequency for both stream flow and sediment in the Daing River watershed.As shown in Fig. 4, the parameters were not uniformly or normally distributed, especially SOL AWC, CANMX and ESCO.ESCO represents the influence of capillarity and soil crannies on soil evaporation in each layer.Therefore, a change in the ESCO value affected the entire water balance component.When there were higher ESCO values, the estimated base flow, tile drainage and surface runoff increased.The greater uncertainty of this parameter indicated that the soil evaporation probably played a greater role in the whole evaporation process, possibly due to the high air temperature in rainy seasons in the TGRA.In comparison, other parameters such as CN2 and SOL K, were close to a uniform distribution while they were also more or less skewed.This non-linearity further implies that the uncertainty in model input did not translate directly into uncertainty in the model outputs, but might rather appear significantly dampened or magnified in the output (Sohrabi et al., 2003).This result demonstrates the important opinion that the model output was influenced by the set of parameters rather than by a single parameter (Beven and Binley, 1992).
Similar to the stream flow simulation, even though many of the parameters were sensitive and affected the sediment simulation, only a small number of the sensitive parameters were identifiable.As shown in Fig. 5, the factors of uncertainty for sediment were CN2, Manning's value for main channel (CH N2), maximum canopy storage (CANMX), base flow alpha factor for bank storage (ALPHA BNK), exp.Re-entrainment parameter for channel sediment routing (SPEXP), lin.re-entrainment parameter for channel sediment routing (SPCON), channel cover factor (CH COV) and channel erodibility factor (CH EROD).Clearly, the parameter samples were very dense around the maximum limit (Fig. 6).Summarizing the information in Figs.3-6, it can be said that the parameters with greater uncertainty of stream flow mainly come from the surface corresponding process and the parameters with greater uncertainty of sediment focused on the channel response process.The results matched well with those of Yang et al. (2011) and Shen et al. (2010).

Conclusions
In this study, the GLUE method was employed to assess the parameter uncertainty in the SWAT model applied in the Daning River Watershed of the Three Gorges Reservoir Region (TGRA), China.The results indicate that only a few parameters were sensitive and had a great impact on the stream Hydrol.Earth Syst.Sci., 16, 121-132, 2012 www.hydrol-earth-syst-sci.net/16/121/2012/ flow and sediment simulation.CANMX, ALPHA BNK and SOL K were identified as identifiable parameters.The values of these parameters could be obtained by using the calibration process without much difficulties.Conversely, there were multiple possible values for CN2 and ESCO.This indicates that calibration of these parameters might be infeasible.These non-identifiabiable parameters even led to equifinality in hydrologic and NPS modeling in the TGRA.It was anticipated that the parameter uncertainties are systematically correlated to the non-identifiable parameters.Under such cases, a user should check if any information related to the watershed characteristics and its underlying hydrologic processes could be used to provide a more precise range for the model parameter.It was anticipated that this study would provide some useful information for hydrological modeling related to policy development in the Three Gorges Reservoir Region (TGRA) and other similar areas.
It is suggested that more detailed measurement data and more precipitation stations should be established in the future for hydrological modeling in the TGRA.In addition, further studies should be continued in the field of model structure and input to quantify hydrological model uncertainty in the TGRA.

Table 1 .
The range and optimal value of model parameter.

Table 2 .
The equifinality of model parameters.