Catchment classification: hydrological analysis of catchment behavior through process-based modeling along a climate gradient

Abstract. Catchment classification is an efficient method to synthesize our understanding of how climate variability and catchment characteristics interact to define hydrological response. One way to accomplish catchment classification is to empirically relate climate and catchment characteristics to hydrologic behavior and to quantify the skill of predicting hydrologic response based on the combination of climate and catchment characteristics. Here we present results using an alternative approach that uses our current level of hydrological understanding, expressed in the form of a process-based model, to interrogate how climate and catchment characteristics interact to produce observed hydrologic response. The model uses topographic, geomorphologic, soil and vegetation information at the catchment scale and conditions parameter values using readily available data on precipitation, temperature and streamflow. It is applicable to a wide range of catchments in different climate settings. We have developed a step-by-step procedure to analyze the observed hydrologic response and to assign parameter values related to specific components of the model. We applied this procedure to 12 catchments across a climate gradient east of the Rocky Mountains, USA. We show that the model is capable of reproducing the observed hydrologic behavior measured through hydrologic signatures chosen at different temporal scales. Next, we analyze the dominant time scales of catchment response and their dimensionless ratios with respect to climate and observable landscape features in an attempt to explain hydrologic partitioning. We find that only a limited number of model parameters can be related to observable landscape features. However, several climate-model time scales, and the associated dimensionless numbers, show scaling relationships with respect to the investigated hydrological signatures (runoff coefficient, baseflow index, and slope of the flow duration curve). Moreover, some dimensionless numbers vary systematically across the climate gradient, possibly as a result of systematic co-variation of climate, vegetation and soil related time scales. If such co-variation can be shown to be robust across many catchments along different climate gradients, it opens perspective for model parameterization in ungauged catchments as well as prediction of hydrologic response in a rapidly changing environment.


Introduction
Catchment classification is an efficient method to synthesize our understanding of how climate variability and catchment characteristics (e.g. vegetation, soils, topography) interact to define hydrological response (McDonnell and Woods, 2004;Wagener et al., 2007). It is also a crucial step in improving predictions in ungauged basins . Differences between the hydrologic responses of catchments can be quantified by means of specific signatures of catchment behavior, such as the runoff coefficient, the flow duration curve or the master recession curve. Gauged catchments can be clustered into separate groups with similar hydrologic signatures and this provides information about similarity of hydrologic responses (Sawicz et al., 2011). Such groups or classes can be regarded as a first step in catchment classification, which offer a catalogue of hydrologic behavior within G. Carrillo et al.: Hydrological analysis of catchment behavior through process-based modeling a region. However, catchment classification is only complete if we understand why certain catchments belong to certain groups of hydrologic behavior, such that we have the means to classify ungauged catchments into their most likely group of behavior.
One way to accomplish catchment classification is to empirically relate climate and catchment characteristics to hydrologic behavior and to quantify the uncertainty of predicting the hydrologic response based on a combination of climate and catchment characteristics. Such a classification system and the related prediction uncertainty will be conditioned by the selection of hydrologic signatures and climate/catchment characteristics, and may result in different classifications depending on the objective of classification (e.g. water balance partitioning, ecological services). In any case, we can call this approach the top-down approach since it is based on measurable hydrologic drivers/responses and landscape features. The measure of uncertainty quantifies the probability of misclassification, and provides insight about how much information is contained in the selected climate and catchment characteristics concerning hydrologic response (Snelder et al., 2005;Oudin et al., 2010). Since there are important surface and subsurface properties that cannot be readily measured or translated into hydrologically relevant information, the uncertainty of classification reflects in part (the lack of) the amount of cross-correlation between observable landscape properties (e.g. vegetation type) and unobservable landscape characteristics (e.g. rooting depth).
An alternative approach, that can partially alleviate the above-mentioned issue of observability, uses our current level of hydrological understanding, expressed in the form of a process-based model, to interrogate how climate and catchment characteristics interact to produce the observed hydrologic response (Sivakummar, 2008). Assuming an appropriate process-based model can be constructed for a wide range of catchments, we can use it to analyze the relationships between hydrologic response and catchment functioning (Samuel et al., 2008). A catchment can be considered as a filter that transforms the climate signal into a hydrologic response by partitioning, storing and releasing incoming energy and water (Black, 1997;Wagener et al., 2007). The different catchment stores (e.g. interception store, root zone store, aquifer store) interact with the different climate fluxes (e.g. rainfall intensity, maximum evapotranspiration) to produce specific time constants of hydrologic behavior (e.g. time to empty root zone store through evapotranspiration). The process-based model can thus be a very useful instrument to analyze different portions of the hydrologic response to identify the important time constants of catchment functioning. For instance, the recession part of a catchment's hydrograph during the dormant season can be used to inform us about the time constant of aquifer release by matching modeled recession flows using lumped aquifer descriptors, such as horizontal hydraulic conductivity or depth to bedrock (Brutsaert and Nieber, 1977;Kirchner, 2009).
Through process-based modeling we can thus obtain estimates of hidden catchment characteristics that are not available in the top-down approach, and ask questions about how these catchment characteristics relate to climate gradients.
Once a sufficient set of catchments across the climatelandscape gradients of a specific region have been analyzed using this bottom-up approach, we can use the model parameters to explain observed hydrologic similarity. Certain model parameters can be prescribed based on observable landscape characteristics (e.g. mean catchment slope, dominant vegetation type). Others cannot be determined a priori and need to be selected during the hydrologic analysis phase. Such hydrologic analysis should not be considered as an automated calibration procedure but rather as a step-bystep methodology to distill relevant information about different catchment functions using appropriate forcing and output variables (Boyle et al., 2000;Yilmaz et al., 2008). The advantage of automated parameter calibration is that it is objective and does not require interaction of the hydrologist with the optimization algorithm . The disadvantage is that typical objective functions used to optimize model performance cannot guarantee that inappropriate combinations of parameter values lead to sets of "behavioral" models (Fenicia et al., 2007), and the functional role of specific parameters is often not preserved .
It is the purpose of this paper to present a general method of hydrologic analysis by means of a process-based model to develop a bottom-up catchment classification system that is compatible with and complementary to top-down classification methods developed elsewhere (Sawicz et al., 2011). In Sect. 2 we present the process-based model to analyze hydrologic response across many catchment in the USA. The model is built around the hillslope-storage Boussinesq (hsB) equation developed by Troch et al. (2003). It uses geomorphologic functions to describe hillslope and channel network topology required to compute subsurface and surface routing. We have chosen this modeling approach because (1) it is parsimonious and thus reduces the problem of equifinality (Beven and Freer, 2001), and (2) it was shown that the hsB equation accurately represents saturated subsurface flow and storage dynamics across complex landscapes . In Sect. 3 we describe a step-by-step procedure to analyze the observed hydrologic response and to assign parameter values related to specific components of the model. It uses different parts of the catchment hydrograph to separate processes in an attempt to reduce parameter uncertainty and to increase the probability to assign a reasonable range of parameter values to different components of the model. In Sect. 4 we apply our hydrologic analysis procedure to 12 catchments selected from the MOPEX (Model Parameter Experiment) database across a climate gradient in the USA, and present a comparison of hydrologic functioning as revealed by our process-based model. In Sects. 5 and 6 we discuss our results and some shortcomings of the bottom-up approach to catchment classification.
2 Process-based model for hydrologic analysis

