Effects of micro-topography on surface–subsurface exchange and runoff generation in a virtual riparian wetland — A modeling study

Advances in Water Resources 33 (2010) 1388–1401

Contents lists available at ScienceDirect

Advances in Water Resources j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / a d v wa t r e s

Effects of micro-topography on surface–subsurface exchange and runoff generation in a virtual riparian wetland — A modeling study S. Frei a, G. Lischeid b, J.H. Fleckenstein a,c,⁎ a b c

Department of Hydrology, University of Bayreuth, Germany Leibnitz Centre for Agricultural Landscape Research, (ZALF), Germany Helmholtz Center for Environmental Research (UFZ), Germany

a r t i c l e

i n f o

Article history: Received 4 December 2009 Received in revised form 20 July 2010 Accepted 20 July 2010 Available online 27 July 2010 Keywords: Micro-topography Integrated surface–subsurface simulation Threshold process Spill and fill Transmissivity feedback Runoff generation in wetlands

a b s t r a c t In humid upland catchments wetlands are often a prominent feature in the vicinity of streams and have potential implications for runoff generation and nutrient export. Wetland surfaces are often characterized by distinct micro-topography (hollows and hummocks). The effects of such micro-topography on surface– subsurface exchange and runoff generation for a 10 by 20 m synthetic section of a riparian wetland were investigated in a virtual modeling experiment. A reference model with a planar surface was run for comparison. The geostatistically simulated structure of the micro-topography replicates the topography of a peat-forming riparian wetland in a small mountainous catchment in South-East Germany (Lehstenbach). Flow was modeled with the fully-integrated surface–subsurface code HydroGeoSphere. Simulation results showed that the specific structure of the wetland surface resulted in distinct shifts between surface and subsurface flow dominance. Surface depressions filled and started to drain via connected channel networks in a threshold controlled process, when groundwater levels intersected the land surface. These networks expanded and shrunk in a spill and fill mechanism when the shallow water table fluctuated around the mean surface elevation under variable rainfall inputs. The micro-topography efficiently buffered rainfall inputs and produced a hydrograph that was characterized by subsurface flow during most of the year and only temporarily shifted to surface flow dominance (N 80% of total discharge) during intense rainstorms. In contrast the hydrograph in the planar reference model was much “flashier” and more controlled by surface runoff. A non-linear, hysteretic relationship between groundwater level and discharge observed at the study site was reproduced with the micro-topography model. Hysteresis was also observed in the relationship between surface water storage and discharge, but over a relatively narrow range of surface water storage values. Therefore it was concluded that surface water storage was a better predictor for the occurrence of surface runoff than groundwater levels. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Riparian zones contain dynamic interfaces between ground- and surface water flowpaths [10,27]. It is important to understand the mechanisms that govern hydrologic flowpaths and stream flow generation in riparian zones because nutrient transformations and export are integrally related to the hydrologic dynamics [8,18,34,53]. However, these dynamics can be quite complex [27,53] and are generally poorly understood [28,49]. In humid temperate climates riparian zones are often occupied by wetlands [27,34,36]. Rapid surface and shallow subsurface flows typically dominate runoff generation in riparian wetlands during rainstorms [8,34]. Gibson et al. [17] showed that runoff dynamics

⁎ Corresponding author. E-mail address: jan.fl[email protected] (J.H. Fleckenstein). 0309-1708/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2010.07.006

highly depend on surface storage and interactions between surface water and shallow groundwater. Kværner and Kløve [27] identified distinctly different runoff generation processes with shifts between subsurface and surface flow dominance for low and high flow events. Non-linear relationships between riparian water table depth and stream flow have often been observed [5,15,38,44]. For catchments dominated by matrix flow, these relationships have been attributed to the transmissivity feedback mechanism [3,4,44]. Stream flow originating from matrix flow increases exponentially when the water table rises into soil layers with progressively increasing lateral hydraulic conductivity [3,44]. In systems where shifts between matrix flow and surface flow dominance occur, additional dynamics and nonlinearities have been observed (e.g. [27]). Peat-forming wetlands are often characterized by a hummocky topography with sequences of high points (hummocks) and depressions (hollows) at the sub-meter scale, which will affect runoff generation during transitions between surface and subsurface flow

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

dominance. Effects of micro-topography on infiltration and runoff generation processes were first investigated by Dunne et al. [9]. They showed that hill slope runoff was controlled by an intricate interplay between rainfall intensity, surface flow depth, vegetation cover and the specific micro-topography of the slope. Micro-topography can attenuate and delay surface flows [27], because surface depressions first need to be filled until a specific surface water storage threshold is exceeded and surface flow towards the stream channel can be initiated [2,12]. Tromp-van Meerveld and McDonnell [51] and [52] termed similar threshold dynamics in the generation of subsurface stormflows on bedrock surfaces with micro-topography the “fill and spill mechanism”. Qu and Duffy [40] reported distinct double peaks in hydrographs from single rainfall events, which they ascribed to complex interactions between small scale micro-topography controlled surface runoff in the wetland and subsurface flow. Several modeling studies have addressed the effects of microtopography on runoff dynamics. Dunne et al. [9] used a conceptual approach to simulate overland flow and infiltration processes for uniform sinusoidal micro-topography. They demonstrated that microtopography resulted in significant spatial variability of infiltration and surface flows. Esteves et al. [11] and Fiedler and Ramirez [12] used finite difference solutions to the two-dimensional depth-averaged dynamic wave equations to simulate overland flow and infiltration processes on small plots with micro-topography. Both studies showed that micro-topography strongly affects flow directions, flow velocities and flow depths and resulted in surface flow along well defined micro-channels. Connectivity indicators for surface flow on plots with micro-topography were systematically investigated with a numerical model by Antoine et al. [2]. Each of the aforementioned modeling studies were restricted to surface flows and infiltration and did not account for feedbacks between surface and subsurface flow, an important process in wetlands [8,17]. An exception was the study by Qu and Duffy [40], who used a finite element coupled surface– subsurface flow model to simulate a series of rainfall events for a 0.08 km2 watershed in Pennsylvania. They demonstrated how small scale topography can control local surface saturation and subsequent connectivity of surface flow paths leading to stream flow generation. However, the spatial resolution of the Qu and Duffy [41] model was too coarse to account for micro-topography on the sub-meter scale. Hopp and McDonnell [20] modeled the effects of bedrock microtopography on subsurface storm flow generation from hillslopes. Our work evaluates the complex hydrologic dynamics of a riparian wetland with micro-topography through a virtual modeling experiment. The purpose of the simulations is to examine process dynamics rather than the calibration of a model to a specific field site. We argue that to accurately describe these dynamics a numerical model has to account for overland flow, variably saturated subsurface flow and complex interactions between the surface and subsurface domains. A fully-integrated modeling approach simultaneously solves all of the equations that govern the complex interactions between surface and subsurface. Efficient numerical models that use the fully-integrated approach have become available in recent years (e.g. [25,50]). The fully-integrated, three-dimensional numerical flow model HydroGeoSphere [50] is used here to examine hydrologic dynamics in a virtual riparian wetland with distinct micro-topography (hummocks and hollows). The micro-topographic relief is geostatistically generated for a 10 m × 20 m area at a resolution of approximately 0.1 m based on surveyed micro-topography in a riparian wetland of the small experimental Lehstenbach catchment located in Germany (Fig. 1). The peat-forming wetlands in the catchment, which have a hummocky surface topography, can be classified as fens. The relative elevation differences between hollows and hummocks range between 0.2 and 0.4 m and the hollows are generally interconnected. The peat is between 1 and 2 m thick. Inflows from deeper groundwater are locally diminished by a basal clay layer. At several locations, lateral inflows from adjacent hillslopes are intercepted by small stream channels

