On the assessment of the impact of reducing parameters and identification of parameter uncertainties for a hydrologic model with applications to ungauged basins

Journal of Hydrology 320 (2006) 37–61 www.elsevier.com/locate/jhydrol

On the assessment of the impact of reducing parameters and identification of parameter uncertainties for a hydrologic model with applications to ungauged basins Maoyi Huang, Xu Liang* Department of Civil and Environmental Engineering, University of California, 631 Davis Hall # 1710, Berkeley, CA 94720-1710, USA

Abstract In this paper, we investigate model parameter uncertainties associated with hydrological process parameterizations and their impacts on model simulations in the Three-Layer Variable Infiltration Capacity (VIC-3L) land surface model. We introduce an alternative subsurface flow parameterization into VIC-3L to reduce the impacts of model parameter uncertainties on model simulations by reducing the number of model parameters that need to be estimated through a calibration process. The new subsurface flow parameterization is based on the concepts of kinematic wave and hydrologic similarity, and has one parameter for calibration. Results from the 12 MOPEX (Model Parameter Estimation Experiment) basins obtained by applying the VIC3L model with the new subsurface flow formulation show that the performance of the new parameterization is comparable to the original subsurface flow formulation, which has three parameters for calibration. In addition, a probabilistic approach based on Monte Carlo simulations is used to evaluate model performance and uncertainties associated with model parameters over different ranges of streamflow. Studies based on the 12 MOPEX watersheds show that compared to the parameter associated with the new subsurface flow parameterization, the VIC shape parameter (i.e. the b parameter that represents the shape of the heterogeneity distribution of effective soil moisture capacity over a study area) has a larger impact on model simulations and could introduce more uncertainty if not estimated appropriately. Furthermore, investigations on the b parameter suggest that the ensembles (i.e. the mean response and its bounds) from the Monte Carlo simulations could provide reasonable predictions and uncertainty estimates of streamflows, which have important implications for applications to ungauged basins. The study also shows that appropriate reduction of the number of model parameters is an effective approach to reduce the impacts of parameter uncertainties on model simulations. This is more so for applications to ungauged basins or basins with limited data available for calibration. The new subsurface flow parameterization and the probabilistic uncertainty analysis approach are general and can be applied to other modeling studies. q 2005 Elsevier Ltd. All rights reserved. Keywords: MOPEX; Monte Carlo; Simulation; Uncertainty; Subsurface flow; Model parameters; Land surface model; Ungauged watersheds

1. Introduction * Corresponding author. Tel.: C1 510 642 2648; fax: C1 510 642 7483. E-mail address: [email protected] (X. Liang).

0022-1694/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jhydrol.2005.07.010

Uncertainties associated with hydrological modeling generally come from two main sources. They are: (1) the uncertainties associated with available data

38

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

(e.g. forcing data, model input data, initial and boundary conditions, etc.); and (2) uncertainties arising from model structure and model parameters. The first type of uncertainty is and will be continually reduced with advancements in information and measuring technology. The second type is related mainly to our understanding of hydrological processes and the way these processes are parameterized. Therefore, reductions of the second type of uncertainty rely on improved understanding of the physics and effective representation of the associated hydrological processes. This can be achieved through improvements in model structure and model parameter estimations. Tremendous efforts have been carried out to reduce uncertainties associated with the second type by finding optimal model parameter sets using model calibration techniques and historical data for a given watershed. Efforts include developments of measures used to evaluate goodness-of-fit between model simulations and observations (e.g. Sorooshian and Dracup, 1980; Sorooshian, 1981), and of global optimization schemes, such as the shuffled complex evolution (SCE-UA) method developed by the University of Arizona (e.g. Duan et al., 1992, 1993; Sorooshian et al., 1993), used to obtain ‘optimal’ parameter sets. These optimization schemes generally assume that an optimal parameter set exists and implicitly ignore the estimation of predictive uncertainties. Also, recent studies by researchers (e.g. Klepper et al., 1991; van Straten and Keesman, 1991; Beven and Binley, 1992; Yapo et al., 1996) suggest that a single optimal parameter set for a hydrologic model may not exist and the uncertainties associated with the optimal parameter sets could be large. A model with the optimal parameter set may have the best fit over the period of the calibration data, but there may exist multiple parameter sets that are as good as the ‘optimal’ set. In addition, using different performance evaluation criteria could result in different optimal parameter sets. These resulted parameter sets are referred to as, for example, ‘equifinality’ (Beven and Binley, 1992), ‘equally probable parameter sets’ (van Straten and Keesman, 1991), and ‘acceptable sets’ (Klepper et al., 1991). One approach to address these limitations is to provide predictions within a range to reflect the fact that an optimal parameter set alone is not enough to

represent the possible uncertainty associated with the model predictions. Thus, the parameter space needs to be sampled to generate realizations of the model simulations so that the prediction range can be estimated based on the model simulations. Several studies use this approach, for example, the generalized likelihood uncertainty estimation (GLUE) method (Beven and Binley, 1992; Freer et al., 1996), and the prediction uncertainty method (PU) (e.g. Klepper et al., 1991). However, these methods have not been widely applied to hydrological modeling due to two reasons: (1) an effective and efficient method to sample the parameter space is needed, and (2) a subjective choice of the model performance evaluation criteria needs to be specified (e.g. in the GLUE method). Kuczera and Parent (1998) pointed out that the above algorithms were sensitive to the space to be sampled. Therefore, extensive sampling in the parameter space is usually required by these algorithms. Comparing the GLUE algorithm to a Markov Chain Monte Carlo (MCMC) (Metropolis et al., 1953; Hastings, 1970) algorithm, Kuczera and Parent (1998) concluded that importance sampling based on the GLUE algorithm is inferior to the MCMC sampling, because the latter is capable of producing reliable inferences of the posterior probability distribution of the parameters with modest sampling. The advantage of the MCMC sampling becomes more significant when the sampling space has a high dimension. Inspired by the promising features of the MCMC sampling method, Vrugt et al. (2003a) recently proposed a Shuffled Complex Evolution Metropolis (SCEM-UA) algorithm that combines the strengths of the classical Markov Chain Monte Carlo (MCMC) methods with those of complex shuffling method (Duan et al., 1992). The SCEM-UA algorithm is an adaptive MCMC sampler that can sample the parameter space effectively and efficiently in the framework of Bayesian inference and is able to estimate uncertainty of parameters based on the inferred posterior distribution. Although the MCMC sampling is better than the basic importance sampling, we investigate the parameter uncertainties in this study by employing a probabilistic approach based on importance sampling, owing to the lowdimensional sampling spaces of only one or two parameters. The MCMC sampler can be employed in

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

future studies for cases with high-dimensional sampling spaces. The other approach argues that significant improvement in the optimized parameter sets can be achieved by using additional information for calibration and/or validation (e.g. De Grosbois et al., 1988; Yan and Haan, 1991a,b) or by exploiting complementary information in the data within a multicriteria framework (Yapo et al., 1998; Gupta et al., 1998; Boyle et al., 2000). The multi-criteria framework is based on the concept of Pareto Optimal Sets. An efficient and effective global optimization algorithm, called multi-objective complex evolution (MOCOM-UA) (e.g. Yapo et al., 1998) is developed to identify the Pareto sets based on the SCE-UA method (Duan et al., 1992, 1993). Gupta et al. (1999) show that the MOCOM-UA method is able to achieve a calibration comparable to the manual calibration by an experienced hydrologist in that it considers multiple requirements in fitting the data. However, it suffers from several weaknesses (Gupta et al., 2003), including the tendency of clustering the Pareto solutions into a central compromise region of the Pareto set, and a tendency of converging prematurely for case studies that involve large numbers of parameters and strongly correlated performance criteria. In recognizing of these weaknesses, Vrugt et al. (2003b) proposed the Multi-objective Shuffled Complex Evolution Metropolis (MOSCEM-UA) algorithm that employs a newly developed and improved concept of Pareto dominance to approximate a fairly uniform Pareto frontier. As pointed out by Gupta et al. (1998), the first and second (i.e. multicriteria) approaches are complementary to each other and need to be used in conjunction with one another in real applications. In this paper, we seek to reduce and quantify uncertainties associated with the hydrologically based Three-Layer Variable Infiltration Capacity (VIC-3L) land surface model simulations from two aspects. The first is from the underlying model structure by modifying the parameterization of the subsurface flow so that the model has fewer parameters for model calibration, but retaining a sound physical basis. The purpose of modifying the subsurface flow parameterization is to reduce potential uncertainties related to parameter estimations, especially when the model is applied to ungauged basins. The second is to quantify

39

the VIC-3L model sensitivity and uncertainty when used with the modified subsurface flow parameterization. By using a probabilistic framework, we analyze the strengths and weaknesses of the modified VIC-3L model in capturing different hydrological events (e.g. floods and droughts) in a statistical sense. Although the subsurface flow parameterization and the uncertainty analysis method are applied to the VIC-3L model, the concepts and methods are general and can be applied to other models and parameter uncertainty analyses. The remainder of the paper is arranged as follows. In Section 2, a brief introduction of the Three-Layer Variable Infiltration Capacity (VIC-3L) land surface model is provided. Problems associated with the ARNO subsurface flow parameterization in the current version of the VIC-3L model are discussed, and a new subsurface flow parameterization is proposed. In Section 3, a probabilistic approach to evaluate and quantify the sensitivities and uncertainties of the VIC-3L model is presented, and the performance of the VIC-3L model is analyzed and discussed. Conclusions are summarized in Section 4.

2. Subsurface flow formulation for the VIC-3L model 2.1. Overview of the VIC-3L model The Three-Layer Variable Infiltration Capacity (VIC-3L) model (Liang et al., 1994, 1996a,b, 1999, 2003; Cherkauer and Lettenmaier, 1999; Liang and Xie, 2001) is a hydrologically based land surface scheme. The spatial variabilities of infiltration, precipitation, and vegetation are considered partially in the model to simulate water and energy budgets at the land surface. The dynamics of runoff generation includes two components, fast runoff and slow runoff. Direct runoff, including both saturation and infiltration excess runoff (Liang et al., 1994; Liang and Xie, 2001), represents the fast component, while the ARNO subsurface flow (Francini and Pacciani, 1991; Todini, 1996) represents the slow component. The upper soil layer of the model is designed to represent the dynamic response of the soil to rainfall events, and the lower layer is used to characterize the seasonal soil moisture behavior through the slower

40

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