Modeling principles
The model we developed for the purpose of this study is based on the following principles: (1) the model should be process-based such that we can use it to analyze catchment behavior derived from routine hydro-meteorological observations at the catchment scale, such as daily discharge, temperature and precipitation; (2) the model should be as parsimonious as possible to avoid problems of overparameterization and equifinality (Beven and Freer, 2001;Wagener and Gupta, 2005) and reduce computer processing time; and (3) the model should be applicable to a wide range of catchments across climate and physiographic gradients. In order to represent the dominant functions of a catchment we consider hillslopes and channel network as fundamental hydrologic units . Hillslope land surfaces interact with the atmosphere and partition water and energy fluxes, and drain surface runoff and subsurface flow into the catchment channel network for routing towards the outlet (i.e. point where discharge is measured). Instead of representing individual hillslopes and how they are connected to the channel network, we adopt the modeling approach of Troch et al. (1994) and use the hillslope width function and the channel width function at the catchment scale to represent the geomorphologic structure of the catchment. Each catchment is thus characterized by a hillslope width function (probability density function of water entering the catchment at a given flow distance from the channel network; see also Bogaart and Troch, 2006) and a channel width function (probability density function of surface and subsurface flow entering the channel network at a given flow distance from the outlet) that are derived from available digital elevation models (DEMs). Important additional terrain properties such as average hillslope/channel slope are also estimated from available DEMs. Other landscape properties, such as land use-land cover and soils, available from various spatial databases are further used to assign initial values to process parameters that control the different catchment functions, such as infiltration and interception.

Hillslope and channel routing
The semi-distributed hillslope-storage Boussinesq (hsB) model, developed by Troch et al. (2003), is used to model perched groundwater dynamics at the hillslope spatial scale:  Paniconi et al. (2003) that this model is an adequate and parsimonious representation of three-dimensional saturated subsurface flow along geometrically complex hillslopes. When saturated storage exceeds the local storage capacity S c (= f W h (x)D, where D is maximum perched aquifer depth) the model produces saturation excess overland flow. The partial differential equation is solved numerically for water table dynamics and outflow rate (see Troch et al., 2003 for details). Some fraction of the total percolation from the root zone (see below) is assigned to enter a fractured bedrock aquifer below the perched groundwater table. We assume the outflow from this bedrock aquifer to sustain drought flow at the outlet, and the aquifer dynamics are represented with a lumped non-linear storage model: where q c (x, t) [s −1 ] is specific discharge resulting from a Dirac impulse input at flow distance x upstream, and G. Carrillo et al.: Hydrological analysis of catchment behavior through process-based modeling where N c (x) is the number of channel links at a given flow distance from the catchment outlet and L T is the total channel length. Interpreting the normalized channel width function as the probability density function of receiving lateral inflow at flow distance x from the outlet, the response of the channel network to an instantaneous unit input of water is: with q c (x, t) defined in Eq.
(3). This parsimonious model of channel routing can be used to compute discharge at the catchment outlet given lateral inflows through either infiltration or saturation excess overland flow (assumed to enter the channel network at time of generation). Shallow subsurface flow above a confining soil/bedrock layer draining from the hillslope perched aquifer and deep fractured bedrock baseflow are produced at the catchment outlet and thus do not need to be routed through the channel network (see below).