1389

bounding the wetlands. Most small streams have their headwaters in the wetlands and practically all the water that reaches the streams either originates in or passes through the wetlands. Previous studies in the Lehstenbach catchment indicated that stream flow generation and solute export are mainly controlled by processes occurring within the riparian wetlands ([1,31–33]). The quick response of catchment discharge to intensive rainfall events occurs via surface and shallow subsurface flowpaths ([34]). Surface runoff within wetland zones is predominantly generated by the saturation excess overland flow mechanism (type Dunne runoff) ([32]). Field observations show a distinct non-linear relationship between discharge in a first-order stream originating in the wetland and the depth to the local water table (Fig. 2). This non-linear relationship suggests that runoff generation in the wetland is controlled by a complex threshold response between the formation of shallow subsurface flows and saturation excess overland flow draining into the channel network. We hypothesize that the interplay between water table depth and surface micro-topography (typical pattern of hollows and hummocks) in the wetland results in distinct shifts between surface and subsurface flow dominance that can explain the observed stream discharge behavior. The main objective of this study is to improve our understanding of the effects of micro-topography on the complex hydrological process dynamics and interactions that govern surface–subsurface exchange and runoff generation in a riparian wetland. Process dynamics are examined in a virtual test case using an integrated surface–subsurface numerical flow model. In particular, the following questions are addressed: 1) What effect does micro-topography in riparian wetlands have on stream discharge generation? 2) Can micro-topography explain the observed non-linear relationship between surface water discharge and the water table depth? 2. Methods The objectives of this study are addressed through a virtual modeling experiment. Geostatistical indicator simulations based on surveyed elevation data are used to represent the micro-topographic structure for a synthetic riparian wetland section typical for the Lehstenbach catchment. Hydrological dynamics in the synthetic wetland are evaluated with the physically based surface–subsurface numerical model HydroGeoSphere [50]. The conceptual idea behind the modeling is similar to the virtual experiments proposed by Weiler and McDonnell [55]. The numerical model is used as a virtual landscape, in which perfect process knowledge is assumed (see e.g. Zehe et al. [57]). 2.1. Representation of micro-topography The spatial structure of the micro-topography for a typical fen in the Lehstenbach catchment was derived from several surveyed transects. A 2D spatial distribution of elevation classes was generated using geostatistical indicator simulations based on Markov Chain models of transition probabilities between categorical data (Transition Probability Geostatistical Software), [7]. The method was originally developed to realistically represent aquifer heterogeneity with discrete transitions between different hydrofacies [7] and it has been widely used for groundwater flow and transport problems (e.g. [14,16,29,56]). This approach provides more options to condition the simulation with field data than traditional geostatistical simulation methods for continuous variables based on variogram analysis. The use of geostatistical simulations provided equally probable generations of different spatial structures of micro-topography and several realizations of each individual structure. The following main steps were used in developing the geostatistical model of micro-topography.

1390

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

Fig. 1. Location of the Lehstenbach experimental catchment (upper panel). The overall hydrology of the catchment is controlled by the structure of the basin (lower plot). Dark grey areas represent forest and light grey areas wetlands, which occupy almost 1/3 of the 4.2 km² catchment area (lower panel). The location of the field site that provided field data is marked with an open circle.

First, the elevation classes representing topographical structures like local depressions (hollows), local maxima (hummocks) or transition zones were defined as indicators. The relative elevations (scaled by the mean value: z−z) were subdivided into five different discrete elevation classes. The upper and lower class bounds and relative frequencies are listed in Table 1. A Markov Chain model of

spatial correlation was developed from the relative frequencies of the different elevation classes, the transition probabilities between them, and the mean length of topographical structures. Estimates of mean structure length and transition frequencies for the micro-topography were obtained from the existing surveyed elevation profiles for the fen area. All profiles were surveyed within a 30 × 30 m plot.

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

1391

Table 2 Transition probability matrix for the five elevation classes in the x- and y-directions. L denotes mean length of the structures (hollows and hummocks), b refers to the background category (transition probabilities for the background category follow from the other entries in the matrix by probability law). x- and y-directions Class

1

2

3

4

5

1 2 3 4 5

2– L 6b 6 6b 6 4b b

b – L 0:25 0:1 0:05

b 0:25 – L 0:1 0:05

b 0:1 0:1 – L 0:25

3 b 0:05 7 7 0:05 7 7 0:25 5 – L

– Note: L = mean length, b = background category/class.

2.2. Surface/subsurface flow simulation

Fig. 2. Observed non-linear relationship between stream discharge at the fen field site (from a first-order stream that has its entire headwater in the fen) and groundwater level (measured in a monitoring well in the peat located about 15 m away from the stream). Data was measured in 10 minute intervals from Dec. 15, 1999 to Aug. 15, 2000. Groundwater levels are scaled to the areas mean surface elevation.

Table 2 shows the embedded transition probability matrices in the horizontal directions used to generate the 2D Markov Chain model within TPROGS [7]. The off-diagonal entries represent transition probabilities from one elevation class to another. All transition probabilities were estimated based on the discretized data from the surveyed elevation data. The entries denoted by b represent transition probabilities for the background category, which are automatically estimated from the transition probabilities for the other classes based on probability laws [7]. Diagonal entries represent the mean length of the structures within the corresponding elevation class. Selftransition rates for each class can be directly derived from their respective mean length [7]. The same transition probabilities were used in the x- and y-directions as no pronounced anisotropy was observed in the field data. Mean lengths were assigned uniformly (L = 0:5 m) for all five elevation classes based on the measured mean length of hummocks and hollows. This model represents a surface with similarly sized hummocks and hollows of different elevations rather than a binary field where each of the depressions and peaks receive a unique elevation value. This was found to be the most representative simulation of the structures observed in the field. Using the sequential indicator simulation routine within TPROGS, 3 realizations of the micro-topography for a 10 × 20 m grid with a spatial resolution of 0.1 m in x- and y-directions were generated. For comparison, a spatial model with a mean length of L = 0:25 m for the structures was generated. This model represents a surface topography with smaller hollows and hummocks.

Table 1 Elevation classes and relative frequencies for defining micro-topography in TPROGS, based on surveyed elevation data from the field site.

Class 1 Class 2 Class 3 Class 4 Class 5

Lower limit [m]

Upper limit [m]

Relative frequency [−]

Topographical structure

−0.1 −0.3 b−0.3 0.1 N 0.3

0.1 − 0.1 – 0.3 –

0.128 0.224 0.224 0.231 0.192

Transition zone Hollow Hollow Hummock Hummock

