Controls on runoff generation and scale-dependence in a distributed hydrologic model

<jats:p>Abstract. Hydrologic response in natural catchments is controlled by a set of complex interactions between storm properties, basin characteristics and antecedent wetness conditions. This study investigates the transient runoff response to spatially-uniform storms of varying properties using a distributed model of the coupled surface-subsurface system, which treats heterogeneities in topography, soils and vegetation. We demonstrate the control that the partitioning into multiple runoff mechanisms (infiltration-excess, saturation-excess, perched return flow and groundwater exfiltration) has on nonlinearities in the rainfall-runoff transformation and its scale-dependence. Antecedent wetness imposed through a distributed water table position is varied to illustrate its effect on runoff generation. Results indicate that transitions observed in basin flood response (magnitude, timing and volume) can be explained by shifts in the surface-subsurface partitioning. An analysis of the spatial organization of runoff production also shows that multiple mechanisms have specific catchment niches and can occur simultaneously in the basin. In addition, catchment scale plays an important role in the distribution of runoff production as basin characteristics (soils, vegetation, topography and initial wetness) are varied with basin area. For example, we illustrate how storm characteristics and antecedent wetness play a dramatic role in the scaling properties of the catchment runoff ratio.                         </jats:p>

section_content: The response of natural catchments to precipitation depends on the mechanisms of runoff generation and their spatial and temporal distribution. Nevertheless, quantitative descriptions of distributed runoff generation are difficult to acquire in field settings. To address this, a number of distributed hydrologic models have been developed to represent multiple runoff mechanisms and their variability in a catchment (e. g., Smith and Hebbert, 1983; Abbott et al., 1986; Ogden and Julien, 1993; Paniconi and Wood, 1993; Wigmosta et al., 1994; Garrote and Bras, 1995; VanderKwaak and Loague, 2001; Ivanov et al., 2004a; Rigon et al., 2006). To date, however, limited attention has been placed on model analyses that identify the spatiotemporal variability of runoff generation and its underlying causes. Few studies have quantified runoff partitioning into the surface and subsurface mechanisms that arise from interactions of a distributed groundwater table and variably-saturated catchment surface. Understanding runoff production in these complex settings is critical for identifying the mechanisms for flood generation, and its features (e. g., magnitude, timing, volume). More importantly, flood characteristics and their scale-dependence may be linked to the runoff mechanisms excited as rainfall properties interact with variable basin characteristics and pre-storm wetness conditions. 

In this study, we investigate the transient basin response to storm events using a coupled surface-subsurface model that accounts for heterogeneities in catchment properties, prestorm hydrologic states and rainfall forcing. Our analysis is based on the physical processes represented in the model and how these interact to generate runoff. Our objective is to E. R. Vivoni et al. : Runoff generation and scale-dependence identify nonlinearity and scale effects in the simulated flood response and their relation to the underlying runoff mechanisms simulated in the basin. The nonlinearity in catchment response refers to the observation that a unit increase in rainfall may not produce an equivalent increase in measured runoff. Numerous studies have recognized that catchment runoff response can be highly nonlinear (e. g., Betson, 1964; Caroni et al., 1986; Troch et al., 1993b; Risbey and Entekhabi, 1996). Nonlinearity is potentially related to the runoff mechanisms operating in a basin, which have been shown to affect flood frequency (Sivapalan et al., 1990; Robinson and Sivapalan, 1997). A possible cause of the observed runoff nonlinearity is the difference in response times between surface and subsurface flows. Surface runoff is a quick response as infiltration is significantly reduced due to soil properties or saturated conditions, while subsurface runoff can be much slower as subsurface flow paths may substantially delay travel toward the channel network. Previous studies have identified that transitions in runoff mechanisms can be controlled by rainfall properties, such as storm duration and intensity (e. g., Larsen et al., 1994; Menabde and Sivapalan, 2001). Nevertheless, antecedent wetness conditions imposed by the water table position after long interstorm periods should also have an important effect on transitions between surface and subsurface mechanisms. 

In addition to influencing runoff nonlinearity, surface and subsurface mechanisms may impart a signature on the runoff scale-dependence, which refers to the effect of catchment area (or size) on runoff properties. Numerous studies have shown that discharge may exhibit stronger or weaker dependence on scale as area increases, depending on the hydrologic processes in the basin (e. g., Goodrich et al., 1997; Morrison and Smith, 2001). Both antecedent wetness and storm properties, through forcing particular runoff mechanisms, may also create differences in flood characteristics at different scales. The scaling properties of interest in the runoff response include the flood peak magnitude, time to peak and volume. In small catchments, the short response times relative to storm duration can lead to contributions from surface and subsurface runoff to the overall basin response. As catchment scale increases, storm duration is shorter than the response time from each mechanism, thus leading to timevarying partial contributions from various mechanisms (e. g., Sivapalan et al., 2002). Thus, flood characteristics are anticipated to result from the interrelationship of runoff processes with multiple time constants and basin size. An important question is whether the catchment runoff ratio exhibits scaledependence and if this is linked to the runoff mechanisms excited at particular basin scales. 

In this study, we utilize the "TIN-based Real-time Integrated Basin Simulator" (tRIBS) (Ivanov et al., 2004a) with spatially-uniform forcing to identify the surface and subsurface runoff mechanisms occurring in a complex basin. The numerical model is used as an interpretive tool for addressing the control exerted by storm properties and antecedent wetness on the distributed basin response. In the experiments, variations in rainfall duration and intensity are used to mimic characteristic storms, while initial wetness conditions are altered by imposing a distributed water table at different baseflow levels. The distribution of initial soil water content in the unsaturated zone corresponds to hydraulic equilibrium and explicitly depends on the initialized water table position. Within the model construct, we seek to identify if exciting a range of different runoff mechanisms with inherent variations in response time, magnitude and spatial organization can explain runoff nonlinearity. Furthermore, we attempt to understand how basin heterogeneities interact with storm characteristics to create favorable sites for runoff production. If persistent catchment niches exist for runoff generation, then the scale-dependence of flood properties may be related to how runoff partitioning changes with basin area. As basins grow in size, the heterogeneous mix of catchment characteristics varies, possibly leading to predictable patterns in runoff production. Understanding the effects of runoff partitioning on catchment response at various scales is important to advance our capability to predict behavior in ungauged basins.
section_name: Introduction