Root zone water balance
The hillslope perched aquifer interacts with the root zone and exchanges recharge and capillary rise fluxes which depend on root zone moisture content and the depth between the root zone and the local water table h(x, t), called the transmission zone. The root zone water balance is given by: where where p t is throughfall rate and i c is infiltration capacity of the soil. If throughfall rate exceeds the infiltration capacity surface runoff is produced, which is instantaneously added to the lateral flow into the channel network. The throughfall rate is computed as: where ω [m] is canopy storage, ω c is canopy storage capacity and p is precipitation rate. The actual canopy storage is computed using a simple canopy water balance that accounts for precipitation rate and evaporation from the wet canopy and is bounded by [0, ω c ]. The canopy storage capacity is related to the leaf area index (LAI) of the catchment vegetation according to Dickinson (1984): ω c = 0.0002 × LAI. The infiltration capacity of the soil is modeled by means of the time compression approximation suggested by Milly (1986): where k v [m s −1 ] is the vertical hydraulic conductivity, s s [m s −0.5 ] is the soil sorptivity and I c [m] is the cumulative infiltration since start of rain/snow melt event.
The rate of capillary rise is modeled according to Gardner (1958) for steady upward flow from a water table: where ψ c [m] is the depth of the capillary fringe, β c [−] is a reduction factor that varies linearly with θ rz between residual moisture content and saturated moisture content, and a and b are parameters that are related to the Brooks-Corey soil water retention parameters (Eagleson, 1978). Z [m] is the depth (distance) between the bottom of the root zone and the local water table, and thus varies along the hillslope. Percolation or recharge from the bottom of the root zone is assumed to be solely gravity driven and is computed as: where θ r is residual moisture content and θ s is saturated moisture content, and B is the Brooks-Corey pore size distribution index. The transmission zone between root zone and perched aquifer transmits water received from the root zone towards the perched aquifer at a rate defined through Eq. (12) with a transmission zone specific vertical hydraulic conductivity and moisture content. It also transmits capillary rise flux from the perched aquifer to the root zone unaltered, without storage of water. The effective depth of the transmission zone is dynamic and depends on the root zone and perched aquifer storage dynamics (Z decreases as S increases). The difference between the recharge flux from the transmission zone and the capillary rise flux, c r , defines the net recharge, N, to the shallow aquifer.

Land surface energy balance
Evaporation from wet canopy and transpiration from vegetation are estimated by means of the land surface energy budget: ] soil heat flux. Net radiation is estimated from the surface radiation budget accounting for incoming and outgoing shortwave and longwave radiation, depending on surface albedo and emissivity. Since outgoing longwave radiation depends on surface temperature, we solve the energy budget iteratively and assume the surface emissivity constant. The latent heat flux can be approximated as (Brutsaert, 2005): where where T a [K] is air temperature. We solve the land surface energy budget for surface temperature at daily time steps such that we can assume the net ground heat flux to be zero. When the canopy is wet (ω > 0) the canopy resistance is zero. Evaporation from wet canopy is then given by: and ω wc is the areal fraction of wet canopy estimated from Deardorff (1978): The transpiration rate removing moisture from the root zone is given by (Teuling and Troch, 2005): where V RF [−] is the vegetation root fraction, µ [−] is the vegetation light use efficiency, E(r c,min ) [m s −1 ] is the potential vaporization rate using a minimal canopy resistance, β is the transpiration reduction coefficient, given by: with θ w soil moisture content at wilting point and θ c the critical moisture content when transpiration reduction starts.

Snow accumulation and melt
We add a simple snow model for catchments with significant snow days (see below). The snow model accumulates all incoming precipitation in a snow pack when the air temperature is below a certain threshold T m . When air temperature rises above this threshold temperature, the snow melt rate is given by: with M [m s −1 K −1 ] a melt coefficient. The daily melt volume is subsequently removed from the stored snow water equivalent in the snow pack and added to the throughfall.

Model forcing
In this study, we run the model at daily time steps, even though it can be run at shorter time steps (e.g. hourly). Required model forcing are daily precipitation, air temperature, downward short-and longwave radiation, relative humidity, atmospheric pressure and wind speed. Other required model inputs include time evolution of catchment-wide leaf area index (LAI) and albedo. We will discuss the different sources of these input variables in Sect. 4. It should be noted that since we use a semi-distributed version of the hsB-SM model the model forcing data is basin-averaged, and soil and vegetation type are effectively uniform, as in Woods (2003). This no doubt will add to modeling uncertainty but is unavoidable in order to keep the number of model parameters to a minimum.

Characteristic time scales and dimensionless numbers
The different components of the process-based model, in combination with catchment-scale climate forcing, reveal characteristic time scales of hydrologic response that are related to catchment hydrologic functions of partitioning, storage, and release of water. Therefore, such characteristic time scales are important indicators of catchment behavior and can help to relate above and below ground landscape characteristics to water balance dynamics. They can also be combined to form dimensionless numbers that can be related to hydrologic regimes through empirical or analytically derived scaling relations (Berne et al., 2005;Harman and Sivapalan, 2009).

Canopy time scales
The time scale associated with filling up the canopy interception storage capacity, ω c , is given by: G. Carrillo et al.: Hydrological analysis of catchment behavior through process-based modeling where p is the average rainfall intensity when it rains. The time scale associated with emptying the interception storage is given by: with e wc the average wet canopy evaporation. Obviously, both the average rainfall intensity and average wet canopy evaporation vary throughout the year, such that the seasonal canopy time scales can be either larger or smaller than the annual averages defined above. In any case, the interception storage capacity is at most a few mm such that in most climates the canopy time scales are of the order of a few days at maximum, and typically less than one day. The time scales are also of same order of magnitude and thus their ratio, reflecting the competition between filling and emptying the interception storage, is close to 1.

Snow pack time scales
The characteristic time scale of snowmelt can be defined as: where s is the average maximum snow accumulation, and Q m is the average snow melt rate during snow melt season. This time scale is important to define what type of runoff generation mechanism is likely to dominate (saturation excess vs. shallow subsurface flow) during snow melt by comparing it with characteristic time scales of root zone and perched aquifer processes (see below).

Root zone time scales
The time scale related to filling the root zone storage by rainfall is defined as: where θ is the average soil moisture content of the root zone and p t is the average throughfall rate when T a > T m . Similarly, the time scale related to filling the root zone by snow melt is given by: It is possible to specify different average soil moisture contents during the rainy season and the snow melt season to reflect different wetness conditions, if necessary. Time scales related to emptying the root zone storage in the absence of capillary rise are: where θ FC is soil moisture content at field capacity, r is the average recharge rate and t is the average transpiration rate. Different combinations of these time scales express competition between different processes affecting the water balance dynamics. For instance, the ratio of the latter two reveals the competition in the catchment between baseflow generation and vegetation water use.