The numerical code HydroGeoSphere [50] is a fully-integrated finite element surface–subsurface flow model. Variably saturated subsurface flow in porous or fractured media is simulated with the Richards equation in three dimensions (3D). Overland- or stream flow in 2D is represented by the diffusion wave approximation to the depth-averaged dynamic wave equations [50]. Surface–subsurface coupling is implemented using the conductance concept. The conductance concept assumes that the exchange flux depends on the gradient across a coupling interface, the thickness of the interface (coupling length), its relative permeability, and the vertical saturated hydraulic conductivity [50]. The governing equations for surface- and subsurface flow are solved simultaneously via a control volume, finite-element approach [50]. HydroGeoSphere has been applied over a range of spatial scales from the plot and river reach scales [6,21] up to the scale of entire catchments [22,30,48]. A 3D triangular finite-element grid was set up using HydroGeoSphere and the preprocessing software GRIDBUILDER [35]. The grid represents a hypothetical 10 m by 20 m section of a riparian wetland with a surface slope of 0.03 m/m that is 2 m thick at the upstream end (slightly thinner at the downstream end due to the inclined surface) and drains to an adjacent channel segment (Fig. 3). The structure and dimensions of the model mimic conditions at a typical fen in the Lehstenbach catchment with a basal clay layer overlain by a 1.5 to 2 m thick layer of peat. The geostatistically generated micro-topography was superimposed onto the surface plane of the grid. As changes in surface elevation are more gradual than the sharp transitions between different elevation classes in the geostatistical model, elevations in the grid were smoothed in GRIDBUILDER. The smoothing algorithm estimates nodal elevations in the grid based on the average elevation of all connected nodes [35]. The final model, subsequently referred to as “micro-topography model”, consists of 210,000 grid nodes with a node spacing of about 0.1 m in the X, Y and Z directions. The channel is represented by a parabolically shaped cross section (1 m in diameter) draining into the positive Y-direction (see Fig. 3) with a constant slope of 0.03 m/m. For comparison, a model with a planar surface was used as a reference model to simulate hydrological dynamics without micro-topographical structures. This model will subsequently be referred to as the “planar model”. The peat body was represented as homogeneous and isotropic. Heterogeneity was intentionally excluded to clearly separate microtopographical effects from effects induced by heterogeneity. Isotropic, saturated hydraulic conductivity for the peat was estimated at 0.2 m/d, which is in the range of values reported for this site [19] and for typical peat soils in general [26,42]. Soil water retention characteristics for the simulation of variably saturated flow in peat were taken from Price et al. [39]. Boundary conditions for the model were assigned so that water can only leave the model domain at the channel outlet (lower right corner of model in Fig. 3). This was realized via a critical-depth boundary condition at the channel outlet. All other boundaries were set to noflow boundaries with the exception of the upper model surface where variable rainfall rates was applied. As an initialization for the

1392

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

Fig. 3. Geometry of the virtual riparian wetland segment: a) planar reference model showing the main drainage direction and channel location; b) smoothed realization of the wetlands micro-topography; c) cross section (Y = 5 m) of the micro-topography model. The location of a virtual monitoring well, channel outlet hydrograph nodes (to track total channel discharge) and surface hydrograph nodes (to track surface flow into the channel) are only shown for the micro-topography model, but are also the same in the planar model. The planar model is about 1.35 m thick at the upstream end of the domain, the micro-topography model about 1.95 m (micro-topographic relief is superimposed onto the planar surface).

subsurface domain the groundwater level was set to 0.5 m above the horizontal base of the model with an equilibrium pressure distribution above the water table. Daily precipitation was applied to the model surface based on the rainfall record from the 2000 hydrologic year (November 1999 through October 2000). The surface domain was initialized with a zero depth of ponded water representing dry initial conditions. The friction slope for surface flow calculations is described using Manning's equation. Manning's roughness coefficients for the peat surface were uniformly assigned as 0.03 m−1/3 s for x and y; a value reported for high grass [45].

ponds at the surface developed connected channel networks, which eventually spilled into the main channel segment around day 45. The subsequent rapid increase in discharge displayed several kinks, which represented the development and maturation of different surface flow networks. The networks eventually all provided water to the channel when equilibrium was reached around day 50. Groundwater levels at equilibrium (evaluated at the location in the upslope center of the domain — see Fig. 3) were about 0.18 m below the land surface for the micro-topography model and at the land surface for the planar model.

3.2. Runoff dynamics and flow components for variable rainfall 3. Results 3.1. Dynamics of runoff generation for steady rainfall To investigate the general dynamics of discharge generation under increasing wetness, a simulation with a constant rainfall rate of 0.008 m/d was run until the steady state discharge at the channel outlet was attained. The rainfall rate represented conditions for a moderate to intense rainstorm (exceeded on about 40 days per year for a typical hydrologic year) and ensured that surface flow networks could develop before the final steady state was reached. Fig. 4 shows the development of channel discharge and the water table (evaluated at a virtual observation well at the upstream end of the model domain — see Fig. 3) for the planar (upper panel) and the microtopography model (lower panel) respectively. In the initial stage of both simulations, channel discharge gradually increased from subsurface inflows caused by increasing hydraulic gradients towards the channel. The increase was more rapid in the planar model compared to the micro-topography model. The slower and slightly undulating increase in the latter case was caused by the progressive formation of ponds when the water table intersects local surface depressions. At this point, the build-up of subsurface gradients towards the channel was slowed. For the same input of rainfall subsurface heads below the ponds increased less rapidly, as if the same amount of water had infiltrated (due to the porosity). Surface flow in the planar model, indicated by a steep increase in channel discharge, occurred after approximately 16 days. Discharge subsequently increased rapidly until the system attained a state of equilibrium with constant discharge around day 24. In the micro-topography model, isolated

Fig. 5 shows the simulated discharge hydrograph at the main channel outlet (lower right corner of the domain — see Fig. 3) for the micro-topography model. Discharge is separated into a surface and a subsurface flow component. The model-forcing daily precipitation record is depicted on the top axis. The separation of flow components was achieved by placing “hydrograph nodes”, which tracked all flow through a node in the grid, along the edges of the channel segment and at the channel outlet (see Fig. 3). The surface flow hydrograph nodes tracked surface flow into the channel at each time step of the simulation. The subsurface flow components tracked all flow that exited the model domain (sum of surface and subsurface flows). The difference between the two components represents subsurface flows into the channel segment. In the initial phase infiltrating rainfall recharged groundwater, the groundwater levels increased, and the hydraulic gradients to the stream increased, resulting in increased subsurface flows. After initial wetting of the system, surface flow via surface channel networks, was initiated on day 125. Maximum discharge was simulated for day 217 after the most intensive rainfall event in the annual record (48 mm/d). Simulated discharge was generated via subsurface flow during most of the year. Only on 52 of the 365 simulated days was surface flow observed in the model. On these 52 days, surface flow accounted for up to 83% of the total channel discharge. Fig. 6 (panel a) shows a typical situation during periods with low to intermediate rainfall intensities. Water was already ponded in local depressions (hollows) at the surface. However, ponded areas are not all interconnected and surface flow into the channel segment was inhibited by the micro-topography. Only during high rainfall

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

1393

structures and could occur as sheet flow as soon as the water table intersected the land surface. In the planar model surface flow contributed up to 95% of the total discharge during individual events. 3.3. Non-linearities and hysteresis