dynamics of inter-storms associated with deeper soil moisture and subsurface flow processes. In this application, we do not explicitly represent groundwater table. In other words, the unsaturated and saturated zones are treated in a lumped sense. However, VIC-3L can explicitly deal with the dynamics of surface and groundwater interactions and calculate groundwater table positions (Liang et al., 2003). Besides the upper and lower layers, VIC3L also has a thin top soil layer, which allows for quick bare soil evaporation following small summer rainfall events. VIC-3L includes soil moisture diffusion processes between the three soil layers (Liang et al., 1996a,b). VIC-3L can be applied under cold climate where snow and frozen soil processes play important roles (Cherkauer and Lettenmaier, 1999). The VIC-3L model has been used in a wide range of applications from soil moisture estimation (e.g. Nijssen et al., 2001a), streamflow forecasting (e.g. Nijssen et al., 2001b; Hamlet and Lettenmaier, 1999a,b), to climate change impact analyses (e.g. Leung et al., 1999), and land data assimilation systems (e.g. Mitchell et al., 2003). One of the main difficulties in applying land surface models or hydrological models is to estimate model parameters appropriately. In VIC-3L model applications, for example, there are seven model parameters (i.e. b, Ds, Dsmax, Ws, two soil depths, and Kroute (the meanings of these parameters are listed in Table 2) that need to be calibrated because of limited available information to quantify them directly. The parameter b is called VIC shape parameter, representing the shape of the heterogeneity distribution of the effective soil moisture capacity over a study area. The distribution of the effective soil moisture capacity within an area can be expressed as w Z wm ½1Kð1KAÞ1=b 

(1)

where w and wm are the effective point and maximum point soil moisture capacity, respectively; A is the fraction of an area for which the effective soil moisture capacity is less than or equal to w. Ds, Dsmax, and Ws are the parameters associated with the ARNO subsurface flow formulation that is currently implemented in VIC-3L. The name of the ARNO formulation originated from its first application to the Arno River in Italy (Todini, 1996). Dsmax

Fig. 1. A schematic representation of the ARNO subsurface flow parameterization.

represents the maximum subsurface flow, Ds and Ws represent, respectively, the fractions of Dsmax and the third layer soil moisture W3,max (see Fig. 1). Readers are referred to Liang et al. (1994) for details. When the three parameters that are related to the ARNO formulation are calibrated with sufficient and representative streamflow data for study watersheds, the VIC-3L model simulated streamflows can reproduce the observations quite well (e.g. Nijssen et al., 2001a,b; Todini, 1996). However, as pointed out by Todini and his co-workers (Todini, 1995; Todini and Dumenill, 1999), the major disadvantage of the ARNO parameterization is its lack of physical grounds, so that these parameters need to be determined through calibration using long-term streamflow records. Thus, such a weakness of the ARNO subsurface flow formulation hinders its applications to ungauged basins. Warrach et al. (2002) also reported similar limitations in their study, in which they compared results obtained from the Surface Energy and Water Balance (SEWAB) land surface model (Mengelkamp et al., 1999; Warrach et al., 2001) by using the VIC runoff parameterizations (i.e. the VIC saturation excess runoff and ARNO subsurface flow formulations) in SEWAB versus the results using the Topographybased Hydrological Model (TOPMODEL) (e.g. Beven et al., 1984, 1995; Beven, 1997) runoff parameterizations in SEWAB. The authors showed that both versions performed well and were computationally efficient, but there were more free

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

parameters in the ARNO parameterization that needed to be calibrated. In this paper, we present an alternative subsurface flow parameterization that has a sound physical basis but has fewer parameters than ARNO. It is worth mentioning that our objective here is to find an alternative subsurface flow formulation to facilitate applications of the VIC-3L model to ungauged basins without sacrificing model performance in general. 2.2. New subsurface flow parameterization Tallaksen (1995) reviewed different ways of characterizing subsurface flow recession rates and pointed out that the recession curves from many catchments could be represented by a general relationship as follows Q Z KSp

(2)

where S is a storage term and K and p are constants. The form of Eq. (2) for describing the subsurface flow is consistent with the findings reported by Zecharias and Brutsaert (1988), Troch et al. (1993), Brutsaert (1994), Brutsaert and Lopez (1998), and Eng and Brutsaert (1999), in which subsurface flows at hillslope and watershed scales are studied from both theoretical and empirical perspectives. In this section, we will show that our new parameterization follows the general form of Eq. (2). Our new subsurface flow parameterization is based on the concepts of kinematic wave and hydrologic similarity, and is directly related to soil properties (e.g. the hydraulic conductivity), topography (e.g. topographic index), and soil moisture of each model grid cell (or a study area). The kinematic wave approach has been widely used in hydrology and water resources research (e.g. Singh, 2001). For example, the subsurface flow formulation (see Eqs. (3) and (4)) of the topography-based hydrological model (TOPMODEL) has been developed based on it, and shown to be a good approximation of the subsurface flow systems based on field and numerical studies (Beven, 1981) and in a number of applications (e.g. Beven et al., 1984, 1995). The transmissivity profile in TOPMODEL is expressed by T Z T0 expðKZ=mÞ

(3)

41

where T0 is the saturated lateral transmissivity, Z is the local storage deficit or depth to water table; m is a scaling parameter describing the transmissivity profile with depth. The subsurface flow [L/T] of TOPMODEL for a watershed with an area of A is expressed by
Qb Z T0 exp½KðZ=mÞKL 1

(4)

where Z is the mean storage deficit or mean depth to water table over the watershed; L1 is the spatial average of logarithmic topographic Ð indices over the watershed defined as L1 Z ð1=AÞ A lnða=tan bÞds; a is the contributing area to a point per unit contour length in the watershed and b is the slope of an aquifer, which is assumed to be the same as the surface slope, and m is assumed to be homogeneous over the watershed. The choice of an exponential transmissivity profile in TOPMODEL was based on observations in upper soil layers (Beven, 1984) as well as its attractive analytical expression. However, such an assumption could have a significant impact on the shape of simulated subsurface flow recession curves. Observing that the simulated subsurface flow obtained by Eq. (4) was not appropriate for the Ringelbach watershed, Ambroise et al. (1996) generalized Eq. (3) to incorporate two alternative forms of the transmissivity profiles, which are linear and parabolic, respectively. Duan and Miller (1997) and Iorgulescu and Musy (1997) further pointed out that a more generalized power function could be used to represent the transmissivity profiles, which is expressed as T Z T0 ð1KZ=MÞn

(5)

where M is the maximum storage deficit, and n is a non-dimensional parameter that ranges from one to infinity. When n goes to infinity, Eq. (5) approaches the original exponential profile expressed by Eq. (3) with mZM/n. Incorporating the transmissivity profile of Eq. (5) into the framework of TOPMODEL, the subsurface flow [L/T] from a watershed then becomes !
n 1K MZ Qb Z T0 (6) L2 Ð where L2 Z ð1=AÞ A ða=tan bÞð1=nÞ ds is a spatial average of the soil-topographic index to the power of 1/n

42

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

over a watershed. It can be seen that Eq. (6) falls into
the form of Eq. (2), with ð1KðZ=MÞÞ to be the storage term, T0 =L1=n and n to be the two different constants. 2 There are two major difficulties in the implementation of Eq. (4) or (6) when applied to a grid-based land surface model. The first is that the contributing area a per hillslope length [L] at a point, in Eqs. (4) and (6), is defined as the area that includes all the upslope points draining to this point on a hillslope, which should be confined by natural watershed or subwatershed boundaries. Thus, the basic computational units of TOPMODEL is represented by sub-watersheds so that the contributing area at any point of the modeling domain can be calculated by taking into account all of the upslope points on the hillslope draining to this point. For details about the calculation of a, readers are referred to Quinn et al. (1995). For a grid-based land surface model, the modeling domains are based on grid cells (e.g. 25 km by 25 km grids), whose boundaries could truncate a hillslope into different grids. Stieglitz et al. (1997) incorporated the TOPMODEL formulations including Eq. (4) into a land surface model and suggested the use of a subwatershed as a basic computational unit in a land surface model. From a hydrological point of view, the use of sub-watersheds to calculate the runoff is more meaningful and accurate. From a climate modeling point of view, however, this would introduce difficulties to link climate and land surface models at their interfaces. The second difficulty, as shown by Eqs. (4) and (6), is that information of the mean depth to water table
over a model grid cell is required. For the new ðZÞ version of VIC-3L, this requirement is not a problem because Liang et al. (2003) introduced a new framework to represent the dynamics of surface and groundwater interactions into the VIC-3L model in which the groundwater table position is calculated explicitly. However, most of the land surface models or hydrologic models do not represent groundwater table explicitly at present, rather the groundwater table is lumped within a soil layer. To facilitate applications to these models, a modified represen
tation of the storage term ð1KðZ=MÞÞ in Eq. (6) is necessary so that the subsurface flow can be calculated based on lumped soil moisture content simulated by a land surface model. In this section, a new subsurface flow parameterization that partially

overcomes the two limitations is presented based on the concepts of TOPMODEL (e.g. Beven et al., 1984; Beven et al., 1995; Beven, 1997) and the Topographic Kinematic Approximation and Integration (TOPKAPI) model by Todini and his co-workers (Todini, 1995; Ciarapica and Todini, 2002; Liu and Todini, 2002). Assuming that a model grid cell is overlaid by a finer resolution DEM, conceptually, we need to sum all of the subsurface flows generated by each of the DEM pixels within the model grid cell to obtain the total amount of subsurface flow for the model grid cell. Since the subsurface flow generated over some DEM pixels within the model grid cell will eventually drain to the DEM pixels at the foot of each hillslope within the same model grid cell, the subsurface flow from a particular model grid cell then becomes the sum of all the subsurface flows leaving the DEM pixels at the foot of the hillslopes that are contained in the model grid cell. By calculating the subsurface flow based on the model grid cell rather than the boundary of a sub-watershed or a watershed, we can easily couple a land surface model with a climate model. Compared to the scale of a model grid cell (e.g. on the order of w30 km or larger in length), the scale of hillslopes (e.g. on the order of w100 m in length) is quite small so that even if a hillslope is truncated by the boundary of a model grid cell, the differences between the calculated contributing areas and the true contributing areas from the DEM pixels for each hillslope within the model grid cell could be small. In addition, the subsurface flow contributed by the areas truncated for the model grid in question will be included in other model grids. Thus, the total amount of subsurface flow for the entire study area will remain to be appropriately represented. The misrepresentation of the subsurface runoff timing for each individual model grid due to the truncation can be partially taken care of by the river network routing. Furthermore, this kind of truncated representation is consistent with the model grid-based framework for surface runoff calculation in which model grids rather than sub-watershed or watershed boundaries are used, where it is assumed that the topography controls the subsurface and surface runoff in a similar way. Therefore, we propose to parameterize the subsurface flow based on the boundaries of each model grid instead of each sub-watershed by integrating

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

a subsurface flow formulation from the hillslope scale to the model grid scale. According to Benning (1994), the transmissivity at a point can be approximated by using an average relative soil moisture over a soil layer that corresponds to the point instead of using the soil depth to ground water table as used in TOPMODEL, and it can be expressed as follows without introducing significant errors Ti ðq~i Þ Z

D ði

kðq~i ðzÞÞdz zksx;i $Di $Qci

(7)

0

0

unit width of contour length [L2/T] from the soil column at point i can be calculated based on the kinematic wave formulation as follows qi Z sin bi $ksx;i $Di $Qci i

(8)

However, in land surface models, Qi at each point is unknown. What is known is the relative soil moisture content averaged over an entire Ð model grid
ð1=A 0 Þ A 0 Qi ds). To cell with an area of A 0 (i.e. QZ relate the vertically averaged relative soil moisture
we apply the hydrologic content Qi at a point to Q, similarity concept proposed by Kirkby (1975), which was later refined by Beven and Kirkby (1979). Specifically, following the assumptions of steadystate and homogeneous recharge made in TOPMODEL, the subsurface flow from a soil layer per unit contour length [L2/T] can be expressed as qi Z ai R

contour length [L]; R is a spatially homogeneous recharge rate [L/T] over a. By combining Eqs. (8) and (9) and the hydrologic similarity concept, and assuming that ksx,i, Di, and ci are uniformly distributed within the model grid cell as made by others (e.g. Ambroise et al. 1996; Iorgulescu and
and c
Musy, 1997), and are represented by k
sx , D, (i.e. assuming effective soil properties for each model grid), we can obtain the following relationship Qi Z

where i is the index of points (e.g. a DEM pixel), Di is the depth of the soil column at point i in which groundwater table is not explicitly expressed; z represents the depth in vertical direction. q~i ðzÞZ ðqi ðzÞKqr;i ðzÞÞ=ðqs;i ðzÞKqr;i ðzÞÞ represents the relative soil moisture content at depth z, where qi(z), qr,i(z), and qs,i(z) are the soil moisture content, residual soil moisture content, and saturated soil moisture content at that depth, respectively; ksx,i is the saturated lateral hydraulic conductivity, ciZ2BiC3, where Bi is the Clapp and Hornberger (1978) pore-size distribution index; and Qi represents the average relative soil moisture content which is expressed as D Ði Qi Z ð1=Di Þ q~i ðzÞdz. Therefore, the lateral flux per

(9)

where ai is, as defined earlier, the contributing area to a DEM pixel within the model grid cell per unit

43

li1=c

Q L

(10)

where liZai/sin bi is a topographic index at a DEM pixel within a model grid cell. Because the slope bi is small, sin bi can be approximated by tan bi, i.e. Ð lizai/tan bi. LZ ð1=A 0 Þ A 0 ðli Þ1=c
ds is a spatial average of the topographic index over a model grid cell with an area of A 0 , instead of over a watershed area or a sub-watershed area A as in L2. The subsurface flow for a particular model grid cell, in a unit of [L3/T], can be then obtained by summing all the subsurface flows over all the pixels at the foot of hillslopes within the model grid cell based on Eqs. (8) and (10) X Qdischarge Z qk !l b k


c
X Q

ak !l Z ksx $D$ L k

(11)

where k is the index of the DEM pixels identified to be at the foot of the hillslopes within the model grid cell, l is the length of a DEM pixel. Based on the fact that subsurface flows from all other pixels in this model grid cell will eventually drain to the pixels at the foot of the hillslopes, it is intuitive that P 0 k ak lZ A . Therefore, Eq. (11) expressed in terms of subsurface flow Qb per unit area per model grid cell (i.e. with a unit of [L/T] can be rewritten as 
c
Q

Qb Z aks $D L

(12)

where a is a parameter introduced by Chen and Kumar (2001) to account for anisotropic properties of saturated hydraulic conductivities between the

44

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

lateral direction ðk
sx Þ and the vertical direction ðk
s Þ, and k
sx Z ak
s . It is important to note that the subsurface flow, Qb, from each model grid cell will be routed through the river network in the same way as the overland flow generated from each model grid cell. Also, it is assumed that there are no lateral flow interactions between model grid cells for the subsurface flow, Qb, as in the case for overland flow in most of the gridbased land surface models. It is worth pointing out that Eq. (12) is similar to Eq. (6), since both of them take the form of Eq. (2). Eq. (12) differs from Eq. (6) in that its storage term is represented in terms of spatially averaged relative soil moisture content instead of the depth to water table; and that its saturated transmissivity T0 is estimated by
By comparing Eqs. (5) and (7), it can T0 Z a$k
s $D. be seen that both of the transmissivity profiles follow the power law relationship. In Eq. (5), the power law coefficient (i.e. n) is a parameter that needs to be calibrated or estimated based on recession curve analyses, while in Eq. (7), the power law coefficient
is characterized by a known soil property (i.e. c) parameter so that no calibration or estimation is needed. Furthermore, Eq. (10) shows that the relative soil moisture content at each DEM pixel in the model grid cell can be related to the spatially averaged relative soil moisture content of the model grid cell through the topographic index, which is analogous to the relationship between the local and mean water table depths in TOPMODEL. Eq. (12) is different from the formulation used in the TOPKAPI model in which Eq. (8) is applied to estimate subsurface flow in a distributed model framework at spatial resolutions ranging from 200 m to 1 km. For applications at larger spatial scales that are commonly used by a land surface model (e.g. on the order of tens of kilometers or larger), Eq. (12) would be more appropriate. The subsurface flow parameterization in Eq. (12) is
and L have
D, physically based, in which k
s , c; physical meanings and can be estimated directly based on soil and topographic characteristics. Q
is a state variable simulated at every time step. With Eq. (12), the number of parameters that need to be calibrated is reduced from three in the ARNO formulation (i.e. Ds, Dsmax, and Ws) to one parameter (i.e. a) in the new parameterization, and the limitations of the ARNO subsurface flow formulation

that are discussed in Section 2.1 are partially overcome. 2.3. Implementation of the new subsurface flow parameterization in the VIC-3L model Eq. (12) is applied, respectively, to layers 2 and 3 in the VIC-3L model. Specifically, the subsurface flow from each of the layers is calculated as 
c
Qj

; (13) Qb;j Z aks;j $Dj L where jZ2,3 denoting layers 2 and 3, respectively; D
j , k
s;j , and Q
j are the soil depth, saturated vertical hydraulic conductivity, and relative soil moisture content of layer j, respectively. Here, the mean topographic index L is derived based on the surface topography and is taken to be the same for layers 2 and 3. The total subsurface/baseflow then becomes Qb Z Qb;2 C Qb;3

(14)

Note that with such a parameterization, the interflow part of the subsurface flow is added into the VIC-3L model through Qb,2, which is absent in the original VIC-3L model framework, although its value is generally quite small. This is a good addition to the VIC-3L model because interflow could contribute significantly sometimes to the total runoff in some watersheds under some conditions (e.g. Beven, 1989). The original ARNO subsurface flow formulation is now replaced by Qb,3.

3. Case studies: evaluation of subsurface flow parameterizations and impacts of parameter uncertainties on model simulations 3.1. Model parameter estimation experiment (MOPEX) data sets Data used in this study are provided by the second MOPEX Workshop held in Tucson, AZ, during April 8–10, 2002. The primary scientific goal of MOPEX is to develop techniques for a priori estimation of model parameters used in land surface parameterization schemes (Schaake et al., 2001). To achieve this goal, MOPEX assembled historical hydrometeorological

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

45

Table 1 Information about hydrological and climatic characteristics of the 12 watersheds provided by MOPEX Tucson workshop USGS Stream gage ID

Basin area (km2)

Gage longitude

Gage latitude

01608500

3820

K78.65

39.45

01643000

2116

K77.38

01668000

4134

03179000 03364000 03451500 03054500 05455500 07186000 07378500 08167500 08172000

1020 4421 2448 2372 1484 3015 3315 3406 2170

PPT/PET

Basin name

562

1.64

39.39

232

1.15

K77.52

38.32

55

1.2

K81.01 K85.93 K82.58 K80.04 K91.72 K94.57 K90.99 K98.38 K97.65

37.54 39.20 35.61 39.15 41.47 37.25 30.46 29.86 29.67

1527 603 1950 1281 633 833 0 948 322

South branch potomac river near Springfield, WV Monocacy river at jug bridge near Frederick, MD Rappahannock river near Fredericksburg, VA Bluestone river near Pipestem, WV East fork white river at Columbus, IN French broad river at Asheville, NC Tygart Valley river at Phillipi, WV English river at Kalona, IA Spring river near Waco, MO Amite river near Denham Springs, IA Guadalupe river near Spring Branch, TX San Marcos river at Luling, TX

data and river basin characteristics for about 200 river basins (500–10,000 km2) over a wide range of climate conditions throughout the world. The MOPEX Tucson data sets include 12 basins from the United States. These basins are screened by the MOPEX organizers so that each basin is covered by an adequate number of rain gauges. Climate conditions of these 12 watersheds vary from arid to wet conditions. Brief information about the 12 watersheds is listed in Table 1, while detailed information can be found from Schaake et al. (2001). Hourly time series of meteorological forcing data for a period from 01/ 01/1960 to 12/31/1998, vegetation and soil data, and daily streamflow are provided by MOPEX. 3.2. Application of the modified subsurface flow parameterization In this section, the version of the VIC-3L model with the new subsurface flow parameterization (i.e. Qb3 or Eq. (12)) is tested and compared to the VIC-3L model simulations with the ARNO subsurface flow parameterization, in which each watershed is treated as one single grid cell due to the limited forcing information provided by MOPEX. Then the model parameter uncertainties are investigated through Monte Carlo simulations for the feasible range of the parameters. The topographic indices needed by Eq. (12) for the 12 watersheds are abstracted from HYDRO1k

Elev. (ft)

1.5 1.21 2.34 1.76 0.89 0.96 1.46 0.45 0.56

(Verdin and Greenlee, 1996), a geographic database developed by the United States Geological Survey (USGS). HYDRO1k provides geo-referenced data sets of aspects, flow directions, flow accumulations, slopes, topographic indices, drainage basin boundaries, and streamlines covering six continents of the world at a spatial resolution of 1 km. Because the VIC-3L model is applied to each watershed as one VIC-3L model grid cell in this study, the contributing area (i.e. a) calculated based on each watershed is applicable so that the topographic index from HYDRO1k can be used. To test the new subsurface flow parameterization, both the original and modified versions of the VIC-3L model are applied to the 12 watersheds by first calibrating the model with the forcing data in the calibration period (i.e. 1979–1983) and then using the calibrated parameters to simulate the streamflow for the entire simulation period (i.e. 1960–1998). Seven parameters (i.e. b, Ds, Dsmax, Ws, two soil depths, and Kroute, as discussed in Section 2.1) in the original VIC3L model are calibrated, while five parameters (i.e. b, a, two soil depths, and Kroute) in the modified version are calibrated. The ranges of the model parameters used in the calibrations and their meanings are summarized in Table 2. The range of the VIC b parameter is set to be [0,5] according to Huang et al. (2003), while the range of a is set to be [0.1, 1000]. Kumar (2004) provided

46

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

Table 2 Six VIC-3L model parameters and one routing parameter that are normally calibrated and their corresponding physical meanings Parameter

Unit

Range

Physical meaning

B Ws

N/A Fraction

0–5 0–1

Dsmax Ds d1 D2 Kroute

mm/day Fraction M M N/A

0–40 0–1 0.1–1 0.5–4 0–10

Soil moisture capacity curve shape parameter Fraction of the maximum 3rd layer soil moisture, which is the threshold of linear or non-linear relationship between the subsurface flow and the 3rd layer soil moisture The maximum baseflow rate of the basin Fraction of Dsmax corresponding to Ws Thickness of the 2nd soil layer Thickness of the 3rd soil layer Routing water storage parameter

a list of values for a to be between 0.003 and 2000. However, Niu and Yang (2003) reported, based on their sensitivity study, that the value of a cannot exceed e5, otherwise the soil would dry out. Therefore, we use [0.1, 1000] for a in this study. A shorter calibration period is used in this study compared to the requirement of the MOPEX. In the MOPEX model inter-comparison experiments, the models are required to be calibrated over a period of 01/01/1979–12/31/1998. However, only five years (01/01/1979–12/31/1983) were used in our study to calibrate VIC-3L parameters (seven with the original version and five with the new version) due to limited computer resources and the limited time available to submit the model results. The VIC-3L model parameters are calibrated using the SCE-UA method (e.g. Duan et al., 1992, 1993) at hourly time step for 10 out of the 12 watersheds. The remaining two watersheds (i.e. the ones associated with USGS stream gauge IDs 08167500 and 08172000) are under relatively arid climate with annual precipitation of 762 and 860 mm, respectively. For the period of 01/ 01/1979–12/31/1983 there were only a few storm events for the two arid watersheds, which do not provide enough information for calibration. Therefore, a longer period (i.e. 01/01/1979–12/31/1998) is used for calibration for the two arid watersheds. The model is warmed up with a period of 1 year by repeating the data of year 1979 during the calibrations. Considering the fact that the calibration period is only a small fraction of the entire simulation period, we decide to present in this paper the results for the calibration period and simulation period (i.e. 01/01/ 1961–12/31/1998, including the calibration period), with the results from the first year discarded as

the warming-up year (i.e. 01/01/1960–12/31/1960). The results for the validation period specifically (i.e. 01/01/1961–12/31/1978) are reported elsewhere and compared with those of other models participating in the MOPEX experiment (Duan et al., 2005). Furthermore, to evaluate the ability of both the original and modified versions of the VIC-3L model to simulate streamflow in ungauged basins, a priori experiments are also carried out for the simulation period (i.e. 01/01/1961–12/31/1998). For the original version, the seven parameter (i.e. b, Ds, Dsmax, Ws, two soil depths, and Kroute) are set to be at their default values, which are 0.10, 0.01, 10.00, 0.75, 0.30, 1.00, and 3.0, respectively, for all the 12 basins. For the modified version, the values of the five parameters (i.e. b, a, two soil depths, and Kroute) are set to be 0.10, 500.00, 0.30, 1.00, and 3.00, respectively, where 500.00 is selected for a because it is the median of the range of a for calibration. Note that with this design, the differences between simulated streamflows of the two a priori experiments for each basin would only be introduced by the difference between the subsurface flow parameterizations of the two versions. It is worth mentioning that the a priori results are presented here to facilitate the comparison between the two subsurface flow parameterizations and are likely to be different from those reported in Duan et al. (2005). The Nash–Sutcliffe coefficient is used as the objective function in the SCE-UA method. The Nash–Sutcliffe coefficient is expressed as PN ðqs ðiÞKqo ðiÞÞ2 E Z 1K PiZ1 N
o Þ2 iZ1 ðqo ðiÞKq

(15)

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

47

Table 3 Comparison of model performances in terms of the Nash-Sutcliffe coefficient of the original and modified versions of the VIC-3L model for applications at the 12 MOPEX watersheds ID

01608500 01643000 01668000 03179000 03364000 03451500 03054500 05455500 07186000 07378500 08167500 08172000 Mean a

New subsurface flow parameterization (i.e. Eq. (13))

ARNO subsurface flow parameterization

a priori (61–98)

Calibration (79–83)

Simulation (61–98)

a priori (61–98)

Calibration (79–83)

Simulation (61–98)

0.33 0.39 0.41 0.35 0.57 0.19 0.23 0.32 0.51 0.63 w0.00 w0.00 0.40a

0.74 0.80 0.69 0.57 0.67 0.69 0.61 0.51 0.64 0.82 0.62 0.60 0.66

0.68 0.73 0.56 0.61 0.63 0.66 0.66 0.51 0.74 0.70 0.45 0.31 0.60

0.31 0.40 0.36 0.30 0.48 0.27 0.21 0.29 0.50 0.54 w0.00 w0.00 0.37a

0.66 0.72 0.65 0.56 0.58 0.80 0.60 0.48 0.67 0.86 0.42 0.67 0.64

0.50 0.68 0.65 0.57 0.64 0.79 0.54 0.50 0.77 0.75 0.27 0.45 0.61

Averaged over the 10 watersheds excluding 08167500 and 08172000.

where qs(i) and qo(i) are the model simulated and observed streamflows at the ith time step; q
o is a mean observed streamflow over a period of length N. Nash-Sutcliffe coefficient is a commonly used coefficient in hydrology to evaluate the performance of a model with an emphasis on high flows. Moreover, since the focus of this study is on subsurface flows, another coefficient, log mean square error (called LOG hereinafter), is also used to evaluate the model performance at the twelve watersheds. This coefficient is essentially the same as the root mean square error (RMSE), which suggests a better model performance when its value is smaller. The LOG coefficient puts more emphasis on the low flow regime and is defined as, vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X LOG Z t (16) ðlog10 qs ðiÞKlog10 q0 ðiÞÞ2 N iZ1 The a priori, calibration, and validation results in terms of the Nash-Sutcliffe cofficient (i.e., equation (15)) and the LOG coefficient (i.e., equation (16)) corresponding to the original and modified versions of the VIC-3L model are given in Table 3 and Table 4, respectively. It can be seen from the tables that the performance of the modified version of VIC-3L is comparable to that of the version with the ARNO formulation. Specifically, the mean Nash–Sutcliffe

(LOG) coefficients of the modified model for the a priori, calibration, and the simulation periods are 0.40 (0.47), 0.65 (0.45), and 0.60 (0.46), respectively, while those of the original version of the model are 0.36 (0.43), 0.64 (0.41), and 0.61 (0.43), respectively. Note that the mean Nash–Sutcliffe coefficients for the a priori runs from both versions are the averages over 10 basins excluding 08167500 and 08172000, since the performances of the a priori simulations of the two basins are not better than the use of the long term average (i.e. the Nash–Sutcliffe coefficients are close to zero) for both versions. Also, the same treatment is applied to the mean LOG coefficients, because the values of the LOG coefficient for those two watersheds (i.e., 08167500 and 08172000) are unreasonably high. It can be seen from Table 3 and Table 4 that for the a priori simulations, the performances of the modified VIC-3L model are close to those of the original version for most of the basins, but with two less parameters to be estimated, which is helpful to mitigate the problems of parameter interactions. The result is also encouraging for the applications to ungauged basins. It is worth mentioning that from Table 4, the LOG coefficients of the calibration and simulation periods are not necessarily improved, for both versions, when compared to those of the a priori runs. This is understandable because in the calibration, equation

48

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

(15) is used as the objective function, which tunes the parameters to match high flows well, while the LOG coefficient focuses on evaluating model performances in the low flow regime. Besides the a priori simulations, both versions of VIC-3L are calibrated to find their ‘optimal’ parameters for each watershed. The calibrated parameter values corresponding to the two versions (not shown) suggest that the subsurface flow parameterizations could have significant impacts on the calibration of model parameters. For most of the watersheds, the calibrated parameters are different between the two versions. For example, for watershed 03054500, the calibrated parameters of the new version are 0.11, 1000.00, 0.37, 0.50, and 0.62 (i.e. b, a, two soil depths, and Kroute, respectively), while those of the original version (i.e. using ARNO parameterization) are 0.54, 0.97, 30.82, 0.28, 0.60, 3.21, and 3.22 (i.e. b, Ds, Dsmax, Ws, two soil depths, and Kroute, respectively). Such differences in the parameters sometimes introduce large differences in their corresponding hydrological processes, which in

turn affect the partition of water in the soil column. Thus, it is possible that the simulated subsurface flow from the two versions of VIC-3L differs significantly. For example, the model performances on streamflow evaluated by Eq. (15) for watershed 03054500 (which has a wet climate as shown in Table 1) with both versions are similar for the calibration period, but, the simulated subsurface flow based on ARNO are mostly zero, while the one based on the new parameterization can be much higher than zero (figure not shown). Parameter interaction seems to play an important role in this example. Although the values of subsurface flow parameters in both versions indicate the subsurface flow component would be significant in this watershed, the calibrated value of the b parameter decreases from 0.54 in the original version to 0.11 in the modified version so that the partition of the water budget changes. A larger value of the b parameter in the original version indicates that a larger fraction of precipitation is partitioned into surface runoff. Hence, there is not much water available from infiltration to support subsurface flow

Fig. 2. Comparison between the simulated subsurface flow from the new parameterization (i.e. dashed line) and that from the ARNO parameterization (i.e. solid line) with calibration parameters for watershed 01608500 during two periods, calibration period (i.e. 01/01/1979– 12/31/1983) and validation period (i.e. 01/01/1971–12/31/1983).

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

generation in the original version. Also, the data used for calibration seem to be pretty noisy for this watershed, in which the precipitation and streamflow events are not consistent with each other sometimes. The new subsurface flow parameterization, on the other hand, seems to be more robust than the original one in this example, so that the problem of parameter interaction is not as significant as that in the original version because of the reduction of number of parameters. In fact, the performance of the new version for the simulation period (i.e. EZ 0.66) is much better than that for the original version (i.e. EZ0.54) for this watershed. On the other hand, there are also cases in which the simulated subsurface flows from both parameterizations are comparable. For example, Fig. 2 shows the comparison between simulated subsurface flows based on the new parameterization and those on the ARNO parameterization for watershed 01608500 for both the calibration period and part of the validation period. The patterns of the subsurface flows from both versions are similar. A closer examination of the performances of the two different subsurface flow parameterizations for individual watersheds shows that the two parameterizations favor different watersheds. For example, the performance of the original version at watershed 03451500 is significantly better than that of the new version; while the performance of the new version at watershed 01643000 is significantly better than that of the original version. The underlying reason is complicated (e.g. may be due to the combination of the subsurface flow parameterization and the interactions among the model parameters that are calibrated) and needs to be further investigated. Special attention is paid to the newly introduced parameter a (i.e. the anisotropic ratio between lateral and vertical hydraulic conductivities). Large variations in the calibrated a values are observed, but for each watershed a varies little (see discussion in the next section). The large variation in a values is consistent with the values reported by Kumar (2004). As suggested by Eq. (12), the subsurface flow is determined by not only a, but also the saturated vertical hydraulic conductivity, the depths of the soil layers (which are calibrated in this study), the relative soil moisture content, and the Brooks–Corey coefficient a. Intuitively, we can expect the three driest

49

watersheds (i.e. 05455500, 08167500, and 08172000) to have the smallest a values, which is consistent with our calibrated a values in this study (i.e. 0.11, 0.60, 1.06, respectively). This implies that the subsurface component has limited contribution to the hydrographs in these watersheds, which is also observed from the original version with the ARNO parameterization. The calibrated values of Ds and Ws in ARNO for these three watersheds are close to 0 and 1, respectively, which implies that the subsurface flow would only be generated under cases when the third layer is saturated, otherwise, it is close to zero. In summary, the preliminary results from this study suggest that the new one-parameter subsurface flow parameterization can capture reasonably well the subsurface flow behavior that is represented by the three-parameter ARNO subsurface flow formulation. The reduction of the number of parameters, the physical soundness, and the comparable and reasonable performance of the new subsurface flow parameterization make it more suitable for applications to both ungauged basins and basins with limited observations. 3.3. Model parameter sensitivity analyses with Monte Carlo simulations As discussed in Section 3.2, the modified version of VIC-3L (i.e. using Eq. (13)) contains five parameters that need to be calibrated. Among the five parameters, the two soil depths have clear physical meanings and can be estimated in principal based on field measurements. The routing parameter, Kroute, could be replaced by a more sophisticated routing scheme. Therefore, the two parameters that require calibration are the effective soil moisture capacity shape parameter b in Eq. (1) and the antistrophic factor of hydraulic conductivity a in Eq. (12). To investigate the sensitivities and uncertainties of model simulated streamflow to the two parameters, we carry out Monte Carlo experiments over three watersheds. We use the Monte Carlo simulation approach to generate N1 successive parameter sets, with p parameters in each set. The N1 parameter sets are drawn from a joint distribution of the p parameters. Then N1 realizations of model simulated streamflows are obtained from VIC-3L.

50

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

Fig. 3. The model parameter response surfaces for different watersheds based on the Monte Carlo simulations. (a) Model parameter response surface of watershed 07378500. (b) Model parameter response surface of watershed 01608500. (c) Model parameter response surface of watershed 03054500.

In this study, p is taken to be two, which are parameters b and a, and they vary over the ranges of [0, 5] and [0.1, 1000], respectively. By dividing the two ranges into 20 and 50 equal intervals, 20 values of b and 50 values of a are generated assuming uniform distributions. Therefore, within the two-dimensional parameter space, 1000 pairs (i.e. N1Z1000) of parameter values are sampled for each watershed. For each of the 1000 pairs of parameter values, a simulation of the modified VIC-3L model is carried out for the period of 01/01/1960–12/31/1998 while the rest of the three model parameters (i.e. the two soil depths and Kroute) are kept at their values obtained from the calibration. The effectiveness of the simulations is evaluated by the Nash–Sutcliffe coefficient (i.e. Eq. (15)). Since the sample size of the two-dimensional parameter space used in this study is large (i.e. 1000 points in which 20 and 50 points are sampled for the b parameter and a, respectively), only three watersheds, i.e. 01608500, 07378500, and 03054500, are used to examine the sensitivity of model simulated streamflows to the parameters. The climate conditions of these three watersheds range from arid to moist

conditions and are representative of the climate variations among the 12 watersheds. The response surfaces generated from the sensitivity analyses corresponding to the three watersheds are shown in Fig. 3, which suggests that the model simulated streamflow is not sensitive to the magnitude of a for a fixed watershed, although the variations of a could be large from one watershed to the other. For example, on the response surface for watershed 07378500, the parameter values associated with the highest Nash– Sutcliffe coefficient (0.74) are 230.08 and 0.13 for a and b parameter, respectively. When the b parameter is fixed to be 0.13, and a varies from 10 to 1000, the Nash–Sutcliffe coefficient changes within the range of [0.62, 0.74]. As a comparison, when a is fixed to be 230.08, and the b parameter varies from 0.13 to 4.88, the Nash–Sutcliffe coefficient changes within the range of [0.31, 0.74]. Similar patterns are observed on the response surfaces for the other two watersheds, as shown in Fig. 3. Compared to the large sensitivities of the model simulated streamflows to the b parameter, the sensitivity to the a parameter is relatively small. Thus, it can be fixed to a preliminary estimated value

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

based on field measurements. This is especially helpful for applications to ungauged basins, where no data or very limited data are available for calibration. However, different values of the b parameter could result in relatively large differences in the model simulated streamflows. Therefore, we focus our investigations on the uncertainties associated with the b parameter in Section 3.4. 3.4. Model parameter uncertainties and their impacts on model simulations In this section, a procedure of quantifying uncertainties associated with the b parameter is described and results are discussed. The VIC-3L b parameter is an important parameter, which controls the shape of the effective soil moisture capacity distribution over a study area. Due to its conceptual nature and limited observations, the b parameter is generally difficult to be estimated directly from available field measurements, and is thus determined in general through model calibration. The role of the b parameter on model simulated streamflows is investigated from two aspects: (1) the accuracy of the estimation of the b parameter using the SCE-UA searching algorithm; and (2) the impacts of uncertainties associated with the estimations of the b parameters on model simulated streamflows. Again, Monte Carlo simulations are carried out by sampling the b parameter within its feasible range of [0, 5], while other parameters are fixed at their calibrated values, including parameter a since it is not sensitive for a fixed watershed. Simulations with Nash–Sutcliffe coefficients over the calibration period less than a prescribed level (e.g. 0.5) are removed from further analyses. If the acceptable sample size of the b parameter is too small, more sampling with smaller increments over [0, 5] will be taken to ensure a large enough sample size (i.e. greater than 100 samples) for analyses. For the selected simulations, the mean VIC-3L simulated streamflow that corresponds to a range of b parameter values is given by

q
is Z

Na 1 X qi ðkÞ Na kZ1 s

(17)

51

where qis ðkÞ is the simulated streamflow at the ith time step based on the kth acceptable parameter set; i2(1, .,N), where N is the length of the simulation period, and k2(1,.,N), where N a is the number of acceptable parameter sets based on the Nash–Sutcliffe coefficients evaluated over the calibration period. Moreover, an upper bound qisðupperÞ and a lower bound qisðlowerÞ associated with the 95% interval can be calculated based on the selected simulations. From a probabilistic point of view, the joint probability of the simulations and the observations of the streamflow can be factored into conditional and marginal probabilities expressed as Pðst ; ok Þ Z Pðst jok ÞPðok Þ

(18)

where ok (kZ1,.,K) denotes the target event (e.g. the observed streamflow is within certain range), K represents the number of different ranges of events that we are interested in; P(ok) is then the probability that the target event occurs, which can be calculated from the relative frequency of the target event in the observed time series of the streamflow; St is a binary variable denoting whether the target event occurs or not (i.e. if the target event occurs, StZ1; otherwise, StZ0) at time t in the corresponding simulation, so that P(stjok) is the conditional probability denoting the chance of occurrence of st in the simulation when the target event ok occurs in the observation. In this paper, we are interested in evaluating P(stZ1jok) based on the times series of the ensemble mean, upper and lower bounds, and optimal simulation, and observation. Since our purpose is to assess the uncertainty introduced by model parameters to model simulations, which is different from the analyses originated from the ensemble weather forecasting (e.g. Wilks, 1995; see Georgakakos et al. (2004) for its application in hydrology), in which the performance of ensemble forecasts is evaluated, we will not use the standard representation format (e.g. the reliability diagram, the Brier score, and the rank histogram) that is often adopted in the ensemble forecast studies. Instead, we carry out our analysis with a procedure described as follows. Let us define Q1 and Q2 to be the prescribed upper and lower threshold values obtained from observed ðiÞ streamflow qðiÞ o , and Qs to be the model simulated

52

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

daily streamflow, where i represents the index in the observed and simulated time series. Under this circumstance, the target event o k becomes Q1 % qðiÞ o ! Q2 , and the corresponding target event in the simulation becomes Q1 % QðiÞ s ! Q2 . By varying Q1 and Q2 through the entire range of observed flows of the study watershed, we will be able to evaluate a model’s ability to simulate observed streamflows in different flow ranges by computing the conditional ðiÞ probability PðQ1 % QðiÞ s ! Q2 jQ1 % qo ! Q2 Þ, which is the same as P(stjok) in Eq. (18). Detailed steps of our procedure are as follows. (1) Calculate the histogram of observed streamflows. (2) Evaluate the observed flow time series to find the 50, 75, 90, and 99 percentiles based on the histogram obtained in step (1). These flows are then used as Q1 and/or Q2 in the analyses. In addition, the minimum and maximum flow values Qmin and Qmax in the time series are also selected as the smallest and largest flow thresholds in this analysis. (3) Compute the conditional probability PðQ1 % ðiÞ QðiÞ s ! Q2 jQ1 % qo ! Q2 Þ by counting the number of streamflow records bounded between Q1 and Q2 in the observed streamflow time series qðiÞ o , and the number of their corresponding model simulated ðiÞ flows QðiÞ s , that are also within this range. Qs can be taken as the mean, upper, and lower bounds of the simulated streamflow based on Monte Carlo simulations over the range of b parameter values, or the simulated streamflow corresponding to the optimized b parameter from the calibration. Q1 and Q2 are the observed streamflows corresponding to two neighboring percentiles of interest (e.g. Q1Z Q0.05, Q2ZQ0.75 where Q0.50 and Q0.75 correspond to the observed streamflows at 50 and 75 percentiles, respectively). Q1 and Q2 are calculated from step (2) above. Impacts and uncertainties of the b parameter on VIC-3L model simulated streamflows are evaluated based on the conditional probability PðQ1 % QðiÞ s ! Q2 jQ1 % qðiÞ o ! Q2 Þ described above. 3.4.1. Uncertainties of model parameters VIC-3L model performances expressed by the Nash–Sutcliffe coefficient are evaluated for

Fig. 4. The Nash–Sutcliffe coefficients versus b parameter values with other parameters fixed at their calibrated values for the 12 watersheds. The symbols ‘B‘ and ‘,‘ represent, respectively, the maximum values of the Nash–Sutcliffe coefficient for the calibration period and simulation period.

the calibration and entire simulation periods based on the different b parameter values across a large range. The robustness of b parameter estimation by the SCE-UA method is first investigated. Fig. 4 shows the VIC-3L model performances as a function of the b parameter sampled within the range of [0, 5]. From Fig. 4, we can see the following: (1) The automatic calibration algorithm may fail to find the optimal b parameter because of the complicated model response within a certain range of the b parameter values. For example, the model performance (i.e. the Nash–Sutcliffe

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

53

Table 4 Comparison of model performances in terms of the LOG coefficient of the original and modified versions of the VIC-3L model for applications at the twelve MOPEX watersheds ID

01608500 01643000 01668000 03179000 03364000 03451500 03054500 05455500 07186000 07378500 08167500 08172000 mean

New subsurface flow parameterization (i.e., equation 13)

ARNO subsurface flow parameterization

a priori (61–98)

Calibration (79–83)

Simulation (61–98)

a priori (61–98)

Calibration (79–83)

Simulation (61–98)

0.52 0.38 0.43 0.44 0.37 0.83 0.49 0.49 0.44 0.36 w w 0.47*

0.28 0.50 0.59 0.45 0.43 0.27 0.30 1.07 0.44 0.29 0.45 0.33 0.45

0.31 0.49 0.56 0.46 0.43 0.28 0.34 1.07 0.43 0.33 0.49 0.36 0.46

0.37 0.24 0.31 0.34 0.33 0.65 0.48 0.51 0.75 0.37 w w 0.43*

0.29 0.66 0.50 0.36 0.37 0.15 0.59 0.52 0.40 0.29 0.51 0.29 0.41

0.32 0.65 0.52 0.40 0.37 0.13 0.60 0.50 0.50 0.33 0.56 0.30 0.43

* Averaged over the ten watersheds excluding 08167500 and 08172000.

coefficient) at watersheds 01668000 (i.e. Fig. 4(c)) shows a large variation for the values of the b parameter within the range of 0–1.5. As a consequence, the automatic calibration algorithm failed to find the optimal b parameter value for this watershed in the calibration period (i.e. 1979– 1983). The optimal b parameter value obtained from the calibration process is 0.4343, but is 0.6813 from the Monte Carlo simulations for the same period. For most of the watersheds, however, the shuffled complex evolution (SCEUA) optimization algorithm (Duan et al., 1992) is generally robust in finding the optimal b parameter values. Table 5 lists the magnitudes of the b parameter obtained from the calibration process and the ones from the Monte Carlo simulations for both of the calibration and validation periods. (2) Multiple ‘optimal’ b parameter values could exist for some watersheds, even when the parameter space is one-dimensional (i.e. the b parameter) as in this study. For example, the response curve of watershed 03179000 (Fig. 4(d)) suggests that over a range of 0.5–2.0, the performances of the model are comparably good, with the Nash–Sutcliffe coefficients close to 0.6 for both the calibration and simulation periods. Similar behaviors are observed in Fig. 4(f). Under such circumstances,

it is hard for the calibration algorithm to identify a unique optimal b parameter value owing to equifinality (e.g. Beven, 2000). (3) The optimal b parameter value for the simulation period could be different from that for the calibration period. For example, the optimal b parameter values corresponding to the calibration and simulation periods for watershed 03179000 (i.e. Fig.4(d)) are 1.24 and 0.87, respectively. Table 5 lists the optimal values obtained from the calibration, Monte Carlo simulation during the calibration and simulation periods, respectively, for the 12 watersheds. Also, the minimum and maximum values of the b parameter among the Na number of simulations for each watershed are listed in Table 5. Our preliminary results suggest that we may not be able to find a unique set of optimal parameters, even with a one-dimensional parameter space as illustrated in this study, based on the current understanding of hydrological processes and available automatic calibration techniques. Also, even if our calibration technique would allow us to find a global optimal for the calibration period, it is likely that we cannot find a time-invariant model parameter set that is optimal for all the different time periods. Therefore, it is critical to recognize and represent the impacts of such

54

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

Fig. 5. Histograms of the observed streamflow for the 12 MOPEX watersheds.

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61 Table 5 Selected number of model simulations for the calculation of mean responses and their confidence intervals for the 12 watersheds Watershed

Na

Watershed

Na

01608500 01643000 01668000 03179000 03364000 03451500

122 322 394 355 177 288

03054500 05455500 07186000 07378500 08167500 08172000

157 181 183 400 100 108

uncertainties that are associated with model parameters on the model simulations in our studies. 3.4.2. Impacts of model parameter uncertainties on model simulations Impacts of model parameter uncertainties on model simulations are investigated by applying conditional probability methodology to model simulations drawn from Na number of simulations. The number, Na, is determined based on the values of the Nash Sutcliffe coefficients. Specifically, the VIC-3L simulations are selected in general for conducting the analyses if their Nash–Sutcliffe coefficients are greater than 0.5 over the calibration period (i.e. 1979–1983). For the three dry watersheds, 08167500, 08172000, and 05455500, Na is determined based on the number of simulations over the calibration period whose Nash–Sutcliffe coefficient are greater than 0.3 for watersheds 08167500 and 08172000, and greater than 0.4 for watershed 05455500, since if the magnitude of 0.5 were used, there would not be enough simulations available for conducting the uncertainty analyses. The number, Na, of accepted simulations based on data from the calibration period (i.e. 1979–1983), for the 12 watersheds are shown in Table 6. In Table 5, the minimum and maximum numbers of Na are 100 and 400, respectively, which are statistically sufficient for conducting our analyses. Table 5 shows a large range of variation for the b parameter values, i.e. the minimum and maximum values of the b parameter that correspond to the Na number of simulations over the calibration period for each watershed are quite large (e.g. from 0.01 to 5.00 for watershed 01668000). With the Na sets of the model parameters, Monte Carlo simulations over the entire period (i.e. 01/01/1961–12/31/1998) are conducted for each watershed, and the mean response and

55

the responses corresponding to the upper and lower bounds can be obtained for every time step over the entire simulation period (i.e. 01/01/1961–12/31/ 1998). The corresponding Nash-Sutcliffe coefficients of the mean responses for the twelve watersheds are shown in Table 8, which suggest that the selected Na simulation provides an ensemble with reasonable performance. The new responses can be used to study the impacts of parameter uncertainties on model simulations based on the conditional probability approach. One critical step in applying the conditional probability approach is to select Q1 and Q2 values ðiÞ used in PðQ1 % QðiÞ s ! Q2 jQ1 % qo ! Q2 Þ objectively. Because the Q1s and Q2s are to be determined based on the percentiles of the observed streamflows, objective selection of the threshold flow values for each watershed is important for conducting meaningful analysis. For example, if the threshold is set to be zero (i.e. to include all of the flows in the analyses), it will then include a significantly large number of small daily flows over the 40-year period. Fig. 5 shows the histograms of the observed streamflows with the threshold set to be zero for the 12 MOPEX watersheds. From Fig. 5, it can be seen that most of the observed streamflow values are very small (i.e. close to zero). From a hydrological modeling point of view, these very small flows are less interesting. Therefore, in this study, we truncate the time series so that flows close to zero are excluded from further analyses. Specifically, the time series of observed streamflows with flow values lower than 75 percentile (i.e. Q0.75) are excluded. The magnitude of Q0.75 of each watershed and its minimum and maximum flows are listed in Table 7. With this operation, the remaining observed and simulated streamflows obtained from the mean response, its bounds, and the simulation with calibrated parameters (called ‘optimal’ hereinafter) form the new time series with approximately 3500 flow values (points) in each case, which are sufficient for statistical analyses. The new time series is then used to calculate the flow bounds (i.e. Q1 and Q2) and the condiðiÞ tional probabilities PðQ1 % QðiÞ s ! Q2 jQ1 % qo ! Q2 Þ following the procedure described earlier. The results of the probabilistic analyses for the 12 watersheds are shown in Fig. 6, in which the horizontal axis represents the magnitudes of flow bounds (i.e. Q1

56

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

Table 6 Flow values corresponding to the minimum, 75 percentile, and maximum of the observed streamflow time series of the 12 watersheds Watershed

Qmin (mm/day)

Q0.75 (mm/day)

Qmax (mm/day)

Watershed

Qmin (mm/day)

Q0.75 (mm/day)

Qmax (mm/day)

01608500 01643000 01668000 03179000 03364000 03451500

0.03 0.02 0.00 0.03 0.05 0.34

1.01 1.20 1.13 1.30 1.12 2.55

93.12 85.57 49.84 38.12 27.12 31.99

03054500 05455500 07186000 07378500 08167500 08172000

0.01 0.00 0.02 0.20 0.00 0.05

2.46 0.65 0.65 1.33 0.29 0.48

52.50 36.77 87.65 77.49 54.96 102.02

and Q2) based on the truncated observed streamflow time series. Values of Qmin (i.e. the 75 percentile of streamflow in the original time series), Q0.50, Q0.75, Q0.90, Q0.99, and Qmax (i.e. the maximum value of streamflow for a watershed during the period of 01/01/ 1961–12/31/1998) are used to serve as the Q1 and Q2 consequently. That is, we use (Q1, Q2) to represent (Qmin, Q0.50), or (Q0.50, Q0.75), or (Q0.75, Q0.90), or (Q0.90, Q0.99), or (Q0.99, Qmax), and use group 1, 2, 3, 4, and 5, respectively, to represent each range. The vertical axis represents the conditional probabilities of ðiÞ PðQ1 % QðiÞ s ! Q2 jQ1 % qo ! Q2 Þ. From the probabilistic analyses and Fig. 6 we can see the following behaviors. (1) Both the optimal simulation and mean response suggest that the VIC-3L model can capture flows in different ranges reasonably well. The probability of the model to capture flows higher than Q0.99 (i.e. group 5) is pretty small for half of the 12 watersheds. The inability of the model to capture extremely high flows is expected because there are very small number of streamflows that fall into group 5 among the 12 watersheds. Therefore, if there are a few model simulated flows that fall outside the range of group 5, then the impact on the conditional probability will be significant. On the other hand, the upper bound responses result in higher probabilities for group 5 than the mean response and optimal simulation (e.g. Fig. 6(b), (c), (e), and (i)) as expected. (2) The mean response and the responses corresponding to the upper and lower bounds shown in Fig. 6 provide good uncertainty estimates of the model simulated streamflows contributed by the uncertainties in the b parameter. Results obtained here suggest that the model simulated streamflow can

be used with reasonable capability of capturing the observed streamflows across a wide range of the b parameter (see Table 5) if the VIC b parameter cannot be well estimated. Caution needs to be paid to the high flows, and the model simulation corresponding to the upper bound may have less uncertainty. Fig. 6 suggests great potential of applying the VIC-3L model to ungauged basins where no data or very limited data are available for estimating the b parameter through model calibration. From Fig. 6, it can also be seen that for most of the watersheds, the probabilities corresponding to the mean responses and ‘optimal’ simulations are bounded within the probabilities corresponding to the responses associated with the upper and lower bounds. For the low flow ranges, the probabilities that the Monte Carlo simulations corresponding to the upper bound (lower bound) to capture the specified observed flow ranges, are lower (higher) than those of the mean responses and the optimal simulations, and the vice versa for the higher flow ranges. Moreover, there are cases where the probabilities corresponding to the optimal simulations are actually lower than those Table 7 Nash–Sutcliffe coefficients of the mean responses for the 12 watersheds during the period of 01/01/1961–12/31/1998 Watershed

N–S

Watershed

N–S

01608500 01643000 01668000 03179000 03364000 03451500

0.6885 0.7348 0.5554 0.6097 0.6331 0.6569

03054500 05455500 07186000 07378500 08167500 08172000

0.6289 0.5437 0.7431 0.7035 0.4580 0.3114

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

57

corresponding to the mean responses and their bounds (e.g. Fig. 6(h) watershed 05455500). This is because (i) the automatic optimization algorithm may fail to find the optimal parameter set over the calibration period, and (ii) the optimal parameter sets for the calibration period are not necessarily the ‘optimal’ sets for the validation period. Therefore, it is very useful to carry out probabilistic analyses as presented here to represent model simulation uncertainties that are associated with parameter uncertainties even if there are data available to estimate model parameter (e.g. the b parameter) through calibrations. (3) A poorly estimated b parameter value may introduce more uncertainty over the high flow ranges than low flow ranges, as illustrated in Fig. 6. For example, for watershed 01608500 (i.e. Fig. 6(a)), the uncertainty bounds for groups 1 and 2 are quite small, while the differences between them for groups 4 and 5 are large. This is partially because b is a parameter that affects surface runoff generation, which contributes mainly to peak flows other than low flows.

4. Conclusions

Fig. 6. Probabilities that the simulated streamflows (i.e. the mean responses, their upper bounds and lower bounds based on the Monte Carlo simulations, the optimal simulations based on calibrated parameters) fall into the intervals when the observations fall into the same intervals for the 12 watersheds. Values corresponding to Qmin, Q0.50, Q0.75, Q0.90, and Q0.99 are indicated on the x-axes.

In this paper, uncertainties associated with the VIC-3L model parameters and the impacts of the parameter uncertainties on model simulated streamflow are investigated. An alternative subsurface flow parameterization based on the concepts of kinematic wave and the hydrologic similarity is presented and implemented into VIC-3L. The new subsurface flow parameterization is applied to the 12 MOPEX watersheds and is compared to the original ARNO formulation used in the VIC-3L model. Sensitivities of the model simulated streamflows to the VIC-3L model parameters with the new subsurface flow parameterization are studied using the Monte Carlo simulations. Impacts of the VIC-3L model b parameter on simulated streamflows are evaluated using a conditional probability method. Preliminary conclusions are summarized as follows.

58

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

(1) The alternative subsurface flow parameterization presented in the study can successfully capture the baseflow component of the streamflow in general with only one parameter for calibration (i.e. parameter a) instead of three (i.e. parameters Ds, Ws, and Dsmax) in the ARNO formulation. This parameter is not sensitive to the model simulated streamflows, and existing values from the literature can be used. In addition, this parameter can also be estimated based on field measurements. Thus, the new subsurface flow parameterization makes it easier for applications to ungauged basins, where no data or very limited data are available for calibration. In addition, this parameterization is general and can be implemented into other models if needed. (2) With the new subsurface flow parameterization, the problem of parameter interactions in the modified VIC-3L model is partially mitigated owing to the reduction of the number of model parameters. Therefore, the modified VIC-3L model may be more robust than the original one, especially when the calibration data are noisy. (3) Results from the probabilistic approach show that the VIC-3L model with the new subsurface flow parameterization can capture flows in different ranges, except for flows higher than the 99th percentile of the observed streamflow. Under such circumstances, the model responses corresponding to the upper bound case of the Monte Carlo simulations provide a better estimation than the optimal simulation or the mean response for each watershed as expected. Thus, the upper bound case can be considered as a good candidate for predictions or simulations for very large flows. (4) For applications to ungauged basins, the mean responses and their uncertainty bounds based on Monte Carlo simulations may provide good streamflow predictions. (5) Sensitivity analyses show that the simulations of the modified VIC-3L model are more sensitive to the VIC shape parameter (i.e. the b parameter) than to the baseflow parameter a. Moreover, impacts of the uncertainties due to the b parameter on simulated streamflows are

generally larger for high flows than for low flows. (6) A probabilistic approach is used to estimate the model performance over different flow ranges in this study. This approach can be used as a complement to other methods in evaluating model performance. To apply this approach, it is crucial to select reasonable thresholds for analyses. Because in this study, the sampling of only one- to two-dimensional parameter spaces is needed, we adopted the importance sampling strategy. However, if the sampling of a higher dimensional parameter space is needed, the MCMC sampling strategy is recommended.

Acknowledgements We would like to thank the two anonymous reviewers for their helpful suggestions. The research reported herein was supported by the National Oceanic and Atmospheric Administration under grant DG133R-03-SE-0680 to the University of California at Berkeley. Partial funding to the first author provided by the Berkeley Atmospheric Science Center is also greatly appreciated.

References Ambroise, B., Beven, K., Freer, J., 1996. Toward a generalization of the TOPMODEL concepts: Topographic indices of hydrological similarity.Benning, R.G., 1994. Towards a new lumped parameterisation at catchment scale. Master Thesis, University of Wageningen, The Netherlands. 10–45. Benning, R. 1994. Towards a new lumped parameterization at catchment scale. Master Thesis, University of Wageningen. Beven, K., 1981. Kinematic Subsurface Stormflow. Water R. Res. 17 (5), 1419–1424. Beven, K.J., 1984. Infiltration into a class of vertically non-uniform soils. Hydrol. Sci. J. 29, 425–434. Beven, K.J., 1989. Changing ideas in hydrology: the case of physically-based models. J. Hydrol. 105, 157–172. Beven, K.J., 1997. TOPMODEL: A critique. Hydrol. Process. 11, 1069–1085. Beven, K.J., 2000. Uniqueness of place and representation of hydrological processes. Hydrol. Earth Syst. Sci. 4 (2), 203–213. Beven, K.J., Binley, A.M., 1992. The future of distributed models: Model calibration and uncertainty prediction. Hydrol. Process. 6, 29–44.

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61 Beven, K.J., Kirkby, M.J., 1979. A physically based variable contributing area model of catchment hydrology. Hydrol. Sci. Bull. 24, 43–69. Beven, K.J., Kirkby, M.J., Schofield, N., Tagg, A.F., 1984. Testing a physically-based flood forecasting model (TOPMODEL) for three UK. Catchments. J. Hydrol. 69, 119–143. Beven, K.J., Lamb, R., Quinn, P., Romanowicz, R., Freer, J., 1995. TOPMODEL. In: Singh, V.P. (Ed.), Computer models of watershed hydrology. Water Resource Publications, Colorado, pp. 627–668. Boyle, D.P., Gupta, H.V., Sorooshian, S., 2000. Toward improved calibration of hydrologic models: Combining the strengths of manual and automatic methods. Water R. Res. 36 (12), 3663– 3674. Brutsaert, W., 1994. The unit response of groundwater outflow from a hillslope. Water R. Res. 30, 2759–2763. Brutsaert, W., Lopez, J.P., 1998. Basin-scale geohydrologic drought flow features of riparian aquifers in the southern Great Plains. Water Resources Research 34 (2), 233–240. Chen, J., Kumar, P., 2001. Topographic influence on the seasonal and inter-annual variation of water and energy balance of basins in North America. J. Climate 14, 1989–2014. Ciarapica, L., Todini, E., 2002. TOPKAPI: a model for the representation of the rainfall-runoff process at different scales. Hydrol. Process. 16, 207–229. Clapp, R.B., Hornberger, 1978. Empirical equations for some soil hydraulic properties. Water R. Res. 14, 601–604. Cherkauer, K.A., Lettenmaier, D.P., 1999. Hydrologic effects of frozen soils in the upper Mississippi River basin. J. Geophys. Res. 104 (19), 599–610. De Grosbois, E., Hooper, R.P., Christopherson, N., 1988. Christopherson, A multisignal automatic calibration methodology for hydrochemical models: a case study of the birkenes model. Water R. Res. 24, 1299–1307. Duan, J., Miller, N.L., 1997. A generalized power function for the subsurface transmissivity profile in the TOPMODEL. Water R. Res. 33 (11), 2559–2562. Duan, Q., Sorooshian, S., Gupta, V.K., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water R. Res. 28 (4), 1015–1031. Duan, Q., Gupta, V.K., Sorooshian, S., 1993. A shuffled complex evolution approach for effective and efficient global minimization. J. Optimization Theory Appl. 76 (3), 501–521. Duan Q., J. Schaake ,V. Andreassian, S. Franks, H.V. Gupta, Y.M. Gusev, F. Habets, A. Hall, L. Hay, T. Hogue, M. Huang, G. Leavesley, X. Liang, O.N. Nasonova, J. Noilhan, L. Oudin, S. Sorooshian, T. Wagener, E.F. Wood, 2005. Model Parameter Estimation Experiment (Mopex): Overview And Summary Of The Second And Third Workshop Results, Journal of Hydrology, x-ref: DOI 10.1016/j.hydrol.2005.07.037. Eng, K., Brutsaert, W., 1999. Generality of drought flow characteristics within the Arkansas River Basin. J. Geophys. Res. 104 (D16), 435–441. Francini, M., Pacciani, M., 1991. Comparative Analysis of several conceptual rainfall-runoff models. J. Hydrol. 122, 161–219.

59

Freer, J., Beven, A.M., Ambroise, B., 1996. Bayesian estimation of uncertainty in runoff prediction and the value of data: an application of the GLUE approach. Water R. Res. 32, 2161–2173. Georgakakos, K.P., Seo, D.-J., Gupta, H.V., Schaake, J., Butts, M., 2004. Towards the characterization of streamflow simulation uncertainty through multimodel ensembles. Journal of Hydrology, 289, 222–241. Gupta, H.V., Sorooshian, S., Yapo, P.O., 1998. Toward improved calibration of hydrologic models: Multiple and noncommensurable measures of information. Water R. Res. 34 (4), 751–763. Gupta, H.V., Bastidas, L.A., Sorooshian, S., Shuttleworth, W.J., Yang, Y.L., 1999. Parameter estimation of a land surface scheme using multicriteria methods. J. Geophys. Res. 104 (D16), 491–503. Gupta, H.V., Bastidas, L., Vrugt, J.A., Sorooshian, S., 2003. Multiple criteria global optimization for watershed model calibration. In: Duan, Q. et al. (Ed.), Calibration of Watershed Models Water Science Applied Series, vol. 6. AGU, Washington, DC, pp. 125–132. Hamlet, A.F., Lettenmaier, D.P., 1999a***. Effects of climate change on hydrology and water resources in the Columbia River basin. J. Am. Water R. Assoc. 35 (6), 1597–1623. Hamlet, A.F., Lettenmaier, D.P., 1999b. Columbia River streamflow forecasting based on ENSO and PDO climate signals. J. Water R. Plan. Manag.-Asce 125 (6), 333–341. Hastings, W.K., 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97–109. Huang, M., Liang, X., Liang, Y., 2003. A transferability study of model parameters for VIC land surface scheme. J. Geophys. Res. 108 (D22), 8864. Iorgulescu, I., Musy, A., 1997. Generalization of TOPMODEL for a power law transimissivity profile. Hydrol. process. 11, 1353– 1355. Kirkby, M., 1975. Hydrograph modelling strategies. In: Peel, R., Chisholm, M., Haggett, P. (Eds.), Processes in Physical and Human Geography. Heinemann, London, pp. 69–90. Klepper, O., Scholten, H., van de Kamer, J.P.G., 1991. Prediction Uncertainty in an Ecological Model of the Oosterschelde Estuary. J. Forecasting 10, 191–209. Kuczera, G., Parent, E., 1998. Monte Carlo assessment of parameter uncertainty in conceptual catchment models: The Metropolis algorithm. J. Hydrol. 211, 69–85. Kumar, P., 2004. Layer averaged Richard Equation with lateral flow. Advances in Water Resources, 27(5) 522–532. Leung, L.R., Hamlet, A.F., Lettenmaier, D.P., Kumar, A., 1999. Simulations of the ENSO hydroclimate signals in the Pacific Northwest Columbia River basin. Bull. Am. Meteorol. Soc. 80 (11), 2313–2329. Liang, X., Lettenmaier, D.P., Wood, E.F., Burges, S.J., 1994. A simple hydrologically based model of land surface water and energy fluxes for general circulation models. J. Geophys. Res. 99 (D7), 415–428. Liang, X., Lettenmaier, D.P., Wood, E.F., 1996a. A onedimensional statistical dynamic representation of subgrid spatial variability of precipitation in the two-layer variable infiltration capacity model. J. Geophys. Res. 101 (D16), 403–421 p. 422.

60

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61

Liang, X., Wood, E.F., Lettenmaier, D.P., 1996b. Surface soil moisture parameterization of the VIC-2L model: evaluation and modifications. Global Planetary Change 13, 195–206. Liang, X., Wood, E.F., Lettenmaier, D.P., 1999. Modeling ground heat flux in land surface parameterization schemes. J. Geophys. Res. 104 (D8), 9581–9600. Liang, X., Xie, Z., 2001. A new surface runoff parameterization with subgrid-scale soil heterogeneity for land surface models. Adv. Water Res. 24 (9–10), 1173–1193. Liang, X., Xie, Z., Huang, M., 2002. A new parameterization for surface and groundwater interactions and its impact on water budgets with the variable infiltration capacity (VIC) land surface model. J. Geophys. Res. 108 (D16), 8613. doi:10. 1029/2002JD003090. Liu, Z., Todini, E., 2002. Towards a comprehensive physicallybased rainfall-runoff model. Hydrol. Earth Syst. Sci. 6 (5), 859–881. Mengelkamp, H.T., Warrach, K., Raschke, E., 1999. SEWAB - a parameterization of the Surface Energy and Water Balance for atmospheric and hydrologic models. Adv. Water Res. 23 (2), 165–175. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E., 1953. Equations of state calculations by fast computing machines. J. Chem. Phys. 21, 1087–1091. Mitchell, K.E., Lohmann, D., Houser, P.R., Wood, E.F., Schaake, J.C., Robock, A., Cosgrove, B.A., Sheffield, J., Duan, Q., Luo, L., Higgins, R.W., Pinker, R.T., Tarpley, J.D., Lettenmaier, D.P., Marshall, C.H., Entin, J.K., Pan, M., Shi, W., Koren, V., Meng, J., Ramsay, B.H., Bailey, A.A., 2003. The multi-institution North American Land Data Assimilation System (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrological modeling system. J. Geophys. Res. 109, D07S90. doi: 10.1029/2003JD003823. Nijssen, B., O’Donnell, G.M., Hamlet, A.F., Lettenmaier, D.P., 2001a. Hydrologic sensitivity of global rivers to climate change. Climatic Change 50 (1–2), 143–175. Nijssen, B., O’Donnell, G.M., Lettenmaier, D.P., Lohmann, D., Wood, E.F., 2001b. Predicting the discharge of global rivers. J. Climate 14 (15), 3307–3323. Niu, G.-Y., Yang, Z.L., 2003. The versatile integrator of surface atmospheric processes (VISA): part 2, evaluation of three topography-based runoff schemes. Global Planetary Changes 38, 191–208. Quinn, P.F., Beven, K.J., Lamb, R., 1995. The Ln(A/Tan-Beta) index-how to calculate it and how to use it within the topmodel framework. Hydrol. Process. 9 (2), 161–182. Schaake, J., Duan, Q., Koren, V., Hall, A., et al., 2001. Toward improved parameter estimation of land surface hydrology models through the model parameter estimation experiment (MOPEX). In: Dolman, A.J. et al. (Ed.), Soil-vegetation atmosphere transfer schemes and large scale hydrological models. IAHS Publications no. 270, IAHS Press, Oxforshire, UK. Singh, V.P., 2001. Kinematic wave modelling in water resources: a historical perspective. Hydrol. Process. 15, 671–706.

Sorooshian, S., 1981. Parameter estimation of rainfall-runoff models with heteroscedastic streamflow errors: The noninformative data case. J. Hydrol. 52, 127–138. Sorooshian, S., Dracup, J.A, 1980. Stochastic parameter estimation procedures for hydrologic rainfall-runoff models: correlated and heteroscedastic error cases. Water R. Res. 16, 430–442. Sorooshian, S., Duan, Q., Gupta, V.K., 1993. Calibration of rainfallrunoff models: Application of global optimization to the Sacramento soil moisture accounting model. Water R. Res. 29, 1185–1194. Stieglitz, M., Rind, D., Famiglietti, J., Rosenzweig, C., 1997. An Efficient Approach to Modeling the Topographic Control of Surface Hydrology for Regional and Global Climate Modeling. J. Climate 10, 118–137. Tallaksen, L.M., 1995. A review of baseflow recession analysis. J. Hydrol. 165, 349–370. Todini, E., 1995. New trends in modelling soil processes from hillslope to GCM scales. In: Oliver, H.R., Oliver, S.A. (Eds.), The Role of Water and the Hydrological Cycle in Global Change Global Environmental Change, NATO ASI Series, Series I, vol. 31. Springer-Verlag, Ney York, p. 317. Todini, E., 1996. The ARNO rainfall-runoff model. J. Hydrol. 175, 339–382. Todini, E., Dumenill, L., 1999. Estimating large-scale runoff. In: Browning, K.A., Gurney, R.J. (Eds.), Global Energy and Water Cycles. Cambridge University Press, Cambridge, pp. 265–281. Troch, P.A., Detroch, F.P., Brutsaert, W., 1993. Effective watertable depth to describe initial conditions prior to storm rainfall in humid regions. Water R. Res.2, 427–434. Warrach, K., Mengelkamp, H.-T., Rachke, 2001. Raschke, Treatment of frozen soil and snow cover in the land surface model SEWAB. Theor. Appl. Climatol. 69, 23–37. Warrach, K., Stieglitz, M., Mengelkamp, H.T., Raschke, E., 2002. Advantages of a topographically controlled runoff simulation in a Soil-Vegetation-Atmosphere transfer model. J. Hydrometeorol. 3 (4), 131–148. Wilks, D.S., 1995. Statistical methods in the atmospheric sciences. Academic Press, San Diego, CA. van Straten, G., Keesman, K.J., 1991. Uncertainty propagation and speculation in projective forecasts of environmental change: A lake-eutrophication example. J. Forecasting 10, 163–190. Verdin, K.L., and Greenlee, S.K., 1996. Development of continental scale digital elevation models and extraction of hydrographic features. In: Proceedings, Third International Conference/Workshop on Integrating GIS and Environmental Modeling, Santa Fe, New Mexico, 21–26, 1996. National Center for Geographic Information and Analysis, Santa Barbara, California. Vrugt, J.A., Gupta, H.V., Gupta, W., Sorooshian, S., 2003a. A Shuffled Complex Evolution Metropolis algorithm for optimization and uncertaintly assessment of hydrological model parameters. Water R. Res. 39 (8), 1201.

M. Huang, X. Liang / Journal of Hydrology 320 (2006) 37–61 Vrugt, J.A., Gupta, H.V., Gupta, L.A., Bouten, W, 2003b. Effective and efficient algorithm for multiobjective optimization of hydrologic models. Water R. Res. 39 (8), 1214. Yan, J., C, T., 1991. Haan, Multiobjective parameter estimation for hydrologic models—Weighting of errors, Trans. Am. Soc. Agric. Eng. 34, 135–141. Yan, J., C, T., 1991. Haan, Multiobjective parameter estimation for hydrologic models—multiobjective programming, Trans. Am. Soc. Agric. Eng. 34, 848–856.

61

Yapo, P.O., Gupta, H.V., Sorooshian, S., 1996. Calibration of conceptual rainfall-runoff models: Sensitivity to calibration data. J. Hydrol. 181, 23–48. Yapo, P.O., Gupta, H.V., Sorooshian, S., 1998. Multi-objective global optimization for hydrologic models. J. Hydrol. 204 (1–4), 83–97. Zecharias, Y.B., Brutsaert, W., 1988. Recession characteristics of groundwater outflow and base flow from mountainous watersheds. Water R. Res. 24, 1651–1658.