Transmission zone time scales
As mentioned earlier, the depth of the transmission zone is time variable as it depends on the soil moisture dynamics in the root zone as well as on storage dynamics in the perched aquifer. Nevertheless, an average transmission zone storage capacity can be numerically derived from the model simulations and used to define the following time scales of transmission zone filling and emptying: In Eq. (28), Z is average transmission zone depth, θ s is saturated moisture content of the transmission zone, θ is average moisture content and r and r t are average recharge rate from root zone and transmission zone, respectively.

Perched aquifer time scales
Much work has been done on defining characteristic time scales of shallow aquifer dynamics (Brutsaert, 1994;Troch et al., 2004;Berne et al., 2005;Harman and Sivapalan, 2009). The characteristic time scale of advection-driven (kinematic) flow in perched aquifers is given by (Berne et al., 2005;Harman and Sivapalan, 2009): where L is hillslope length (maximum flow distance between divide and nearest channel), pD is average saturated thickness, and a c is the rate of con/divergence of the hillslope width function. Likewise, the characteristic time scale of diffusion-driven flow is given by: Their ratio, τ K /τ U , defines the hillslope Péclet number (Pe; Berne et al., 2005) and high values of Pe indicate that shallow subsurface flow is mainly dominated by gravity drainage. Harman and Sivapalan (2009) extended the similarity framework of Berne et al. (2005) to account for the responsiveness of the hillslope subsurface flow to temporal variability of the recharge events, as well as for the effects of lower boundary condition of hillslope drainage. They used the concept of hydrologic regimes of Robinson and Sivapalan (1997) to develop a hillslope subsurface flow classification system based on the Pe number and the dimensionless characteristic time of recharge events: where τ r is the average storm duration and τ hc is the concentration time of the hillslope. Either the advection or the diffusion time scale defined above can be used to estimate the hillslope concentration time. Their classification system defines slow/fast, advection/diffusion dominated subsurface flow, depending on the numerical value of Pe (below 1: diffusion; above 1: advection) and π r (below 1: slow; above 1: fast, although the separation between fast and slow flow in the diffusion dominated case depends on the boundary condition assumed: fixed (small) flow depth vs. kinematic).

Fractured bedrock time scales
Time scales for non-linear reservoirs representing baseflow dynamics have been proposed by Woods (2003). In many cases, the master baseflow recession curve of a given catchment converges to a straight sloping line in semi-logarithmic plots of ln(Q b ) versus time, indicating that most deep aquifer dynamics are best represented by a linear reservoir equation with b = 1. In that case, the characteristic time scale of deep (fractured bedrock) aquifer dynamics is given by 1/a, the reservoir time constant.

Channel network time scales
The advective characteristic time scale of channel flow is given by: where L c is flow length along the channel network from the centroid to the outlet and V is average flow velocity. Obviously, the channel flow Froude number is an appropriate dimensionless number to characterize the flow regime.

Linking parameter values to dominant process behavior
The above-described hydrologic model is one of many alternative process-based models that can be formulated to describe different surface and subsurface stores and their interactions that generate streamflow (Jothityangkoon and Clark et al., 2008). Within the context of such models, routine hydro-meteorological observations can be analyzed to inform us about the different catchment functions of partitioning, storage and release of incoming water and energy fluxes. During different parts of the hydrologic response not all components of the model are equally active, such that one can link parameter values to specific storage dynamics to avoid unwanted parameter interactions often encountered in automatic calibration procedures. In the following we describe a step-by-step procedure of linking model parameters to specific hydrologic responses generated by the proposed model. This procedure can easily be modified when other process-based models are used.