Fig. 4. Stream flow hydrographs and development of the local groundwater level for a simulation with constant rainfall (0.008 m/d). Results for the planar model are shown at the top and for the micro-topography model on the bottom. Stream flow and groundwater level are evaluated at the channel outlet and in a virtual observation well (as shown in Fig. 3). The thicknesses of the soil profiles at the location of the virtual well in the micro-topography and planar models are different as the micro-topography was superimposed onto the planar surface. Hence the initial depths to the water table are different for both models.

rates (panel b) did ponded areas start to become interconnected and form extended surface flow networks and micro-channels. Under these conditions, a large fraction of the wetland surface drained into the adjacent channel. Drainage into the channel occurred at two distinct locations (Fig. 6). Similar patterns were observed in the fen located at the field site during a rainstorm in the spring of 2009 (Fig. 7). Fig. 8 shows the simulated discharge hydrograph for the planar model. Compared to the micro-topography model, the hydrograph generally showed higher peak discharges. Surface flows were generated much earlier (around day 55) and occurred more frequently compared to the micro-topography model (75 of 365 simulated days). During the relatively dry summer period between days 150 and 217, rainfall intensities during flows the six different events were high enough to generate surface drainage. The micro-topography model, in comparison, did not show any surface flows during this period. In the planar model surface flow was not inhibited by micro-topographic

No unique groundwater level (at the virtual observation well) or rainfall rate could be associated with the development of surface flow networks and the onset of surface flows. In contrast, the amount of ponded surface water, necessary to initiate flow to the channel via the surface flow networks and micro-channels, was narrowly defined. Fig. 9 (upper plot) shows the relationship between surface discharge and surface water storage (total amount of ponded surface water in m³ stored in local depression and flow networks). Fig. 9 (upper plot) summarizes results for the 365 day simulation for the microtopography model. The different loops represent different trajectories for single rainfall events. The trajectory for the most intense rainstorm of the simulated year is discussed in more detail (marked by a line in Fig. 9 upper plot). The precipitation event (48 mm/d) occurred right after an extended drier period (days 150 to 217, see Fig. 5) on day 217 followed by only 7.2 mm/d on day 218. In the beginning of the rain storm, infiltrating rainwater exclusively recharged groundwater (no surface discharge). With rising groundwater levels, local depressions were filled with water and increasingly more water was stored on the wetland surface (increasing surface storage without surface flow to the channel). Later, the filled depressions started to interconnect until a critical surface storage value (~ 5.6 m³) was exceeded. The resulting flow network was subsequently large enough to provide first surface flow to the channel. The surface flow rapidly increased until it reached a stable rate of ~3.8 m³/d. After that surface flow abruptly increased as surface storage exceeded another critical value (~7.2 m³). This was caused by a second surface flow network that developed and drained independently from the first one. That is illustrated in Fig. 10 by different snapshots, taken for five different time steps. The lines on the model surface in Fig. 10 delineate different surface flow networks, determined by an analysis of overland flow stream traces. The first snapshot shows the situation right at the beginning of the rainstorm on day 217. Due to the preceding drier period only a few, isolated depressions were filled with water. After the onset of rainfall, additional depressions were filled with water and the isolated ponded areas began to interconnect (second snapshot). After 12 h and 10 min (third snapshot) networks 1 and 2 reach their maximum extent, although surface drainage to the channel was not yet initiated. Although the third snapshot seems to suggest that there were small ponded flow bridges connecting the two networks, an analysis of stream traces showed that there was no surface water exchange. After 12 h and 20 min (fourth snapshot), the first flow network started to drain into the channel causing the first increase in total channel discharge (Fig. 9 upper panel). 7 h and 10 min later (fifth snapshot), the second flow network was activated and started to spill into the channel resulting in the second rapid increase in discharge (Fig. 9 upper panel). Consequently peak discharge (at the end of day 217) occurred when both networks were connected to the channel. At that time, the third zone was still not connected to either of the two other networks. This area was characterized by depressions, which remained isolated from the flow networks and where the ponded surface water was immobile during events. The process of growing (during rainfall events) and shrinking networks (during flow recessions) was responsible for the different clockwise loops that can be seen in the relationship shown in Fig. 9 (upper panel). In contrast to the wetting process, the drying cycle during flow recessions proceeded much more uniformly (Fig. 9 upper panel) because it was not characterized by the same stepwise threshold behavior as the wetting process. This distinctly different behavior of the system during wetting and drying resulted in the observed hysteresis loops. The same behavior was also evident in the

1394

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

Fig. 5. Simulated, yearly hydrograph for the micro-topography model. Precipitation at the field site for the hydrologic year 2000 (10/31/1999–11/1/2000) is shown on the top. Surface and subsurface fractions of total channel discharge are shown in green and black respectively.

development of subsurface and surface flow throughout the event (Fig. 12). Subsurface flow showed a gradual increase during wetting of the system, whereas surface flows were initiated at distinct thresholds when specific flow network started to spill into the channel. In contrast the recession of surface flow was much more gradual. Fig. 9 (upper panel) also suggests that there is not one critical value of surface storage that, when exceeded, results in surface flow. The figure indicates that there is a range in the surface storage that defines when surface flow is likely to occur. This range was much more closely defined than the range of groundwater levels at a specific location that were associated with the presence of surface flow. For the microtopography model, the minimum (Scritmin) and maximum (Scritmax) surface storage values for generation of surface flow were 4.7 m³ and 7.2 m³, respectively. Within this range, every flow event displayed a unique hysteretic trajectory, as seen in Fig. 9 (upper panel). The surface storage–discharge relationship for the planar model is shown in the lower panel of Fig. 9. Maximum surface storage (Scritmax) before the initiation of surface flow is about 0.41 m³, which is almost 20 times smaller than for the micro-topography model. A very mild form of hysteresis can also be seen in the relationship for the planar model,

but in the other direction (counter clockwise) compared to the microtopography model (clockwise). Hysteretic loops in the planar model only occurred over a range of 0.3 m³ of surface storage, whereas the range in the micro-topography model was approximately ten times higher. The mild hysteretic behavior in the planar model was attributed to the curved free-water table intersecting the straight planar land surface differently during wetting and drying of the system. Fig. 11 depicts the simulated relationship between groundwater level (again evaluated at the observation well shown in Fig. 3) and channel discharge. It shows a non-linear relationship that is very similar to the one observed at the field site (see Fig. 2 — total discharge volumes are different because of different sizes of the contributing areas). When the system was drained by subsurface flows only (indicated by the circles in Fig. 11) large groundwater fluctuations were accompanied by relatively small changes in flow to the channel (see also Fig. 12). When surface flow was initiated by the activation of surface flow networks (indicated by crosses in Fig. 11), drastic increases in channel discharge were realized and the system shifted from pure subsurface flow to surface flow dominance. This behavior can best be visualized for the most intense rainstorm on day 217 and

Fig. 6. Snapshots of the evolving surface flow networks for a) moderate flow conditions (day 180) and b) during peak flow (day 218). Blue zones indicate ponded surface water and yellow arrows stream traces in the surface flow networks. Snapshots show simulated results.

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

1395

Fig. 7. Picture of the field site taken during a storm flow event in spring 2009. Channel location is marked by a line.