section_content: tRIBS is a continuous, physically-based, fully-distributed model designed for hydrologic research and forecasting (Ivanov et al., 2004a). The model provides a spatiallyexplicit treatment of basin heterogeneities in topography, soils and aquifer properties, vegetation and atmospheric forcing. As in grid models, runoff production and routing are tracked over complex terrain captured by individual nodes. In tRIBS, Voronoi polygons (control volume associated with a given TIN node) uniquely associated with a triangulated irregular network (TIN) are used as a finite-volume domain for mass balance and flux computations (Tucker et al., 2001). In the following, we present a brief discussion of the model physics, emphasizing those components most relevant to understanding the coupled surface-subsurface basin response. The reader is referred to Ivanov et al. (2004a) for additional details and discussions of model limitations.
section_name: Coupled surface-subsurface distributed model

section_content: A catchment is represented in tRIBS through a TIN consisting of elevation, channel, and basin boundary nodes (Vivoni et al., 2004). Triangulated irregular networks are a piecewise linear interpolation of a set of points, sampled from a digital elevation model (DEM), resulting in triangular facets of varying size. The triangulation represents topographically complex surfaces that include hillslopes, valleys, floodplains and ridges. The stream network is composed of a set of channels ranging from headwater tributaries to large, meandering rivers. The channel cross section is established through geomorphic relations to contributing area (Ivanov et al., 2004b). The soil profile and shallow aquifer are bounded by a spatially distributed bedrock assumed to be an impermeable surface. A shallow aquifer interacts with the stream network and land surface to produce saturated areas that expand and contract (e. g., de Vries, 1995; Lamb et al., 2000).
section_name: Model domain representation

section_content: Basin hydrologic response requires an appropriate depiction of the two-way interaction between surface and subsurface processes. The model accounts for moving infiltration fronts, water table fluctuations and moisture losses due to evapotranspiration and groundwater drainage. Each element consists of a sloped column of heterogeneous, anisotropic soil with an exponential decrease in saturated hydraulic conductivity (e. g., Beven and Kirkby, 1979; Beven, 1982) :

where K Si (z) is the saturated hydraulic conductivity at depth z in the normal or parallel directions (i=n or p), K oi is the saturated hydraulic conductivity at the soil surface (z=0), and f is a hydraulic conductivity decay parameter. A kinematic approximation for unsaturated flow is used to compute infiltration and propagate moisture fronts in the soil column (Cabral et al., 1992; Garrote and Bras, 1995; Ivanov, 2002). 

The unsaturated moisture profile is determined from hydrostatic equilibrium using the Brook and Corey (1964) parameterization as:

where θ (z) is the soil moisture at depth z, θ r and θ s are the residual and saturation soil moisture contents, N wt is the depth to the local water table, b is the air entry bubbling pressure and λ o is the pore-size distribution index (Ivanov et al., 2004a). Coupled to the vertical dynamics is lateral moisture redistribution in the vadose zone and shallow aquifer driven by gradients in surface and groundwater topography. In the unsaturated zone, horizontal flow between contiguous elements is computed over the saturated wedge and along the steepest direction. In the shallow aquifer, a quasi-three-dimensional model based on the Dupuit-Forchheimer approximation redistributes groundwater from recharge zones to discharge areas. Lateral exchanges between elements are controlled by hydraulic gradient as:

where Q S is the groundwater outflux, w is the flow width, tan β w is the local water table slope and T is the depth averaged aquifer transmissivity:

where D is the bedrock depth and a r is the anisotropy ratio (K op /K on ). Water table dynamics are computed from groundwater fluxes, vertical recharge and exfiltration. Overall, the water table position anchors the soil moisture profile and determines regions of saturation prior to a storm.
section_name: Coupled unsaturated and saturated dynamics

section_content: The coupled nature of the unsaturated and saturated processes results in a robust set of runoff mechanisms. Four basic runoff types are simulated in the tRIBS model: infiltration-excess runoff (R I ) (Horton, 1933), saturationexcess runoff (R S ) (Dunne and Black, 1970), groundwater exfiltration (R G ) (Hursh and Brater, 1941), and perched return flow (R P ) (Weyman, 1970). Total runoff (R) is composed of the four production mechanisms:

where R I +R S andR P +R G are the surface and subsurface components. Infiltration-and saturation-excess runoff are rapid surface responses as infiltration is limited by soil conditions, while perched return flow and groundwater exfiltration are slower mechanisms as subsurface flow delays the response to rainfall. Over complex terrain, the occurrence and frequency of runoff generation depends on the spatiotemporal characteristics of catchment topography, soils, climate, rainfall and antecedent wetness. Given this variability, it is recognized that watershed response can correspond to runoff production from multiple mechanisms arranged in spatially distinct areas or possibly due to a single dominant type in the basin (e. g., Freeze, 1974; Dunne, 1978; Smith and Hebbert, 1983). Runoff production from multiple mechanisms will vary with the rainfall and landscape factors influencing the coupled unsaturated-saturated dynamics.
section_name: Runoff generation processes

section_content: The soil moisture state resulting from the interaction of infiltration, runoff and subsurface flows is coupled to losses from evaporation and plant transpiration. The surface energy balance, R n -G=λE+H, is solved as function of surface temperature using parameterizations for net radiation (R n ), latent (λE), sensible (H ) and ground (G) heat fluxes. Total evapotranspiration (ET) is determined from moist bare soil, intercepted water and plant transpiration based on soil and vegetation parameters that include vegetative cover (v), surface albedo (a), canopy height (h), stomatal resistance (r s ) and an optical coefficient (K t ), in addition to atmospheric conditions (e. g., air temperature, relative humidity, pressure, wind speed) and solar radiation. Moisture in the top surface layer, root zone and canopy storage play a key role in limiting ET when the atmospheric demand is high (Ivanov et al., 2004a). Conversely, evapotranspiration impacts the unsaturated-saturated processes and leads to differences in runoff production.
section_name: Surface energy balance and evapotranspiration