Dormant vs. growing season
First, we divide the hydrologic year into two periods: one when the vegetation is dormant and one when the vegetation is active (growing season). This decision is based on analyzing the average leaf area index (LAI) curve derived from several years of remote sensing observations at the catchment scale. In this study we use MODIS (Moderate Resolution Imaging Spectroradiometer; http://modis-land. gsfc.nasa.gov/lai.htm) data and more specifically the LAI product available at https://lpdaac.usgs.gov/lpdaac/products/ modis products table from 2000 to 2008. From the annual signals of LAI the average LAI curve is derived and subsequently rescaled using the minimum and maximum average LAI. The hydrologic year is then separated into the dormant season and growing season using the time instances when the rescaled LAI curve crosses the 50 % cut-off level (Fig. 1). This method is similar to the phenology model for monitoring vegetation responses developed by White et al. (1997), and seems to be able to capture the inflexion points of the average LAI curve well.

Step 1: baseflow recession and aquifer dynamics
An obvious starting point for hydrologic analysis of catchment response is when the catchment is non-driven and relaxes from previous hydro-meteorological fluxes that have replenished some/all stores. In order to isolate several possible release fluxes from the catchment it is best to start focusing on baseflow recessions during the dormant season. Such recession hydrographs will be minimally affected by root water uptake and subsequent transpiration losses and thus can be considered mainly controlled by aquifer properties. Our process-based model considers two separate aquifer stores: the near-surface perched aquifer that develops during wet period above a confining layer (i.e. fractured bedrock with reduced vertical hydraulic conductivity), and a deep aquifer that receives a fraction of all percolation water from the root zone (i.e. a fractured bedrock aquifer). To relate baseflow recessions to these aquifer stores we perform a baseflow separation as follows:    with Q(t) total streamflow at time t, Q b the computed baseflow contribution to total streamflow (Q b ≤ Q), and ε a low-pass filter parameter (Arnold and Allen, 1999;Eckhardt, 2005). The filter parameter ε is set for all catchments at 0.925. Since the purpose of the study is to address hydrologic similarity across a climate gradient, the selection of a different cut-off level would not change the relative differences between the catchments (a desired characteristic of the data manipulation), but obviously will affect to some degree the absolute values. Next, all recession periods during the dormant season are selected for recession curve analysis (Fig. 2). The catchment master recession curve (MRC) is constructed by time shifting individual recession curves to match the lower end of the baseflow values, and progresses from low to high baseflow values. This procedure is described in more detail in Posavec et al. (2006). Subsequently, the MRC is defined as the smoothed lower envelope of all observed recession curves. According to our conceptual model of baseflow generation, we can consider the early part of the MRC as being composed of both perched and deep aquifer contributions while the late part of the MRC is solely composed of deep aquifer contributions. Therefore, starting from the low flow end of the MRC, the deep aquifer parameters are estimated to match that part of the MRC. In all applications of the model to our study sites (see Sect. 4) we have observed that the lower end of the MRC can be approximated by means of a linear reservoir model, characterized by a time constant of storage release given by the reciprocal value of the slope of the linear regression line through the lower end of the MRC (Fig. 2). Parameter values are estimated using the downhill Simplex method (Nelder and Mead, 1965) with least square error objective function. The inset of Fig. 2 shows a Brutsaert-Nieber plot of recession rates versus baseflow of binned observations and MRC. The lower end reveals the linear reservoir response of the deep aquifer whereas the upper end shows the non-linear recession characteristics against which the hillslope-storage Boussinesq equation is calibrated.
Using the deep aquifer model we can now identify the perched aquifer contributions to the early part of the MRC. Once isolated from the deep aquifer contributions, the perched aquifer recession curve is used to estimate the parameters controlling release from the hillslope-storage Boussinesq model (viz. horizontal hydraulic conductivity, k h , and drainable porosity, f ). The maximum perched aquifer baseflow contribution is used to define the steady-state recharge rate required to generate this amount of drainage. This recharge rate is then applied to the hsB model to bring it to steady-state, after which recharge is set to zero and the model parameters are estimated such that the time history of relaxation from the maximum baseflow matches the observed recession. Since these parameters also define the total storage during steady state, this procedure is repeated until no further improvements, measured by means of least square error, are obtained using the downhill Simplex parameter estimation algorithm (Nelder and Mead, 1965). The maximum water table depth during steady state is next used to define the upper boundary of perched aquifer storage capacity, expressed as maximum perched aquifer depth, D.
Other conceptualizations of observed baseflow dynamics could have been proposed to capture the early-time nonlinear behavior, such as the transmissivity feedback mechanism (Bishop, 1991). Given the size of the selected catchments and the lack of biogeochemical data it is very difficult to unambiguously decide which subsurface flow mechanism is responsible for the observed baseflow dynamics and both conceptualizations (the one used in this study and the one based on transmissivity feedback) are equally likely.

Step 2: streamflow generation during dormant season
The total amount of baseflow produced by our model does not depend on the parameters assigned during the previous step, but on the total amount of infiltrated water that percolates down to the perched water table and the deep aquifer. Likewise, total streamflow generated by our model during the dormant season will include direct runoff produced either through infiltration excess or saturation excess. The next step therefore is to assign values to parameters controlling the infiltration and percolation processes in the root zone. From available soil databases, such as STATSGO and SURGGO, we select the dominant soil type within a given catchment. From this soil type we assign values of total porosity and residual porosity, θ s and θ r , using look-up tables from Clapp and Hornberger (1978). Other soil hydraulic parameters, viz. sorptivity and vertical hydraulic conductivity, are estimated by means of the downhill Simplex algorithm using a multiobjective function that accounts for the absolute values of normalized residuals between modeled and observed baseflow, direct runoff, and total streamflow volumes. In this way we select infiltration and percolation parameters that match all runoff generation mechanisms active in the catchment during the dormant season. Parameters that control root water uptake are set to typical values from look-up tables associated with dominant vegetation type. Once reasonable parameter values for the hydraulic properties of the root zone soil are obtained, other critical processes such as deep aquifer percolation and snow melt, are added to the list of parameters to be optimized. The fraction of total percolation that enters the deep aquifer will control late time recession dynamics. Snowmelt during the dormant season may or may not be an active process, depending on the climate of the basin. In any case, we test whether better modeling performance can be achieved by adding these three parameters (fraction of total percolation rate, melt rate M, and threshold temperature T m ). Since we use basin-average and daily averaged temperature to force the snow melt model, the value of the temperature threshold and melt rate should be interpreted with care.

Step 3: streamflow generation during growing season
During the growing season, parameters that control root water uptake and vegetation transpiration will have an important effect on hydrological partitioning of incoming water and energy fluxes. These parameters include soil and vegetation parameters such as critical moisture content, θ c , wilting point moisture content, θ w , vegetation root fraction, V RF , vegetation light use efficiency, µ, as well as aerodynamic parameters, such as zero plane displacement height, d, and roughness length, z 0 . These aerodynamic parameters are related to the vegetation height through (Brutsaert, 2005): and therefore vegetation height, H v , is used during the parameter estimation procedure. The five parameters are estimated using the same procedure as described above (downhill Simplex). Once reasonable parameter values are obtained, the snowmelt parameters are revisited to investigate if better model performance can be obtained by means of modified values from previous iterations.

Step 4: channel network routing
The next step takes the daily-generated surface runoff (both infiltration excess and saturation excess) and uses Eq. (6) to route these volumes to the catchment outlet. These routed volumes are added to the daily subsurface flow from the perched aquifer and fractured bedrock aquifer. The two routing parameters, c and F , are estimated by maximizing the  Nash-Sutcliffe efficiency (Nash and Sutcliffe, 1970) measure for streamflow values.

Matching hydrologic signatures
The final step in our model identification procedure is to compare modeled and observed hydrologic signatures, such as the annual runoff coefficient, annual baseflow index and the slope of the flow duration curve Yilmaz et al., 2008). The annual runoff coefficient for any given hydrologic year is defined as: where t is day in hydrologic year (1 October-30 September). Similarly, the annual baseflow index is defined as: The slope of the flow duration curve is defined as (Yadav et al., 2007;Sawicz et al., 2011): where Q 33% and Q 66% are the flow values exceeded 33 % and 66 % of the time, respectively. Discrepancies between modeled and observed hydrologic signatures are used to repeat the parameter estimation procedure after Step 1 until no further improvements in reproducing these signatures are obtained.