the subsequent flow recession on day 218 (Fig. 12 and solid line in Fig. 11). At the beginning of the rain storm, the peat groundwater level was low (approximately 0.3 m below surface), due to the dry antecedent conditions. After the onset of rainfall groundwater levels increased to a depth of around 0.08 m below the surface (Figs. 11, 12) with a minimal increase in channel discharge (Fig. 11). Then the first surface network started to spill into the channel (at time 217.52 in Fig. 11), causing a tripling of channel discharge (from b 1 m3/d to more than 3 m3/d) over a very short time period (see Fig. 12). The first surface flow network matured and surface and subsurface flow to the channel gradually increased (see Fig. 12). At time 217.81 (Fig. 11), the next surface storage threshold was reached and the same process was repeated for the second surface flow network. During the flow recession (starting on day 218), the relationship shown in Fig. 11 took

a different trajectory. The filled surface channel network was rapidly drained into the channel accompanied by a rapid decline in the water table (Fig. 12). The difference between the threshold type behavior of the system with a stepwise activation of the different surface flow networks during wetting and the more gradual flow recession during drying causes the hysteretic relationship between groundwater level and channel discharge shown in Fig. 11. The same hysteretic loops in this relationship were seen for other rainstorms during the wetting and drying of the system (Fig. 11). The general form of this pronounced hysteresis was similar to the one seen in the field data (see Fig. 2) and the one evident in the surface storage–discharge relationship. The planar reference model does not show the same consistent hysteretic loops as the micro-topography model. The transition

Fig. 8. Simulated, yearly hydrograph for the planar reference model. Precipitation at the field site for the hydrologic year 2000 (10/31/1999–11/1/2000) is shown on the top. Surface and subsurface fractions of total channel discharge are shown in green and black respectively.

1396

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

realization. However, in the second model a small surface flow network close to the channel developed more rapidly and started to spill into the channel when the upstream water table was still 0.2 m below the land surface. This stresses the importance of connectivity between the surface depressions that form the surface flow networks, which can be quite different between different realizations of one geostatistical model. The micro-topography model, with a reduced mean length, showed a less pronounced hysteresis and a shape in the groundwater level–discharge relationship that was somewhere between the planar model and the original micro-topography model. 4. Discussion 4.1. Effects of micro-topography on runoff generation

Fig. 9. Relationship between surface storage and channel discharge for the microtopography model (upper panel) and the planar model (lower panel). The black line depicts the peak flow event around day 218. Scrit(max)–Scrit(min) (upper and lower panels) represents the critical range of surface water storage, within which surface flows occur in the yearly simulations.

between subsurface and surface flow dominance in the planar model was characterized by a very sharp and defined break in the relationship, in contrast to the smoother transition in the microtopography model. The relationship between groundwater level and channel discharge for the flow event on days 217 to 218 in the planar model is shown in Fig. 13 (lower right panel) in addition to the same relationship for the micro-topography model (upper left panel). Another geostatistical realization of the micro-topography model (upper right panel) and a model with micro-topographic structures that have half the mean length of the original model (lower left panel) are shown for comparison. The planar model showed no hysteresis in the relationship between groundwater level and channel discharge. The relationship followed the same trajectory during wetting and drying. The second realization of the original micro-topography model showed a similarly pronounced hysteresis as the first

4.1.1. Threshold behavior and flow components Simulation results suggest complex dynamics of runoff generation in riparian wetlands with a distinct micro-topography that experience frequent shifts between subsurface and surface flow dominance. For a steady input of rainfall, stream flow was initially generated from subsurface flow into the stream channel. Generation of surface flow was controlled by a threshold process of surface ponding and the subsequent stepwise formation of different surface flow networks. For variable rainfall input, these networks dynamically expand and shrink and the system frequently shifted between subsurface and surface flow dominance. Threshold controlled dynamics of runoff generation are not uncommon in hydrologic systems [60] and have been described for macro-porous soils [57], cracking clay soils [54,58,59], complex glaciated landscapes [46,47] and subsurface storm flow on bedrock with micro-topography [51,52]. As in the case presented here, non-linear response of threshold systems to external forcing (e.g. rainfall inputs) is typically controlled by distinct shifts in dominant flow regimes (e.g. [57]) and/or a stepwise activation or deactivation of different hydrologic compartments in the system (e.g. [47]). Spence and Woo [47] described such a system as a spill and fill runoff system. The same concept was later used by Tromp-van Meerveld and McDonnell [52] to describe the generation of subsurface storm flow on slopes with a micro-topographic bedrock structure. The depth to groundwater (assessed at the virtual well or averaged over the model domain) proved to be a relatively poor indicator for the development of surface flow to the channel. Depending on the history of the system (e.g. wetting or drying phase) quite different surface flow responses were associated with similar water table elevations (Fig. 11). This behavior is typical for threshold systems. Their response is often only poorly described by the macroscopic average states (e.g. average water table) of the system [58,59]. The amount of water stored on the surface was found to be a better indicator for the occurrence of surface flows (Fig. 9). A narrow range of surface storage values defined when surface flow was initiated (see Fig. 9). One possible reason for the absence of one single threshold value could be the dynamics of re-infiltration in the transition zone between the peat and the channel. In this zone re-infiltration rates were usually high because of steep gradients between the groundwater and the channel, which varied with the dynamics of the system. Therefore, different amounts of water stored on the surface might be needed to generate sufficient surface flow to exceed re-infiltration rates in this zone. The simulated stream flow hydrograph in the micro-topography model was significantly moderated compared to the much flashier response in the planar model. Micro-topography functioned as a filter that buffers rainfall inputs, facilitated dynamic exchange between surface and subsurface and delayed runoff response. Similar moderations of stream flows (e.g. retardation of peak flows due to temporal surface storage of water) have been observed in field settings, e.g. in peat-forming wetlands [13,27] or in other landscapes with complex

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

1397

Fig. 10. Six consecutive snapshots of the evolving surface flow networks during the largest flow event of the year (day 217 to day 218). The red lines separate different flow networks (1–3) that developed independently from each other.

surface topography [46]. As a consequence of this flow moderation subsurface flow was the dominant runoff generating process throughout most of the year (during dry to moderately wet conditions). Total surface flow contributions per year were only 25% in the micro-topography as opposed to 42% in the planar model. Surface flow only became the dominant runoff producing process during intense rainstorms (on 52 out of 365 days in our simulations), when it was responsible for up to 83% of the total discharge generated on a given day. In contrast, in the planar reference model surface flow contributions occurred much more frequently (on 75 days out of 365) and provided up to 89% of the total daily discharge. 4.1.2. Structure and dynamics of the surface flow networks In the simulations with micro-topography several surface flow networks developed that operated independently with no exchange of water between them. Flow depths and velocities in the simulated surface flow networks were found to be highly variable and discontinuous in space and time, which is consistent with findings

Fig. 11. Simulated relationship between groundwater level and channel discharge for the micro-topography model. Circles represent times when no surface flow occurs, crosses represent conditions when surface flow is being generated, the sequence of days 217 to 219, representing an intense rain storm, is depicted by a line.

from other field and simulation studies [11,12]. The specific structure of the micro-topography controlled the timing and magnitude of surface water spills into the channel. These topographic controls and the resulting spill and fill dynamics [51] caused stepwise increases in channel discharge (e.g. Fig. 4). Similar results were presented in other studies on runoff generation in systems with micro-topography [2]. Simulation results further showed that the structure of microtopography resulted in isolated depressions that were never connected to one of the main surface flow networks (Fig. 10). Certain flow networks may only be activated during very large rainstorms and remained stagnant during smaller rainfall events. Similar findings were presented by Antoine et al. [2], who performed numerical surface flow simulations on synthetic micro-topographies with different structures. They concluded that on a larger scale, realistic micro-topographies lead to conditions where connected and unconnected depressions coexist. For variable rainfall inputs the ratio