section_content: Runoff generated at each element is routed across an individual hillslope overland flow path and then through the channel network. The hillslope paths are defined over the edges of the triangular facets that connect a node to the closest downstream stream node (Tucker et al., 2001). A nonlinear relation is used to determine velocity over a hillslope path (Ivanov et al., 2004a) :

where v h is the hillslope velocity, A h is the upslope contributing area, Q is the discharge at the downstream channel node, and r and c v are spatially-uniform parameters of the velocity relation. Thus, overland travel time (t h =l h /v h ) is a function of discharge (Q) and hillslope path length (l h ). 

Overland flow from multiple hillslope nodes serves as lateral inflow into a kinematic wave, one-dimensional routing scheme solved in the channel network (Ivanov et al., 2004a). Channel travel time (t c =v c / l c ) depends on the channel link distance (l c ) and the discharge (Q=v c A c ) through each link. 

For a wide, rectangular channel (A c =bH), discharge for each link is:

where n is the Manning coefficient, S is the channel slope, b is the channel width, and H is the water depth. As overland travel time is faster than groundwater pathways, the partitioning of precipitation into surface and subsurface flow is critical for determining the basin response.
section_name: Hillslope and channel flow routing

section_content: We investigate the effects of storm properties and antecedent catchment wetness on basin response through the use of tRIBS in a complex, humid watershed in northeastern Oklahoma. The basin is well-suited for investigating the relation between runoff mechanisms and flood response due to a spatially-variable groundwater table that reacts quickly to rainfall (Imes and Emmett, 1994; Sloan, 2000). Furthermore, the Springfield Plateau aquifer provides baseflow discharge throughout the year. Large flood occurrences are typically due to mid-latitude frontal storms during early fall and late spring (Bradley and Smith, 1994; Michaud et al., 2001).
section_name: Hydrometeorological observations and catchment simulations

section_content: tRIBS is applied to the Baron Fork at Eldon, OK (BF, 808 km 2 ) that includes two gauged subbasins, the Peacheater Creek at Christie, OK (PC, 65 km 2 ) and the Baron Fork at Dutch Mills, AR (DM, 107 km 2 ). Figure 1a presents the basins, gauging stations, and stream network overlaid on a TIN derived from a USGS 30-m DEM. Parts of the basin are rugged and heavily dissected while others are flat or gently sloping. The basin is composed of a mixture of forest (52. 2%), croplands (46. 3%), and towns (1. 3%) (Fig. 1b ), Table 1. Catchment runoff ratio ( =R/P, dimensionless) and specific discharge (q s =Q/A, in mm) for selected events (labeled 1, 2 and 3 in Fig. 2 ) during April-May 1999 for the three stream gauges in the Baron Fork watershed. Basin-averaged rainfall depths in the Baron Fork are 34. 05 mm, 35. 11 mm and 60. 51 mm for events 1, 2, and 3, respectively. BF DM PC (808 km 2 ) (107 km 2 ) (65 km 2 ) Runoff Ratio ( ) Event 1: 2-5 April 1999 0 while surface soil texture is silt loam (94%) and fine sandy loam (6%). The channel network has a maximum length of 67. 3 km and a mean drainage density of 0. 86 km -1. Channel geometry is parameterized using geomorphic relations from Carpenter et al. (2001). In addition to the gauging sites, twelve (12) ungauged interior basins were delineated to represent a range of catchment area, A=0. 78 to 808 km 2 (Fig. 1b ). Several studies have focused on the Baron Fork due to its unregulated nature, high stream gauge density and long time series of radar rainfall data (e. g., Johnson et al., 1999; Smith et al., 2004a).
section_name: Study catchments

section_content: Hydrologic measurements in the catchment consist of three USGS gauges (BF, PC, DM), overlapping NEXRAD radars and one meteorological station (Fig. 1b ). In the basin, gaugecorrected radar is a reliable source of hourly (4-km by 4-km) precipitation data (e. g., Smith et al., 1996; Young et al., 2000; Grassotti et al., 2003). The mean annual rainfall of 1240 mm is distributed in two wet periods (March to June; September to November). Over the years 1993 to 2000, the mean areal rainfall in the Baron Fork varied from light rain (e. g., less than 1 mm/h for one hour) to intense storms (e. g., greater than 20 mm/h over 6 h). The variability in rainfall duration, intensity and spatial distribution, in addition to prior wetness in the basin, leads to a complex runoff response during storms (Smith et al., 2004b). Multiple runoff mechanisms occur due to the heterogeneity in basin properties and lead to flood hydrographs of varying magnitudes (Finnerty et al., 1997; Carpenter et al., 2001; Ivanov et al., 2004b). Large flood events in the basin have a tendency to occur in early fall and late spring due to frontal storms (Bradley and Smith, 1994). As an example, Fig. 2 shows rainfall and runoff observations for April-May 1999. For each indicated event, Table 1 presents the runoff ratio ( =R/P, where R and P are the event runoff and rainfall volumes). The observations depicted here provide an indication of the runoff magnitudes and the effects of catchment scale on the basin response. For example, the maximum flood for 1999 occurred on 4-5 May (event 3) with discharges of 351, 206 and 12 m 3 /s in the BF, DM and PC basins. The high flood magnitude resulted from wet antecedent conditions, as evidenced by the high baseflow prior to the event. Pre-storm wetness is also reflected in high runoff ratios ( =0. 41 to 0. 68). In contrast, a storm of similar magnitude in 22-24 April (event 2) led to a weaker runoff response and lower runoff ratio ( =0. 06 to 0. 17). Prior to this event, the antecedent wetness and baseflow in the basin were low due to a long interstorm period. Similarly, the runoff response during 2-5 April (event 1) was weaker than the annual flood. Nevertheless, the runoff ratio for each basin ( =0. 08 to 0. 21), was higher than for event 2 despite having similar rainfall volumes and a lower rainfall peak. As illustrated in this example, the interactions between storm properties and antecedent wetness are expected to partially explain the observed nonlinearity in flood response with respect to rainfall. Systematic variations in the runoff ratio among the nested basins also reveal scale-dependence due to internal differences in basin runoff dynamics. These limited observations, however, do not allow rigorous study of the runoff nonlinearity and scale-dependence as a function of storm and initial wetness conditions. In this study, we address this observed behavior through numerical experiments designed to identify the effects of storm conditions and antecedent wetness on the simulated rainfall-runoff processes.
section_name: Rainfall and streamflow observations