Study sites across climate gradient
We applied the above described hydrologic analysis procedure to 12 MOPEX catchment east of the Rocky Mountains, USA. These catchments were previously used in van Werkhoven et al. (2008) to study SAC-SMA (Sacramento Soil Moisture Accounting) model parameter sensitivities across a hydroclimate gradient using multiple time periods between 1980-1989. As can be seen from the listed wetness indices and runoff coefficients in Fig. 3, these catchments represent a wide range of climate and hydrologic regimes. Table 1 lists some catchment characteristics of our 12 study sites. Catchment area ranges from 1000 km 2 to 4500 km 2 . Mean catchment elevation ranges from about 100 to 800 m a.s.l. The mean annual precipitation ranges from 750 mm to 1500 mm, and the mean annual potential evapotranspiration ranges from 1500 mm to 700 mm.

Forcing data
The model uses the following eight variables as input time series: precipitation, land surface albedo, air temperature, long and short wave downward radiation, atmospheric pressure, actual vapor pressure and wind speed. Daily precipitation data is provided through the MOPEX website (ftp://hydrology.nws.noaa.gov/pub/gcip/mopex/US Data/) (Duan et al., 2006). The other seven variables are derived from the 3-h, 1/8 degree hydroclimate data set developed by Maurer et al. (2002), and available at http://   www.hydro.washington.edu. The 3-h data are converted to daily averages and then spatially averaged over the catchments using a weighted averaging procedure that accounts for complete or partial coverage of data grid and catchment boundaries.

A priori parameter assignments
For each basin, the MOPEX database provides fractional spatial coverage of each of the 16 USDA soil types, as well as the fractional spatial coverage of vegetation type according to the University of Maryland vegetation classification system (see also http://www.geog.umd.edu/landcover/ global-cover.html). From this information, the dominant soil type and vegetation type is selected and typical parameter values are selected from Clapp and Hornberger (1978) for total soil porosity, and from the North American Land Data Assimilation System -NLDAS (http://ldas.gsfc.nasa.gov/nldas/ NLDASmapveg.php) database for initial values of root zone depth and vegetation height.

Modeling results
Figure 4 compares observed and modeled average runoff coefficient for the period 1990-1999 for all 12 catchments. We used 1990-1994 to calibrate the model and ran the calibrated model for 1990-1999. As can be seen, the model has captured very well the average annual water balance, and similar results were obtained for the inter-annual variability of hydrologic partitioning (not shown). From Fig. 5 we can see that the model also captured very well the fraction of total streamflow that is generated as baseflow. The observed baseflow indices in Fig. 5 are computed after baseflow separation, as described in Sect. 3, while the modeled baseflow indices are computed from the generated baseflow volumes from the perched and deep aquifer in the model. There is a slight tendency to underestimate the baseflow contribution to streamflow but the differences between observed and modeled average baseflow index are not statistically significant.  In order to evaluate the model performance at daily time steps, Fig. 6 shows the observed and modeled flow duration curves for the period 1990-1999 for all catchments. Even though the model efficiency to reproduce observed hydrographs is moderate (see inset values of Nash-Sutcliffe efficiencies in Fig. 6), the match with observed flow duration curves is remarkable at all flow levels (with a few exceptions). This suggests that the model captures the dynamic transformation of climate forcing into streamflow rather well but that timing of individual storm events may not be modeled accurately. For the purpose of this study we consider it more important to be able to reproduce the different modes of response (in terms of frequencies of low, medium and high flow) given certain climate forcing than to match/overparameterize the model to fit hydrographs. Figure 7 compares the monthly regime curves of precipitation, evapotranspiration and discharge for two catchments in different climate settings. San Marcos catchment in Texas (left panel of Fig. 7) is a water-limited catchment, whereas Amite catchment in Louisiana (right panel of Fig. 7) is a more energy-limited catchment. The model reproduces the discharge regime curve for both catchments remarkably well, illustrating that the model is capable of filtering different climate signals in ways that are comparable with the real catchment filters. Similar results were obtained for the other 10 catchments (not shown).

Model parameters and time scales
Table 2 lists all model parameters for all 12 catchments, together with catchment characteristics derived from available geographic information, such as drainage area, mean catchment slope and mean channel slope. Total porosity was selected from look-up tables (Clapp and Hornberger, 1978) based on dominant soil type. All other parameters were obtained using the methods described in Sect. 3. From these model parameters we have computed the different time scales discussed in Sect. 2.4 (see Table 3). Many different dimensionless numbers can now be formulated as ratios of time scales. In the next section we relate these time scales and dimensionless numbers to hydrologic signatures to reveal scaling relationships that could be used to determine hydrologic similarity between different catchments.
An attempt to perform an automated parameter sensitivity analysis failed due to the highly coupled and non-linear character of the model equations, which caused instabilities in the numerical solution of Eq. (1). In future work we will reformulate the presented model and replace the dynamic groundwater equation with derived storage-discharge relationships. This will remove most of the issues of numerical stability and will allow testing of the parameter uncertainty and sensitivity   Table 2 are.

Regionalization and scaling relationships
We regressed all readily available catchment characteristics, such as drainage area and mean catchment slope, to the different model parameters, in an attempt to reveal regionalization patterns. Not many linear regressions between catchment characteristics and model parameters were statistically significant at 95 % confidence limits. Table 4 shows all regression relationships that were significant with p < 0.05 (some were significant at p < 0.01, indicated by * ). Figure 8 shows some of these statistically significant relationships for the no-snow dominated catchments. Only very few significant relationships showed up for all 12 catchments or for the 6 snow dominated catchments, indicating that the parameters of the snow dominated catchments were not related to catchment characteristics and therefore could not be regionalized. The remaining regression relations for the 6 no-snow catchments appear to be rather strong. In particular, information of minimum and maximum LAI can be translated to root zone and vegetation parameters in the model quite reliably. Mean elevation of the 6 catchments seems to be strongly related to the saturated hydraulic conductivity of the transmission zone,  while catchment slope defines vegetation height. The latter relationship is most likely caused by the a priori choice of vegetation height from land cover databases that show a similar relation between vegetation height and catchment slope (K. Sawicz, personal communication, 2010). However, these regressions are mainly the result of our hydrograph analysis to inform model parameters rather than from regression catchment characteristics and hydrologic response. Obviously more work is needed to define the robustness of these relationships, their physical meaning (why is hydraulic conductivity of the transmission zone related to mean elevation?), as well as answering the question why no significant relationships showed up for the snow dominated catchments. It is possible that the soil parameters defined during the dormant season are affected during the snow accumulation period while in fact most partitioning processes controlled by these soil parameters are inactive. Another possible explanation is that different parameters of the model control different runoff generation mechanisms in different climates (van Werkhoven et al., 2008). We will address this issue in our on-going research. Next, we regress the model time scales to hydrologic signatures (runoff coefficient, baseflow index and slope of the flow duration curve; Table 5). Figure 9 shows some of these regression relationships. Even though our initial regression analysis is based on simple linear regression, some  regression functions were altered when non-linear relations were apparent from the data trends. Again, no strong regression relationships showed up for the snow-dominated catchments, so Fig. 9 shows only significant relationships for no snow catchments. It should be kept in mind that the reported R 2 values apply to the initial linear regression analysis (as well as the significance levels), even though after inspection of the data trends it was clear that non-linear (power law) regressions better represent the patterns. The runoff coefficient is related to 3 time scales of the models: the time scale related to emptying the canopy store (function of potential evaporation, and thus strongly related to climate), the time scale associated with emptying the root zone by transpiration (a function of actual evapotranspiration and thus part of the competition between ET and drainage), the time scale related to emptying the transmission zone through drainage (clearly defining the generation of slow flow at the expense of ET). RC is clearly also affected by the mean storm duration, a climate time scale and not a model time scale. Zooming in on the latter two model time scales and how they relate to RC, it is interesting to note that RC increases when it takes longer for the root zone to be emptied by ET, indicating that in such situations water can move through the root zone to become baseflow or is more likely to generate quick flow. Likewise, RC is higher when it takes less time to empty the transmission zone through drainage, indicating the (lack  of) feedback between subsurface moisture storage and ET through capillary fluxes. Baseflow index is strongly and linearly related to the time scale associated with filling up the root zone by rainfall. When that time scale is short, less water will leave the catchment as baseflow and more as surface runoff. This does not seem to be affected much by climate, since even in our most arid catchments the baseflow index can be high (e.g. GUA which is a karst dominated catchment).
The only relationship that holds for all 12 catchments is between the slope of the FDC and the time scale of the deep aquifer. This relationship is especially robust because of the physical link between short time scales of the linear reservoir behavior and the release of water from the catchment, as captured by the slope of the FDC.
Finally, we regressed hydrologic signatures with different dimensionless numbers, created as ratios of the model time scales. We used π i-j to indicate the dimensionless number created by time scale i over time scale j , with i and j indices referring to time scales listed in Table 3. For instance, π 2-8 is defined by the ratio of the time scale to empty canopy storage by potential ET and the time scale of the perched aquifer advection. It is clear from Fig. 10 that the aridity index is a strong control on the runoff coefficient, for all 12 catchments (but somehow more strong for the no-snow catchments). This is no surprise since the left panels of Fig. 10 are nothing but the Budyko curve for our catchments. However, for the no-snow catchments, a dimensionless number defined by the time scale to empty the canopy storage (linked to PET) and the diffusion time scale of perched aquifer drainage (linked to catchment early-stage drainage), does equally well to explain the observed variance compared to the aridity index (R 2 of 0.935 at p < 0.01 vs. 0.926 at p < 0.05, respectively). Figure 11 suggests that the same dimensionless number (π 15-12 : ratio of interstorm duration to deep aquifer time constant) explains both the baseflow index and the slope of the FDC for all 12 catchments, even though stronger relationships are possible when separating snow from no-snow catchments and when using different dimensionless numbers. In any case, the time scales of aquifer drainage always play an important role to explain these hydrologic signatures, indicating that subsurface catchment characteristics are controlling release of water stored in the saturated zone.

Co-variation of climate, vegetation and soil time scales
Finally, we investigated how different ratios of time scales are related to each other. If significant (linear or non-linear) relationships exist between different dimensionless numbers characterizing the catchments, this could indicate that different time scales interact to create systematic emerging patterns of hydrologic partitioning across the climate gradient. For example, Fig. 12 shows a strong linear relationship (R 2 = 0.86; p < 0.0001) between π 14-9 and π 2-8 . Catchments with low values of these two dimensionless numbers have climates characterized by high PET and short storm durations, have low canopy storage capacity (low LAI), slowly drain the root zone and transmission zone, and have low perched aquifer storage capacity and hence high transmission zone storage capacity. The opposite is true for catchments with high values for these 2 dimensionless numbers. The data suggest either a linear trend between those two extremes or a non-linear (sigmoid) trend. The latter has a lower mean absolute error regarding the data points. All this indicates that climate, vegetation and subsurface characteristics of root zone, transmission zone and perched aquifer somehow co-evolve along the climate gradient.
To test whether the observed trend can be explained solely by trends in mean storm duration (time scale 14 in y-axis) and mean potential evapotranspiration (denominator of time scale 2 in x-axis), we plotted these two climate variables against each other in the bottom panel of Fig. 12. It is clear that the trend of the top panel cannot be explained by climate only, and that we have to take vegetation and subsurface characteristics identified through our process-based hydrograph analysis into account. This reinforces the notion that climate, vegetation and subsurface storage and release properties of these basins co-vary systematically across the climate gradient.   Fig. 11. Significant linear relationships between baseflow index and slope of the FDC and different dimensionless numbers related to model time scales (4: root zone filling by rain; 9: transmission zone emptying; 11: Boussinesq aquifer diffusion; 12: deep aquifer; 15: mean inter-storm duration). Triangles indicate no snow catchments and dots represent snow catchments.
Apparently, our hydrograph analysis with the aid of the process-based model resulted in a systematic variation of subsurface properties, expressed as time scales related to the root zone and transmission zone, between dry less vegetated and wet more vegetated catchments. At the far left in Fig. 12 appear the catchments situated in Texas (GUA and SAN) and Missouri (SPR) and at the far right are catchments situated in West Virginia (POT and TYG). If such relationship between climate, vegetation and soil time scales can be shown to hold for other catchments along similar climate gradients, it can provide guidance for catchment model parameterization that would apply to ungauged basins. Obviously more research is required to support this conclusion.

Limitation of bottom-up modeling approach to explain hydrologic similarity
There are a number of disadvantages associated with the procedure outlined in this paper. First, model construction is to some degree subjective and different hydrologists will develop different generic catchment models with the same purpose of capturing hydrologic response. Therefore, model time scales derived from individual model components are not universal and will depend on the model construction.
Model inter-comparison is needed to check to what degree different model formulations will lead to different conclusions about the cause of hydrologic similarity. The observed scaling relations between model time scales and hydrologic signatures presented in this study should therefore be interpreted with care, as they are probably unique to the modeling procedure used in this study, and more work is needed to test their robustness. Figure 12: Top: relationship between two dimensionless numbers characterizing co-variation between climate forcing, canopy storage and belowground storage and release characteristics of the 12 catchments. The data suggest either a linear or a sigmoid functional relationship, with MAE (mean absolute error) smallest for the latter; Bottom: Lack of significant relationship between mean storm duration and potential evapotranspiration illustrating that trend in top panel is not due to climate gradient only.
58 Figure 12: Top: relationship between two dimensionless numbers characterizing co-variation between climate forcing, canopy storage and belowground storage and release characteristics of the 12 catchments. The data suggest either a linear or a sigmoid functional relationship, with MAE (mean absolute error) smallest for the latter; Bottom: Lack of significant relationship between mean storm duration and potential evapotranspiration illustrating that trend in top panel is not due to climate gradient only.

Fig. 12.
Left panel: relationship between two dimensionless numbers characterizing co-variation between climate forcing, canopy storage and belowground storage and release characteristics of the 12 catchments. The data suggest either a linear or a sigmoid functional relationship, with MAE (mean absolute error) smallest for the latter; right panel: lack of significant relationship between mean storm duration and potential evapotranspiration illustrating that trend in top panel is not due to climate gradient only.
Second, the described hydrograph analysis to select appropriate model parameters is time consuming. Some of the hydrograph analysis can be performed with the aid of computer scripts, but still requires supervision of a skilled hydrologist. Applying our procedure to all 280 catchments used by Sawicz et al. (2011) is therefore beyond the scope of this work, but will ultimately be required to test the robustness of our preliminary results.
Third, our method requires daily observations of precipitation, temperature, streamflow, and other hydrometeorological variables. These are, by definition, not all available in ungauged basins. Even though several model parameters can be selected a priori from available databases and remote sensing products, it is unclear whether this can lead to the construction of behavioral models that can guide catchment classification methods in ungauged basins. However, the observed co-variation between model time scales along a climate gradient is an encouraging result for application in ungauged basins.

Conclusions
In this study, we developed a parsimonious process-based modeling procedure to investigate hydrologic similarity across catchments. The basic idea behind this approach is that we use the model to interrogate hydrologic behavior manifested in streamflow dynamics that are the result of how catchment properties, such as soils, aquifers, geomorphology and vegetation filter available water and energy fluxes. Different parts of the hydrograph reflect different catchment functions (e.g. baseflow recession during dormant season) that can be captured in individual model components through parameter selection informed by careful hydrograph analysis. The resulting parameter values reveal characteristic model time scales of partitioning, storage and release of water at the catchment scale. These model time scales can be grouped as dimensionless numbers that serve as similarity indices to explain specific hydrologic behavior.
We applied this procedure to 12 catchments across a climate gradient in the eastern US. The process-based model is capable of representing accurately observed hydrologic responses at annual, seasonal and daily time scales. Some model parameters are related to specific catchment properties, which offer potential for regionalization. At the same time, we show that inter-catchment variability of three hydrologic signatures (runoff coefficient, baseflow index and slope of the flow duration curve) can be explained by variability in model time scales and their dimensionless ratios.
Perhaps the most intriguing result of our study is shown in Fig. 12. Figure 12 suggests that climate, vegetation and soil storage and conductivity co-vary predictably across a climate gradient. Apparently, available energy and storm characteristics interact with catchment properties, such as vegetation cover and belowground water storage and release capacity, and result in specific water balance partitioning. It is well known that local vegetation and soil properties vary systematically along climate gradients in similar geologic settings (Rasmussen et al., 2011;Anderson and Goulden, 2011). It stands to reason that co-evolution of climate, vegetation and soils is also present at larger scales, and that such co-evolution of catchment properties manifest itself in how catchments partition incoming water and energy fluxes. Obviously, at regional scales the initial conditions set by geology and tectonics can strongly control evolutionary trajectories and can result in complicated patterns that are difficult to unravel. Our preliminary results suggest that such co-evolution of catchment properties can be revealed through process-based model interrogation of observed hydrologic behavior, confirming a similar experience in Australian catchments (Jothityangkoon and , which highlights the diagnostic role that process models can be expected to perform in the future. If we go further and develop ways to understand how such co-evolution came about, how it is manifested in hydrologic response, and how it is affected by geologic and tectonic processes, we can make important progress in our ability to predict hydrologic response in ungauged basins as well as in our ability to predict how hydrologic systems will evolve in a changing environment (Wagener et al., 2010).