Fig. 12. Simulated development of subsurface, surface flow and groundwater level during the largest rainstorm (days 217–218). Vertical dashed lines mark the different periods when stream flow is characterized by subsurface flow only, subsurface and surface flow with either network 1 or networks 1 + 2 being active and by flow recession.

1398

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

Fig. 13. Simulated relationships between channel discharge and groundwater level (evaluated in the virtual observation well) for two different geostatistical realizations of the micro-topography model (upper panels), a micro-topography model with reduced mean length of the structures (1/2 of the original model) (lower left panel), and the planar reference model (lower right panel).

between active and evolving networks and stagnant pools may frequently shift, resulting in a highly non-linear flow response of the system. This resonates with the conclusion by Fiedler and Ramirez [12] that representing complex slopes as planar surfaces and ignoring small scale dynamics and interactions between surface and subsurface may lead to a false representation of the hydrological system.

induced dynamics will dominate the upper stream flow range when surface flow has been initiated, because then flow to the stream is no longer constrained by porous media flow. Should the uppermost soil layers be so conductive to allow very large lateral subsurface flows to the stream, it is unlikely that ponding will occur, as infiltrating water is quickly moved away laterally. In this case the effects of micro-

4.2. Micro-topography and non-linear response The observed non-linear relationship between groundwater level and discharge (Fig. 2) could be adequately represented with the micro-topography model (Fig. 11). Similar non-linear relationships have frequently been observed in the field [13,15,23,28,37,38,41,43]. A common explanation for such relationships is the transmissivity feedback effect [4]. It is based on an exponential decay of lateral hydraulic conductivity in the soil with depth. As a result, transmissivity and consequently the rate of lateral water movement increase exponentially as the groundwater level rises into superficial layers [4]. In this study a similar non-linear relationship was generated with a uniform distribution of hydraulic conductivities in the subsurface solely by effects of the micro-topography (Fig. 11). In contrast, the planar model showed a linear increase of discharge with rising water table and a distinct and much more abrupt shift between subsurface and surface flow dominance (Fig. 13). This suggests that in systems that experience frequent shifts between subsurface and surface flow dominance additional mechanisms might be at work to generate the non-linear relationships between the riparian water table depth and discharge commonly observed in the field. In many field settings probably both mechanisms will operate jointly. We suspect that in such cases transmissivity feedback will control non-linear response in the low to moderate stream flow range, whereas micro-topography

Fig. 14. Relationship between discharge and groundwater level for two peak flow events observed for a small catchment located in British Colombia, Canada (modified after Fitzgerald et al. (2003) [13]).

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

topography on the non-linear runoff response will be minimal. Which mechanism dominates will be dependent on the specific field conditions, such as magnitude and spatial variability of subsurface hydraulic conductivity, rainfall duration and intensity, initial water table depth and the structure of the surface micro-topography. Similar clockwise hysteresis as seen in the simulated relationship between water table depth and discharge has been observed in peatforming wetlands ([5,13]) as well as in riparian zones characterized by mineral soils ([23]). A very pronounced hysteresis with very similar characteristics to the ones obtained from our micro-topography simulations was observed by Fitzgerald et al. [13] during a rainstorm in a headwater swamp in British Columbia, Canada (Fig. 14). Similar to our simulated relationship, their study also showed distinctly different hysteresis loops for individual peak flow events. Whereas Branfireun et al. [5] and Kendall et al. [23] attributed the observed hysteretic behavior to different compartments of the system dominating flow to the channel during rising flows and flow recession, Fitzgerald et al. [13] provided no specific explanation for the observed relationship. In their study channel discharge only increased significantly after extended surface ponding was observed, which suggests similar dynamics as observed in our micro-topography simulations. Hysteresis in our model was caused by a distinct sequence of the threshold controlled, stepwise increase in surface flow followed by a more gradual draining of the surface flow networks associated with a rapid decline in groundwater levels. In systems that frequently shift between subsurface and surface flow, which is not uncommon in riparian wetlands in humid climates [13,27], the described micro-topography induced dynamics will superimpose non-linearities caused by other mechanisms such as transmissivity feedback. Whether one process is dominant or several mechanisms operate jointly to produce the observed non-linear response will depend on the specific site characteristics. 4.3. Limitations and constraints The virtual systems simulated here are a simplification of most real field situations. Peat soils are usually quite heterogeneous and hydraulic conductivities are typically non-uniform. Retention characteristics of peat soils can be hysteretic. Hence the non-linear response of real systems may be caused by several reasons. Furthermore higher effective hydraulic conductivities in some peat soils (e.g. caused by preferential flow) may so efficiently drain the wetland that for realistic rainfall rates surface ponding never occurs. Significant inflows from adjacent hillslopes or from deeper groundwater could potentially affect the non-linear response. Most of these aspects were intentionally excluded from this study to specifically highlight the effects of the micro-topography in wetlands that experience subsurface as well as surface flows. This may limit the degree to which the results can be generalized, but it allows insights into the process dynamics caused by distinct surface micro-topography during the transition between subsurface and surface flow dominance. Hummocky topography is not uncommon i peat-forming wetlands and riparian wetlands in humid climates typically experience shifts between subsurface and surface flows. The structure of microtopography, hydraulic conductivities of the peat and rainfall rates were taken from a riparian fen in an experimental watershed in Germany and are representative for other hummocky riparian wetlands in humid climates. Surface ponding and surface flows have also frequently been observed at the field site. Therefore the simulation results can provide valuable insights into the dynamics of runoff generation in such systems that may help to explain observed nonlinear responses in similar field settings (e.g. [13]). The number of test cases that was simulated in this study was limited due to very long simulation times (up to 7 weeks for one year of hydrology). However, a parallel version of the used code has just become available and will facilitate the simulation of more test cases

1399