section_content: Hydrologic model calibration was accomplished by adjusting soil, vegetation and routing parameters to match the observed hydrograph at the basin outlet over the 1993-2000 period. The multi-step calibration procedure, fully described in Ivanov et al. (2004b), ensures that the model performance is reliable at the element, hillslope and catchment scales. Initial parameter estimates were based on physical relationships to soils and vegetation types (e. g., Rawls et al., 1982; Bras, 1990). Calibration is focused on a limited number of parameters to which the model is most sensitive within narrow, physically plausible ranges. Table 2 presents the set of calibrated model parameters derived from Ivanov et al. (2004b). Spatial parameter variability is captured by the use of soil and vegetation classes where within-class parameter variations are not allowed. For this study, as in Ivanov et al. (2004b), the spatial variability of soils and vegetation are overlapping and correspond to forest, grassland and urban classifications. In this manner, the potential for over-parameterization is reduced and the overall strengths of the distributed approach (e. g., capturing spatial variability) are highlighted. Ivanov et al. (2004a, b) present an analysis of the model performance in the Baron Fork during a long-term simulation in terms of discharge at the gauging stations and the distribution of hydrologic states in the basin. NEXRAD rainfall data from Smith et al. (2004a) and weather data from Maurer et al. (2002) were used to force the distributed model over separate calibration and verification periods (see Smith et al., 2004a). Results indicate that model simulations capture the flood response at the gauging stations and reproduce the observed nonlinearities in the rainfall-runoff dynamics via the multiple runoff mechanisms. Spatial distributions of runoff production also illustrated the topographic, soils and veg-etation controls on basin response, although data to verify the predicted patterns are unavailable for the area. Overall, Ivanov et al. (2004b) obtained a robust parameter set for the basin through the multi-year application with radar rainfall forcing.
section_name: Hydrologic model calibration and verification

section_content: The spatiotemporal variability in catchment response and its relation to runoff generation mechanisms is explored in this study through a set of numerical simulations. Given the calibrated model parameter set, the sensitivity of the basin runoff response is assessed as a function of the initial distributed water table position and storm event properties (duration, t r, and intensity, i). For the experiments, a set of discrete, spatially-uniform storm events are used as model forcing as our objective is to identify nonlinearity and scale effects in the flood response and its relation to underlying runoff mechanisms. Our analysis considers the total runoff response consisting of both the baseflow and stormflow contributions. Selected storm durations and intensities represent conditions leading to flooding in the Southern Great Plains (Michaud et al., 2001). We focus on flood events in the fall period to minimize effects from snowmelt that may induce antecedent wetness not related to the water table position. In addition, the impact of evapotranspiration is reduced during the autumn due to tree leaf-fall and low atmospheric demand. Moisture conditions in the basin during the fall period are a strong function of the interannual fluctuations in summer rainfall. We capture this interannual variability through the different initial water table states corresponding to pre-storm baseflow levels. In this manner, rainfall forcing in the fall is applied to the expected baseflow conditions given possible amounts of summer precipitation.
section_name: Numerical simulations of catchment response

section_content: Numerical simulations of basin response to uniform rainfall are discussed with respect to the variations in the antecedent wetness and storm properties. Model predictions of hydrologic response are related to the underlying runoff mechanisms and to catchment scale. Both spatially-averaged and distributed metrics are used to quantify surface-subsurface interactions and their impact on the basin response. Prior to describing the response to storm forcing, we present the antecedent conditions imposed by a set of initial groundwater distributions.
section_name: Results and discussion

section_content: A key to partitioning runoff into surface and subsurface components in complex terrain is the position of the water table depth and its control on surface soil moisture state (e. g., Troch et al., 1993a; Salvucci and Entekhabi, 1995) To initialize the model, both the saturated zone thickness and soil moisture profile must be specified. As shown in Eq. ( 2 ), the initial moisture profile θ(z) is determined entirely by the soil hydraulic properties and the depth to water table (N wt ). Ivanov et al. (2004a) describe a method for constructing an initial groundwater distribution based on the topographic or wetness index following Sivapalan et al. (1987). As an alternative, we determine a model-based initial water table position through a basin drainage experiment as in Vivoni et al. (2005). By allowing a fully-saturated catchment to drain for a long period of time (∼10 years), the transient readjustment of the subsurface head field occurs in the context of the basin geomorphology following subsurface flow pathways and gradients. In the absence of rainfall and evapotranspiration, the discharge (or baseflow) produced in the catchment is due exclusively to drainage from the saturated zone. Baseflow drainage is governed by the calibrated model parameters as well as by the geometry of the hillslope and channel system. VanderKwaak and Loague (2001) used a similar strategy to initialize a surface-subsurface model prior to a storm event simulation. 

Figure 3 presents the basin response to the long-term drainage experiment in the form of the groundwater rating curve, which relates subsurface discharge (Q b ) to the water table position (N wt ) (Eltahir and Yeh, 1999). The depth to the water table is captured by the first and second spatial moments of the N wt field (mean, µ and standard deviation, σ ). Note the nonlinearity in the Q b -µ [N wt ] relation arises due to the interaction of the groundwater level with the surface topography and aquifer thickness (Marani et al., 2001). As an aquifer initially drains, sharp baseflow decreases occur as the saturated zone is rapidly depleted along seepage faces. As drainage continues, Q b declines at a slower rate in response to a deeper N wt with limited surface interactions. Interestingly, the spatial variability in the water table σ [N wt ] exhibits an inflection point as the water table deepens (Q b ∼1 m 3 /s). This inflection is an indication of the water table disconnecting from the surface along the stream network, resulting in a reduction of σ [N wt ]. In light of this, the catchment response to rainfall should be a strong function of the nonlinear groundwater rating curve. To capture this nonlinearity, we selected three initial states from the observed discharge record: wet (Q b =5 m 3 /s), medium (Q b =2 m 3 /s) and dry (Q b =0. 5 m 3 /s). These baseflow levels represent prestorm basin conditions with exceedence probabilities of 0. 48 (wet), 0. 75 (medium), and 0. 99 (dry) derived from USGS gauging records. 