in the future. The relative importance of different mechanisms (e.g. micro-topography induced dynamics versus transmissivity feedback) under different model configurations (different distributions of hydraulic conductivity, different rainfall patterns etc.) could now be evaluated in additional modeling scenarios. 5. Conclusions Hydrologic systems typically show complex stream flow response to rainfall inputs. Deciphering the processes that cause the observed response is usually difficult due to the strongly non-linear behavior of hydrologic systems [57]. Using physically based numerical models as controlled replicates of natural systems to conduct “virtual experiments” [55] can be a useful tool to elucidate individual processes and their interdependencies (see also Zehe et al. [57]). This approach was used here to investigate the effects of surface micro-topography on runoff generation in a virtual riparian wetland in a humid climate. The virtual system mimics a hummocky riparian wetland that experiences frequent shifts between subsurface and surface flows. Simplifying assumptions (e.g. homogeneous hydraulic conductivity) were made to specifically address the effects of micro-topography on runoff generation. Simulation results revealed complex threshold processes with stepwise expansions and contractions of surface flow networks that govern stream flow generation. Distinctly different behaviors of the system during wetting and drying resulted in a pronounced clockwise hysteresis in the non-linear relationship between stream flow and riparian groundwater level that resembled similar relationships observed in the field. Simulations for different microtopographies and for a planar reference model showed clear differences in the shape of the non-linear relationship and demonstrated how stream flow was moderated by the micro-topography. The planar model did not show significant hysteresis in the stream flow–water table relationship. Results from a model with smaller mean length of the micro-topographic structures (1/2 of the mean length of the original model) suggested that for the decreasing size of the structures the response of the system approaches that of the planar model. A comparison of the model results with results presented by Fitzgerald et al. ([13]) from a field study in a humid riparian wetland in Canada, suggested that the simulated dynamics might provide a consistent explanation for the observed behavior of the system. Due to long simulation times only a limited number of model configurations could be tested in this study. However, a newly released parallel version of the numerical flow model will provide new capabilities to test additional model configurations and to assess the relative importance of different mechanisms (e.g. micro-topography induced dynamics versus transmissivity feedback) that cause non-linear runoff response. We hypothesize that the simulated hydrologic dynamics in wetlands with a defined micro-topography will result in a large range of subsurface residence times and will facilitate enhanced mixing between surface and subsurface water of different ages with potential implications for water quality. Preliminary particle tracking simulations, which will be presented in a follow-up paper, support this hypothesis. To what degree the simulated dynamics could provide a new framework to interpret the common variability in stream water chemistry during events that is described in Kirchner's double paradox [24], remains to be investigated. Future work will also address to what degree simplified conceptual representations of surface structures in numerical models (e.g. by defining a rill storage height for larger model cells) can mimic the effects of the microtopography on surface flow and surface–subsurface exchange. Acknowledgements The authors would like to thank the anonymous reviewers for constructive comments, which helped to improve the final

1400

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401

manuscript. This study was funded by the German Research Foundation (DFG, grant FL 630/6-2). Their financial support is greatly appreciated. The authors also thank Rob MacLaren, Young-Jin Park, Andrea Brookfield and Ed Sudicky at the University of Waterloo, Canada for their invaluable help with the ins and outs of the numerical code HydroGeoSphere.

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.advwatres.2010.07.006.

References [1] Alewell C, Paul S, Lischeid G, Storck FR. Co-regulation of redox processes in freshwater wetlands as a function of organic matter availability? Sci Total Environ 2007;404(2–3):335–42, doi:10.1016/j.scitotenv.2007.11.001 [2] Antoine M, Javaux M, Bielders C. What indicators can capture runoff-relevant connectivity properties of the micro-topography at the plot scale? Adv Water Resour 2009;32(8):1297–310, doi:10.1016/j.advwatres.2009.05.006. [3] Bishop K, Seibert J, Köhler S, Laudon H. Resolving the double paradox of rapidly mobilized old water with highly variable responses in runoff chemistry. Hydrol Process 2004;18:185–9. [4] Bishop, KH. Episodic increase in stream acidity, catchment flow pathways and hydrograph seperation, Ph.D. thesis, Cambridge University, Cambridge, UK 1991. [5] Branfireun BA, Roulet NT. The baseflow and storm flow hydrology of a Precambrian shield headwater peatland. Hydrol Process 1998;12(1):57–72. [6] Brookfield AE, Sudicky EA, Park YJ, Conant Jr B. Thermal transport modelling in a fully integrated surface/subsurface framework. Hydrol Process 2009;23:2150–64, doi:10.1002/hyp.7282. [7] Carle SF, Fogg GE. Transition probability-based indicator geostatistics. Math Geol 1996;28(4):453–76. [8] Devito KJ, Hill AR. Sulphate dynamics in relation to groundwater–surface water interactions in headwater wetlands of the southern Canadian Shield. Hydrol Process 1997;11(5):485–500. [9] Dunne T, Zhang W, Aubry BF. Effects of rainfall, vegetation, and microtopography on infiltration and runoff. Water Resour Res 1991;27(9):2271–85. [10] Duval TP, Hill AR. Influence of stream bank seepage during low-flow conditions on riparian zone hydrology. Water Resour Res 2006;42(10):10425, doi: 10.1029/2006WR004861. [11] Esteves M, Faucher X, Galle S, Vauclin M. Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values. J Hydrol 2000;228(3–4):265–82. [12] Fiedler FR, Ramirez JA. A numerical method for simulating discontinuous shallow flow over an infiltrating surface. Int J Numer Methods Fluids 2000;32(2):219–39. [13] Fitzgerald DF, Price JS, Gibson JJ. Hillslope–swamp interactions and flow pathways in a hypermaritime rainforest, British Columbia. Hydrol Process 2003;17(15): 3005–22, doi:10.1002/hyp. 1279. [14] Fleckenstein JH, Niswonger RG, Fogg GE. River–aquifer interactions, geologic heterogeneity, and low-flow management. Groundwater 2006;44(6):837–52. [15] Fraser CJD, Roulet NT, Moore TR. Hydrology and dissolved organic carbon biogeochemistry in an ombrotrophic bog. Hydrol Process 2001;15(16):3151–66, doi:10.1002/hyp. 322. [16] Frei S, Fleckenstein JH, Kollet SJ, Maxwell RM. Patterns and dynamics of river– aquifer exchange with variably-saturated flow using a fully-coupled model. J Hydrol 2009;375(3–4):383–93. [17] Gibson JJ, Price JS, Aravena R, Fitzgerald DF, Maloney D. Runoff generation in a hypermaritime bog-forest upland. Hydrol Process 2000;14:2711–30, doi:10.1002/ hyp. 1279. [18] Gilliam JW. Riparian wetlands and water quality. J Environ Qual 1994;23(5):896. [19] Hauck, A. Hydrological characterization of the Lehstenbach catchment, unpublished diploma thesis, Department of Ecological Modelling, University of Bayreuth, Bayreuth 1999. [20] Hopp L, McDonnell JJ. Connectivity at the hillslope scale: identifying interactions between storm size, bedrock permeability, slope angle and soil depth. J Hydrol 2009;376(3-4):378–91, doi:10.1016/j.jhydrol.2009.07.047. [21] Jones JP, Sudicky EA, Brookfield AE, Park YJ. An assessment of the tracer-based approach to quantifying groundwater contributions to streamflow. Water Resour Res 2006;42:W02407, doi:10.1029/2005WR004130. [22] Jones JP, Sudicky EA, McLaren RG. Application of a fully-integrated surface– subsurface flow model at the watershed-scale: a case study. Water Resour Res 2008;44:W03407, doi:10.1029/WR005603. [23] Kendall KA, Shanley JB, McDonnell JJ. A hydrometric and geochemical approach to test the transmissivity feedback hypothesis during snowmelt. J Hydrol 1999;219 (3–4):188–205, doi:10.1016/S0022-1694(99)00059-1. [24] Kirchner JW. A double paradox in catchment hydrology and geochemistry. Hydrol Process 2003;17(4):871–4. [25] Kollet JS, Maxwell RM. Integrated surface–groundwater flow modeling: a freesurface overland flow boundary condition in a parallel groundwater flow model. Adv Water Resour 2006;29(7):945–58.

.

[26] Kruse J, Lennartz B, Leinweber P. A modified method for measuring saturated hydraulic conductivity and anisotropy of fen peat samples. Wetlands 2008;28(2): 527–31. [27] Kværner J, Kløve B. Generation and regulation of summer runoff in a boreal flat fen. J Hydrol 2008;360:15–30, doi:10.1016/j.jhydrol.2008.07.009. [28] Laudon H, Seibert J, Kohler S, Bishop K. Hydrological flow paths during snowmelt: congruence between hydrometric measurements and oxygen 18 in meltwater, soil water, and runoff. Water Resour Res 2004;40(3):W03102, doi:10.1029/ 2003WR002455. [29] Lee S-Y, Carle SF, Fogg GE. Geologic heterogeneity and a comparison of two geostatistical models: sequential Gaussian and transition probability-based geostatistical simulation. Adv Water Resour 2007;30(9):1914–32. [30] Li Q, Unger AJA, Sudicky EA, Kassenaar D, Wexler EJ, Shikaze S. Simulating the multi-seasonal response of a large-scale watershed with a 3D physically-based hydrologic model. J Hydrol 2008;357:317–36, doi:10.1016/j.jhydrol.2008.05.024. [31] Lischeid G, Lange H, Moritz K, Büttcher H. Dynamics of runoff and runoff chemistry at the Lehstenbach and Steinkreutz catchment. In: Matzner E, editor. Biogeochemistry of forested catchments a. changing environment: a German case study. Ecological Studies. Heidelberg: Springer; 2004. p. 399–436. [32] Lischeid G, Bittersohl J. Tracing biogeochemical processes in stream water and groundwater using non-linear statistics. J Hydrol 2008;357(1–2):11–28, doi: 10.1016/j.jhydrol.2008.03.013. [33] Lischeid G, Kolb A, Alewell C. Apparent translatory flow in groundwater recharge and runoff generation. J Hydrol 2002;265(1–4):195–211. [34] Lischeid G, Kolb A, Alewell C, Paul S. Impact of redox and transport processes in a riparian wetland on stream water quality in the Fichtelgebirge region, southern Germany. Hydrol Process 2007;21(1):123–32, doi:10.1002/hyp.6227. [35] McLaren RG. GRID BUILDER — a pre-processor for 2D, triangular element, finiteelement programs. University of Waterloo, Waterloo: Groundwater Simulations Group; 2008. [36] Metcalfe RA, Buttle JM. Semi-distributed water balance dynamics in a small boreal forest basin. J Hydrol 1999;226(1–2):66–87. [37] Moldan F, Wright RF. Episodic behaviour of nitrate in runoff during six years of nitrogen addition to the NITREX catchment at Gårdsjön, Sweden. Environ Pollut 1998;102(1S1):439–44, doi:10.1016/S0269-7491(98)80066-3. [38] Molenat J, Gascuel-Odoux C. Modelling flow and nitrate transport in groundwater for the prediction of water travel times and of consequences of land use evolution on water quality. Hydrol Process 2002;16(2):479–92, doi:10.1002/hyp. 328. [39] Price JS, McLaren RG, Rudolph DL. Landscape restoration after oil sands mining: conceptual design and hydrological modelling for fen reconstruction. Int J Min Reclam Environ 2009;99999(1):1–15. [40] Qu Y, Duffy CJ. A semidiscrete finite volume formulation for multiprocess watershed simulation. Water Resour Res 2007;43(8):W08419, doi:10.1029/2006WR005752. [41] Rice KC, Bricker OP. Seasonal cycles of dissolved constituents in streamwater in two forested catchments in the mid-Atlantic region of the eastern USA. J Hydrol 1995;170(1–4):137–58, doi:10.1016/0022-1694(95)92713-N. [42] Schlotzhauer SM, Price JS. Soil water flow dynamics in a managed cutover peat field, Quebec: field and laboratory investigations. Water Resour Res 1999;35(12): 3675–83. [43] Seibert J, Bishop K, Rodhe A, McDonnell JJ. Groundwater dynamics along a hillslope: a test of the steady state hypothesis. Water Resour Res 2003;39(1): 1014, doi:10.1029/2002WR001404,2003. [44] Seibert J, Grabs T, Köhler S, Laudon H, Winterdahl M, Bishop K. Linking soil- and stream-water chemistry based on a iparian flow-concentration integration model. Hydrol Earth Syst Sci 2009;13:2287–97. [45] Shen HW, Julien PY. Erosion and sediment transport. Handbook of Hydrology; 1993. p. 12.1–12.61. [46] Spence C. On the relation between dynamic storage and runoff: a discussion on thresholds, efficiency, and function. Water Resour Res 2007;43(12):W12416. [47] Spence C, Woo M. Hydrology of subarctic Canadian shield: soil-filled valleys. J Hydrol 2003;279(1–4):151–66. [48] Sudicky EA, Jones JP, Park YJ, Brookfield AE, Colautti D. Simulating complex flow and transport dynamics in an integrated surface–subsurface modeling framework. Geosci J 2008;12(2):107–22. [49] Taylor CH. Runoff processes in temperate headwater wetlands. In: Majumdar SK, Miller EW, Brenner FJ, editors. Ecology of wetlands and associated systems; 1997. p. 169–81. [50] Therrien R, McLaren RG, Sudicky EA, Panday SM. HydroGeoSphere a threedimensional numerical model describing fully-integrated subsurface and surface flow and solute transport (manual). University of Waterloo: Groundwater Simulations Group; 2008. [51] Tromp-van Meerveld HJ, McDonnell JJ. Threshold relations in subsurface stormflow: 1. A 147-storm analysis of the Panola hillslope. Water Resour Res 2006;42(2):W02410. [52] Tromp-van Meerveld HJ, McDonnell JJ. Threshold relations in subsurface stormflow: 2. The fill and spill hypothesis. Water Resour Res 2006;42(2):W02411. [53] Vidon PGF, Hill AR. Landscape controls on nitrate removal in stream riparian zones. Water Resour Res 2004;40(3):W03201, doi:10.1029/2003WR002473. [54] Vogel HJ, Hoffmann H, Roth K. Studies of crack dynamics in clay soil: I. Experimental methods, results, and morphological quantification. Geoderma 2005;125(3–4):203–11. [55] Weiler M, McDonnell J. Virtual experiments: a new approach for improving process conceptualization in hillslope hydrology. J Hydrol 2004;285(1–4): 3–18. [56] Weissmann GS.Toward new models of subsurface heterogeneity: an alluvial fan sequence stratigraphic framework with transition probability geostatistics, unpublished Ph.D. thesis, Hydrologic Sciences Graduate Group, Universtiy of California, Davis 1999.

S. Frei et al. / Advances in Water Resources 33 (2010) 1388–1401 [57] Zehe E, Becker R, Bárdossy A, Plate E. Uncertainty of simulated catchment runoff response in the presence of threshold processes: role of initial soil moisture and precipitation. J Hydrol 2005;315(1–4):183–202. [58] Zehe E, Blöschl G. Predictability of hydrologic response at the plot and catchment scales: role of initial conditions. Water Resour Res 2004;40(10):W10202.

1401

[59] Zehe E, Elsenbeer H, Lindenmaier F, Schulz K, Blöschl G. Patterns of predictability in hydrological threshold systems. Water Resour Res 2007;43(7):W07434. [60] Zehe E, Sivapalan M. Threshold behaviour in hydrological systems as (human) geo-ecosystems: manifestations, controls, implications. Hydrol Earth Syst Sci 2009;13:1273–97.