The spatial organization of the groundwater table recession is further explored in Fig. 4 as the time evolution of the water table drop ( N i wt ) at different basin locations. Model elements are classified according to an index of hydrologic similarity:

where A c is the upslope area and tan β is the terrain slope. Large λ occur for flat regions along channels, while small λ depict high slope upland areas. λ is similar to the topographic index of Beven and Kirkby (1979) total drainage period ( N T wt ) are strongly related to basin location. Larger decreases in N wt are observed for regions with λ=9 to 18, suggesting these areas dominate baseflow production. Total drainage, however, is composed of separate periods with different spatial contributions. From the initial to the wet state ( N 1 wt ), water table decreases are larger in steep upland regions (λ=7. 5 to 15). This pattern continues for N 2 wt suggesting that upland areas are initially responsible for baseflow. Flat lowland regions (λ=20 to 27) do not exhibit water table decreases until N 3 wt, indicating these areas contribute to baseflow during drier states. Note the shift in the peak of the spatial distribution of N wt to larger λ values as the water table becomes progressively more disconnected from the land surface. 

Given the spatial heterogeneity in the groundwater recession, it is possible that catchment scale influences the antecedent wetness conditions. Figure 5 presents the scaledependence of the mean depth to groundwater, µ [N wt ], for the three initial states (wet, medium, dry). Catchment scale (A) variation is captured by sampling fifteen subbasins ranging in area from 0. 78 to 808 km 2 (see Vivoni et al., 2006 becomes progressively shallower (closer to the land surface, z=0) as more low lying areas near the stream network are sampled. Since lowland regions have less effective drainage in natural settings, inclusion of these areas in the basin average reduces µ [N wt ]. Furthermore, the spatial variability of the depth to groundwater, σ [N wt ], increases with A (not shown), suggesting that heterogeneity in pre-storm wetness increases with catchment scale. Thus, the initial condition in larger basins appears to be wetter than in smaller basins, on average, but exhibits a higher spatial variation. While a limited sampling of the internal basins is performed, the results indicate the potential scale-dependence in the initial groundwater conditions. As storage capacity depends on N wt, the spatial distribution of the antecedent wetness should play an important role in the catchment runoff response.
section_name: Role of antecedent groundwater conditions in nonlinearity and scale-dependence

section_content: The variation of catchment runoff response with storm forcing and antecedent wetness is explored by altering the spatially-uniform rainfall properties (duration, t r, and intensity, i) and the distributed groundwater table position. Rainfall properties are varied over a range of intensities (ten values from i=1 to 40 mm/h) for three durations (t r =1, 6, 12 h) to obtain thirty (i,t r ) storm pairs for each initial groundwater condition (90 total runs). The selected storm properties mimic the rainfall attributes of autumn events in the region with respect to storm duration and intensity (Bradley and Smith, 1994; Michaud et al., 2001). For each simulation, the uniform and stationary storm event arrives 120 h into the simulation period which lasts for 15 days, sufficiently long to capture the flood recession in the basin. As indicated previously, the initial water table distributions vary significantly in the mean groundwater depth: µ [N wt ]=-4. 57 m (wet), -5. 49 m (medium) and -6. 79 m (dry). In terms of the pre-storm surface saturation, however, the three initial conditions vary only slightly, with saturated surface area fractions (A s /A) of 0. 06 (wet), 0. 03 (medium) and 0. 01 (dry), where A s is the basin area with saturated conditions (θ=θ s ) in the top 10 cm. 

Figure 6 presents the simulated flood hydrographs for two selected combinations of storm duration and intensity at multiple basins in the Baron Fork. The first storm type (i=1 mm/h, t r =12 h) corresponds to a long-duration, lowintensity drizzle, while the second storm type (i=40 mm/h, t r =1 h) is more reflective of a short-duration, high-intensity thunderstorm. Storm types represent fall frontal storms (drizzle) and squall lines (thunderstorm) in the area (Grassotti et al., 2003; Van Horne et al., 2006). Differences in the storm properties impact the flood response as rainfall interacts in complex ways with the coupled surface-subsurface system. For example, the peak magnitude, time to peak and recession behavior vary dramatically for the two storms. Note, however, the storms do not have the same rainfall volume and are selected simply to show the range of basin responses. Antecedent wetness also has important effects on the flood discharge. Under dry conditions, the basin storage capacity leads to low runoff generation; whereas under wet conditions, the shallower water table promotes a more intense flood response. Differences in flood response between the dry and wet conditions are not constant within the storms, suggesting that interplay between rainfall and initial wetness is responsible for runoff nonlinearity. 

Figure 7 further explores the variation of the peak discharge (q p ) and the time to peak (t p ) with storm properties and antecedent wetness. The changes in q p and t p represent the full range of simulated conditions. With respect to storm intensity, two hydrologic regimes can be identified at low and high rainfall intensities, with a transition near i=10 to 20 mm/h. Note the variation in the rate of increase in q p and in the rate of decrease in t p with increasing rainfall intensity (i) in the two regimes. This variation is indicative of a shift in the underlying runoff mechanisms. The transition is potentially due to the saturated hydraulic conductivity (K s ) and its control on runoff generation. For the Baron Fork, the basin-averaged K s =19. 7 mm/h agrees well with the transition in hydrologic regime. For interior basins, the transition may occur at values near the areal-averaged K s which varies with subbasin properties. Interestingly, the two regimes respond differently to changes in the antecedent wetness. At low i values, the variation of q p and t p with µ [N wt ] is more pronounced than at high rainfall rates for all storm durations. (a, c, e ) depicts the variation in peak discharge (q p in m 3 /s) with rainfall intensity (i) and duration (t r ) for the wet, medium and dry conditions. The bottom row (b, d, f) presents the variation in the time to peak discharge (t p in h) for the same conditions. For each case, the hydrograph characteristics are presented for the Baron Fork outlet (808 km 2 ). For example, the ratio of q p (wet) to q p (dry) varies from 14. 2 at i=1 mm/h, to 1. 6 at i=40 mm/h for t r =1 h. This suggests that the initial water table position plays a more critical role for the low intensity rainfall regime, whose upper limit is set to some extent by the basin-averaged infiltration characteristics. 

The effects of storm properties and antecedent wetness on runoff nonlinearity are shown in Fig. 8 through use of the runoff ratio ( =R/P ). For a linear system, is a con-stant, while for a nonlinear response, depends on rainfall amounts (Risbey and Entekhabi, 1996). Clearly, the runoff response exhibits a large variation of with storm properties and is modulated by pre-storm wetness. Interestingly, three general types of nonlinearity are exhibited: (1) asymptotically increasing (i) (e. g., dry, t r =12 h); (2) asymptotically decreasing (i) (e. g., wet, t r =1 h); and (3) a (i) function with a minimum value (e. g., medium, t r =1 h). For each type, a transition is present between the low and high rainfall intensity regimes. Higher nonlinearity is observed for low i across all wetness states as the runoff production is limited. In general, increasing rainfall intensity leads to an asymptotic rise in, which under certain circumstances can exhibit linearity (e. g., constant for wet, t r =1 h, i=20 to 40 mm/h). The minima observed in the (i) relation at various wetness conditions suggests a transition from subsurface to surface regimes that may be associated with the relative importance of pre-storm and flood contributions. A closer look at the various runoff mechanisms can reveal their potential linkage to the runoff nonlinearity.
section_name: Effect of storm properties and antecedent wetness on runoff response and nonlinearity

section_content: Multiple runoff mechanisms arise within a particular basin location due to the interaction between infiltration fronts, the water table position and lateral moisture transport (e. g., Ivanov et al., 2004a). The superposition of different runoff mechanisms generated at individual catchment sites leads to the flood hydrograph in the channel network. Thus, the effects of storm properties and antecedent wetness conditions on runoff nonlinearity should be mediated by transitions in runoff generation. To illustrate this, Fig. 9 presents the total runoff response decomposed into infiltration-excess (R I ), saturation-excess (R S ), perched return flow (R P ) and groundwater exfiltration (R G ) for the drizzle and thunderstorm events. Total discharge from the drizzle event is dominated by R S during the flood, with relatively high R G throughout the period. In contrast, the thunderstorm event is composed primarily of R I with minor contributions from R S and R G. The different runoff mechanisms result in a ten-fold variation in the runoff ratio: =0. 033 (drizzle) and =0. 336 (thunderstorm). These two examples clearly indicate the potential for process controls on runoff response and nonlinearity, as explored further in the following. 

Figure 10 presents a more direct comparison of the runoff ratio and the individual runoff generation mechanisms to in-vestigate the process controls on nonlinearity. To facilitate the comparison, the runoff ratio ( ) and runoff components (R I, R S, R P and R G ) are normalized as:

where max is the maximum value of (i) over the rainfall intensity interval (i=1 to 40 mm/h), shown as connected circular symbols, and:

depicted as shaded regions corresponding to the fraction of the total runoff volume (R). As noted previously, the total surface runoff consists of F I +F S (dark shaded area), while the subsurface component is F P +F G (light shaded region). 

Clearly, a strong relation exists between the runoff nonlinearity and the surface-subsurface partitioning, as the shape of the * (i) function is tied to changes in the shaded regions. For instance, * (i) decreases with i when subsurface runoff (F G +F P ) is the primary mechanism (e. g., wet, t r =1, 6, 12 h). The transition in runoff nonlinearity, characterized by a minima in * (i), occurs when surface runoff exceeds F I +F S >0. 25. Note the position of the * (i) minima varies along i according to the runoff fractions. For cases with high surface runoff, * (i) increases with i (e. g., dry, t r =1, 6, 12 h). This comparison suggests that surface-subsurface interactions control the degree and functional form of the nonlinearity. 

It is interesting to consider the spatial distribution of runoff generation to further explore its link with runoff nonlinearity. As an example, Fig. 11 shows the time-averaged runoff rate from each component as a function of the topographic index (λ) for the drizzle and thunderstorm events. Runoff rate distributions indicate the basin locations that produce large runoff peaks via a particular mechanism (Ivanov et al., 2004b). The drizzle exhibits a response dominated by R S at high λ locations in flat, near channel areas. As wetness increases, saturated areas contributing to R S expand towards lower λ and produce higher runoff. In addition, R G is observed at high λ for the wet condition, indicating subsurface stormflow. For the thunderstorm event, R I is the primary runoff mechanism and occurs at a range of intermediate locations, λ=10 to 22. Interestingly, R I decreases for wetter conditions and shifts toward lower λ, indicating that the shallower water table positions decrease infiltration capacity. At the same time, smaller contributions from R S occur at high λ and expand toward low λ with increased wetness. These examples show that storm properties and antecedent wetness interact with the basin to simultaneously generate runoff mechanisms in specific catchment niches. As the distribution of niches varies with scale, process controls on runoff nonlinearity are analyzed with regard to basin area in the next section.
section_name: Process controls on runoff response, nonlinearity and spatial distribution

section_content: The variation of runoff with catchment scale captures the complex interactions between storm properties and internal hydrologic processes. Runoff statistics, such as the peak discharge of a specific return period, have been shown to follow scaling laws with basin area (e. g., Gupta et al., 1994; Menabde et al., 2001; Ogden and Dawdy, 2003). For individual events, the flood hydrograph peak and time to peak may follow scaling relations of the form:

where c q, c t, α and β are parameters. Studies from numerous basins have found Eq. ( 11 ) to hold, with 0. 5<α<1, depending on return period (Gupta and Dawdy, 1995; Blöschl and Sivapalan, 1997). In addition, Robinson and Sivapalan (1997) derived a response time scaling similar to Eq. ( 12 ) with β=0. 5. Figure 12 shows the variation of q p and t p with A for a limited set of subbasins (0. 78 to 808 km 2 ). Note that a power law relation is apparent in q p and t p for the two storm types. Scaling exponents vary with wetness and storm properties in the range of 0. 74<α<0. 88 and 0. 34<β<0. 86. An increase in initial wetness leads to a higher α, indicating that scale-dependence is stronger as the basin becomes wetter. Initial wetness also has a larger effect on the drizzle event, composed primarily of R S and R G, suggesting the scaling relations may be related to the runoff mechanisms. 

As storm properties are varied, changes are anticipated in the scaling parameters. To further explore the runoff scale-dependence, Fig. 13 shows the runoff ratio ( ) as a function of catchment area (A). Despite the variability in with A, scale-dependence is observed in the runoff ratio in the limited set of subbasins. For the drizzle event, decreases with A, with a stronger scale-dependence in the wet condition. In this case, small basins have high subsurface fractions (F G +F P ) (not shown), thus leading to high runoff ratios ( ∼0. 8). As A increases, decreases as the relative amount of subsurface runoff diminishes. This suggests that larger basins are less sensitive to the initial water table position for drizzle events. Consistent with this, the dry condition, where F G +F P are small, has a runoff ratio which does not vary with A ( ∼0. 04), indicating that the impact of the initial condition has been effectively reduced across all basin scales. In contrast, the thunderstorm event exhibits opposing runoff ratio scaling behavior. For the three wetness conditions, is generally increasing with A, with greater scale-dependence for the dry condition. As this event is dominated by surface runoff (F I +F S ), variations of with A depend on internal surface properties. Small forested basins have a lower due the low F I +F S relative to the rainfall (not shown). As A increases, surface runoff increases in basins with large fractions of non-forested area, generally leading to higher. Note that a maximum value in is observed at A=20 to 60 km 2, depending on the initial wetness. The peak identifies basins with surface properties that promote higher runoff (e. g., lower forest fractions). Interestingly, a slight decrease in occurs for the largest basins, which results from lower runoff production as more permeable regions are sampled. Since this occurs for all wetness conditions in the thunderstorm event, it is likely due to the fraction of forested areas rather than initial wetness. Overall, this preliminary evidence suggests that runoff mechanisms are responsible for the variation in runoff ratio with catchment scale and that storm forcing and antecedent wetness play important roles. 

Figure 14 presents a more in-depth analysis of the relationship of the runoff ratio and the underlying runoff mechanisms at three catchment scales (A=0. 78, 65. 06 and 808. 39 km 2 ) and for the three initial conditions (wet, medium and dry). 

Here, the normalized runoff ratio ( * = / max ) is presented as a function of the peak discharge (q p ) resulting from each storm condition (e. g., thirty pairs of i and t r ). Peak discharge is used as a surrogate for event intensity such that * (q p ) reflects runoff nonlinearity (e. g., variable * (q p ) implies nonlinearity). The normalized runoff ratio is compared to the underlying mechanisms represented by the space-time averaged runoff fractions ( I, S, P, G ) calculated over the full q p range in each basin. Time-averaged runoff production at a basin location (r) is computed over the simulation duration (T d ) as:

where r j is the runoff rate at the j th time step. The spacetime averaged runoff production, r, is estimated from the time-averaged runoff rates at all basin locations (total of N elements) as:

where a k is the area of the kth element and A is the total basin area. For each mechanism (x=I, S, P and G), the space-time averaged runoff fraction ( x ) is estimated as:

where r T is the total space-time averaged runoff rate. Surface runoff consists of I + S (dark shaded area), while the subsurface component is P + G (light shaded region). 

A strong relation exists between the * (q p ) form and the runoff partitioning into surface or subsurface components. Note the existence of two regimes in each catchment: (1) decreasing * (q p ) for low values of q p, where total runoff production is primarily subsurface flow ( P + G ); and (2) increasing * (q p ) for high values of q p, where surface runoff ( I + S ) is the principal mechanism. The minima in the * (q p ) relation indicates a transition from subsurface to surface runoff as event intensity (q p ) increases. This transition occurs at different of q p for each basin suggesting that scale-dependence influences process controls on runoff nonlinearity. For example, the smallest basin exhibits a minima in * (q p ) at low q p (∼10 -2 -10 -1 m 3 /s) as the fraction of subsurface runoff decreases sharply for larger events. Furthermore, the dry condition exhibits an increasing * (q p ) throughout the q p range due to low subsurface runoff fractions. For the larger basins, the minima in * (q p ) occurs at higher q p (10 -1 -10 m 3 /s) when surface runoff increases relative to subsurface contributions. While shifts in the runoff fractions are not smooth due to the effects of varying storm forcing along the q p range, there is a general increase in surface runoff for more intense events (e. g., from I + S =∼0. 2 to ∼0. 8 for the largest basin). 

It is also interesting to note how runoff fractions change with catchment scale and initial wetness. At one extreme, the small basin under the dry condition has primarily surface runoff ( I + S ), indicating a relatively deep groundwater with negligible subsurface runoff. At the other extreme, the largest basin under the wet condition has a mixture of runoff mechanisms, implying a strong surface-subsurface interaction. Under all initial conditions, each basin exhibits a nonlinear runoff response as * (q p ) varies with q p, except possibly for cases where surface runoff is overwhelmingly dominant (e. g., small basin, wet condition, q p >0. 1 m 3 /s). Runoff production variations with catchment scale and initial wetness are linked to the interaction between the water table and land surface. As basin area increases, a higher fraction of low-lying locations exist in the catchment which can interact with the shallow water table. As basin wetness increases, the water table is closer to the land surface, leading to higher amounts of subsurface contribution. Thus, increases in either catchment scale or initial wetness promote the existence of a diverse set of runoff production mechanisms.
section_name: Catchment scale-dependence of runoff response, nonlinearity and process controls

section_content: The coupled surface-subsurface model utilized in this study has afforded the opportunity to examine the links between hydrograph characteristics and the underlying mechanisms leading to runoff production. In particular, the observed nonlinearity in the rainfall-runoff transformation has been tied directly to the partitioning of a storm event of particular intensity and duration into surface and subsurface contributions as simulated in the distributed model. Runoff mechanisms can vary in their spatiotemporal distribution as a function of storm properties and initial wetness and may occur simultaneous and preferentially in different catchment locations. An analysis of runoff production in a limited set of subbasins also illustrated the scale-dependence of the flood response (magnitude, timing, volume) and its relation to the underlying runoff mechanisms. In particular, we identified variations in the event runoff ratio with catchment scale. In the coupled system, the interaction of the water table position with the surface topography was identified as an important determinant of runoff characteristics and influences the shift in hydrologic regimes from surface to subsurface dominance. 

Several important conclusions arise from the analyses of the fully-distributed watershed model response to a series of storm events in a complex, humid basin in Oklahoma. The events under study are representative of storm characteristics in the fall period and pre-storm baseflow conditions arising during the preceding summer. Results from the study indicate the following:

1. A strong relation was found between runoff nonlinearity and the partitioning of total flood response into surface and subsurface components. Transitions in runoff nonlinearity are clearly due to shifts in the dominance of runoff mechanisms. This behavior is robust as it occurs across a wide range of storm properties, initial wetness conditions and catchment scales. 

2. The water table position relative to the surface topography dictates the temporal and spatial distribution of runoff production and flood response in the basin. Depth to groundwater varies with catchment scale leading to different initial conditions in particular subbasins. 

In limited cases, high initial wetness can lead to linear basin response for intense storm events. 

3. Runoff production mechanisms and groundwater dynamics can preferentially occur in specific catchment niches related to a range of topographic indices. Storm properties and the water table position influence the spatial dynamics of runoff production. In particular, saturated regions near channels expand with initial wetness and can contribute to high runoff production. 

4. Flood magnitude, time to peak and the runoff ratio were found to be scale-dependent and approximately follow power law relations over a range of catchment areas for two selected storm types. Scale-dependence of the runoff ratio varies with initial wetness and storm properties and is related to the fraction of forested basin area at each scale. 

5. Scale-dependence was identified in the internal runoff production and nonlinearity for three basin sizes. Increases in catchment scale promote the existence of a diverse set of runoff mechanisms as greater complexity is present in surface-subsurface interactions. Initial conditions modulate runoff production and may lead to runoff linearity for wet cases and large flood events. 

Insights from the modeling experiments reveal that the basin flood response is related to the runoff mechanisms excited as storm properties interact with particular catchment locations and their wetness state. Identification of the intimate link between runoff response characteristics and the underlying mechanisms provides a process-based explanation for nonlinear responses in gauged and ungauged basins. An important result emerging from our modeling exercise is that large changes in basin response occur when the dominant mechanism transitions between surface and subsurface runoff. If runoff partitioning can be properly captured in numerical models, there is the possibility of reproducing observed nonlinear responses across a range of real watersheds. Distributed modeling results also provide a physical explanation for the scale-dependence of runoff generation in complex basins. A testable hypothesis arising from our experiments is that the scale-dependence of the runoff ratio exhibits different regimes which vary according to the underlying mechanisms. Our results indicate that spatial heterogeneities in landscape and initial wetness interact with storm forcing to produce runoff generation patterns that exhibit variations with aggegration scale. Capturing surface-subsurface dynamics in numerical models of gauged and ungauged basins may allow understanding of process controls on runoff scaledependence. 

Understanding the nonlinearity and scale-dependence of hydrologic response and its relation to the underlying runoff mechanisms is important to advance our capability to predict behavior in ungauged basins, an important challenge in catchment hydrology. While our study is focused on a single basin, the watershed exhibits similar hydrologic behavior to other regions in the Great Plains (Garbretch et al., 2004). As a result, basin response characteristics identified in this study may be applicable to similar settings in the broader region. It is important to consider, however, that the behavior in other basins with different climate and surface characteristics may vary from our study results. Despite this, the metrics introduced here to assess the mechanistic causes of catchment response can be useful tools for detailed investigations in other basins or with different numerical models. Examples of these metrics include the linkage between runoff partitioning and nonlinearity over a range of storm conditions and the scaledependence of the runoff ratio for specific storms. Testing the robustness of our results in alternative settings or with other model structures would be a fruitful avenue that may lead to generalizable conclusions on the role played by runoff mechanisms on basin response nonlinearity and scale-dependence. 

The results of this study are also limited to spatiallyuniform storms with no account made for spatial patterns or within-storm rainfall variability. This assumption allowed us to focus on the process controls on runoff production for simple, rectangular rainfall pulses. As a result, the spatial patterns and scale-dependence exhibited in surface and subsurface runoff are due exclusively to the interaction of storm properties with the basin characteristics, including its distributed water table position. Nevertheless, the distributed hydrologic model could be used for exploring basin response under more complex storm forcing. For example, model forcing can consists of either observed radar rainfall or a stochastic rainfall model conditioned on regional observations (e. g., Ivanov et al., 2004b Ivanov et al.,, 2007)). The spatial and temporal variability of the precipitation forcing is expected to resonate with the runoff production mechanisms and the distribution of travel times to influence basin response and its scaling behavior, a topic of current investigation. Ultimately, the distributed hydrologic model can be used as an interpretive tool to assess the surface-subsurface processes that control runoff production resulting from a range of possible forcing conditions. Furthermore, this approach is a useful means to identify the effects of catchment scale on hydrologic response.
section_name: Discussion and conclusions

section_content: Acknowledgements. This work was supported by the U. S. Army Research Office ( ARO DAAD19-00-C-0114 ; ARO W911NF-04-1-0119 ), NOAA ( NA97WH0033 ), and NASA ( NAG5-7475 ). We appreciate the comments of M. Sivapalan and an anonymous reviewer that helped improve the manuscript. 

Edited by: P. Molnar