Predicting subsurface stormflow response of a forested hillslope – the role of connected flow paths

Rapid flow processes in connected preferential flow paths are widely accepted to play a key role in the rainfall–runoff response at the hillslope scale, but a quantitative description of these processes is still a major challenge in hydrological research. This paper investigates the approach of incorporating preferential flow paths explicitly in a process-based model for modelling water flow and solute transport at a steep forested hillslope. We conceptualise preferential flow paths as spatially explicit structures with high conductivity and low retention capacity, and evaluate simulations with different combinations of vertical and lateral flow paths in conjunction with variable or constant soil depths against measured discharge and tracer breakthrough. Out of 122 tested realisations, six set-ups fulfilled our selection criteria for the water flow simulation. These set-ups successfully simulated infiltration, vertical and lateral subsurface flow in structures, and allowed predicting the magnitude, dynamics and water balance of the hydrological response of the hillslope during successive periods of steadystate sprinkling on selected plots and intermittent rainfall on the entire hillslope area. The number of equifinal model setups was further reduced by the results of solute transport simulations. Two of the six acceptable model set-ups matched the shape of the observed breakthrough curve well, indicating that macrodispersion induced by preferential flow was captured well by the topology of the preferential flow network. The configurations of successful model set-ups suggest that preferential flow related to connected vertical and lateral flow paths is a first-order control on the hydrology of the study hillslope, whereas spatial variability of soil depth is secondary especially when lateral flow paths are present. Virtual experiments for investigating hillslope controls on subsurface processes should thus consider the effect of distinctive flow paths within the soil mantle. The explicit representation of flow paths in a hydrological process model was found to be a suitable approach for this purpose.


Introduction
Understanding how the internal architecture of hillslopes controls subsurface flow and transport processes and predicting this interplay with models "that work for the right reasons" are still unsolved problems in hillslope hydrology, but also of considerable importance for hydrological predictions at larger scales.
Structures and patterns play a key role in the organisation of hydrological processes across scales (Vogel and Roth, 2003;Schulz et al., 2006;McDonnell et al., 2007).In the context of soil hydrology, it is well known that structural features like pipes and macropores generated by plant roots and animals, or soil cracks from desiccation, offer much less resistance to gravity-driven flow than the surrounding soil matrix, and hence allow rapid flow and transport rates, which has led to the term "preferential flow" (Beven and Germann, 1982;Flury et al., 1994).Together with the bedrock and the soil matrix, these preferential flow pathways determine the subsurface flow characteristics of a hillslope (Peters et al., 1995;Buttle and McDonald, 2002;Uchida et al., 2005;Kienzler and Naef, 2008).Connected networks of preferential flow paths facilitate rapid vertical and lateral transport of water and solutes in the subsurface over considerable distances (Sidle et al., 2001;Anderson et al., 2009a;Wienhöfer et al., 2009a;Baram et al., 2012).This occurrence of preferential flow is bound to the existence of distinctive void structures (Sanders et al., 2012), although the actual preferential flow path may involve the soil matrix around active macropores (Lamy et al., 2009).These rapid flow processes pertain to the runoff mechanism termed subsurface stormflow (Weiler et al., 2006), which dominate the processes involved in runoff generation in response to heavy rainfall at the scale of hillslopes and small catchments, especially at steep hillslopes in temperate humid climates (Bonell, 1993;Uchida et al., 1999;Weiler and McDonnell, 2007;Zehe and Sivapalan, 2009;Jones, 2010).Vertical preferential flow is furthermore an important determinant for leaching and fate of agrochemicals through the vadose zone and related TO soil and groundwater pollution (Flury et al., 1995;Zehe and Flühler, 2001;Clothier et al., 2008).Fast lateral flow has also been related to slope stability (Uchida et al., 2001;Lindenmaier et al., 2005;Hencher, 2010;Wienhöfer et al., 2011;Krzeminska et al., 2012).Consequently, it is highly relevant to incorporate preferential flow processes in models for predicting flow and transport through the vadose zone.Our understanding and conceptualisation of preferential subsurface flow processes, however, is still incomplete.Representing preferential flow processes is challenging from the profile to the hillslope and catchment scales, because both an adequate physical theory linking all types of flow and the observational techniques that could provide the required scale dependent parameterizations are still lacking (Beven and Germann, 2013).Virtual experiments that combine computer modelling with available field evidence are a promising way to advance research on this topic (Weiler and McDonnell, 2004), especially to gain a predictive understanding of the interactions of structures such as macropores, soil layers, and bedrock topography, and the role of this internal architecture for determining the subsurface flow response at the hillslope scale.Nevertheless, current physically based and conceptually based models often avoid the challenges of conceptualizing and parameterizing the effects of lateral preferential flow on gauged and ungauged hillslopes (Weiler and McDonnell, 2007).
Deterministic, spatially explicit models are widely used for simulation of water and solute transport in soils from the core to the field scale.The representation of macropores, that is, pores with equivalent diameters of more than 1 mm or even much larger structures (Beven and Germann, 1982;Luxmoore et al., 1990), and their hydraulic effects has been the subject of numerous studies.As a result, a variety of modelling concepts have been proposed, which are covered in detail by excellent review articles (Šimůnek et al., 2003;Jarvis, 2007;Gerke, 2006;Köhne et al., 2009).The model concepts to account for preferential flow range from alterations of the classical Darcy-Richards model by modification of the unsaturated hydraulic conductivity functions (Durner, 1994;Zehe et al., 2001;Kelln et al., 2009), to dualdomain models that conceptually split the soil into a matrix and a preferential flow domain (Gerke and van Genuchten, 1993;Tsutsumi et al., 2005;Stadler et al., 2012).
Generally, spatially explicit approaches are based on the concept of a Representative Elementary Volume (REV) as the basic element of the model.The parameterization of the REV reflects a representative spatial average of the actual structure of the pore-space, which has a length scale typically much smaller than the spatial discretization of the model.In the same way, the spatial configuration of macropores is represented implicitly in most models, even when flow processes in micropores and macropores are conceptually separated in different domains.Some studies, however, have incorporated preferential flow structures explicitly as discrete fine-scale elements within a spatially explicit model in order to geometrically separate preferential flow paths from the micro-structure of the soil.This strategy has been adopted in numerical experiments to investigate the role of soil pipes for subsurface stormflow (Nieber and Warner, 1991), the role of earthworm burrows for dissipation of free energy (Zehe et al., 2010a), and the effect of disconnected macropores on preferential flow (Nieber and Sidle, 2010).Moreover, the approach has been successfully tested against experimental data; for example for modelling lab-scale experiments with soil cores containing artificial vertical macropores (Allaire et al., 2002b, a;Castiglione et al., 2003;Lamy et al., 2009) and sloping soil blocks containing artificial lateral pipes (Kosugi et al., 2004;Tsutsumi et al., 2005).
Although an explicit consideration of macro-structures is conceptually appealing and instrumental in investigating hydrological processes across scales (Vogel and Roth, 2003), the approach has been rarely tested with field experiments, because detailed information on subsurface flow paths is typically not available.It has to be acquired in the field by marking flow paths with dye or other substances and carefully excavating the soil (Noguchi et al., 1999;Anderson et al., 2009b;Abou Najm et al., 2010).These methods are destructive and require tremendous efforts when applied at scales larger than the profile scale.The application of non-invasive geophysical imaging techniques is promising (Samouelian et al., 2003;Tabbagh et al., 2007), but these are currently not able to resolve preferential flow paths in the field (Moysey and Liu, 2012;Greve et al., 2010;Bievre et al., 2012).It has been shown that random placement of structures (Weiler and McDonnell, 2007) and genetic modelling of structure formation (Vogel et al., 2006) are promising ways of representing structures in process-based models when direct information is not available.This route was recently followed by Klaus andZehe (2010, 2011) for modelling a field-scale transport experiment at a tile-drained site.They tested different realisations of stochastically generated structures to represent vertical earthworm burrows in a 2-D model.Several of these realisations performed equally well in simulating the flow response (Klaus and Zehe, 2010), and tracer transport was acceptably reproduced by a subset of these behavioural model architectures (Klaus and Zehe, 2011).
In this paper we adopt and refine the modelling concept of Klaus andZehe (2010, 2011), and examine its applicability for modelling a hillslope-scale sprinkling and tracer experiment at a natural forested site (Wienhöfer et al., 2009a), where plot-scale observations of prominent vertical and lateral macropores have been linked to very fast hillslope-scale transport of water and solutes.We conceptualise these flow structures as elements with high hydraulic conductivity and low retention capability at a fine spatial resolution, and implement different combinations within the Richards-based CATFLOW model to simulate the hydraulic response of the hillslope to steady-state sprinkling and transient natural rainfall, as well as tracer transport.
The general objective of this study is -to further explore the approach of explicit representation of structures for physically based modelling at the hillslope scale, while the specific objectives are -to model the hillslope-scale tracer and sprinkling experiments, and -to investigate the hydrological functioning of preferential flow paths at the study site.
For these purposes, it is essential to compare the results of the simulations with measured data, and to critically examine the underlying perceptual, conceptual and numerical models and the relations between them.

Study site and relevant field observations
The focus of this paper is on hillslope-scale modelling of rapid flow and transport processes observed at a natural forest site (Fig. 1).The hillslope we seek to model is located in the study area Heumöser in the Vorarlberg Alps (Austria), for which a short overview is given below; further information is provided by Lindenmaier et al. (2005), Lindenmaier (2008) and Wienhöfer et al. (2009aWienhöfer et al. ( , 2011)).The Heumöser belongs to the headwater catchments of the Ebniterach, which is the main tributary of the river Dornbirnerach that drains into the Lake Constance parallel to the Rhine, and is situated 10 km south-east of the city of Dornbirn and 0.5 km south of the village of Ebnit (47  (Schneider, 1999).The friable Amden marl shows no significant fracturing, and thus presents an aquitard with a low hydraulic permeability (Lindenmaier, 2008).
The hillslope represents a well-defined subcatchment of 1232 m 2 on the steep side slopes in the south-western part of Heumöser.The vegetation consists of loose stands of common spruce (Picea abies) and sycamore (Acer pseudoplatanus), and herbaceous understorey.Slope angles vary between 18 and 54 • (median: 30 • ).The subcatchment is the source area of a perennial spring, and was considered a key area for understanding subsurface flow processes that possibly influence slope movement in the central parts (Lindenmaier, 2008;Wienhöfer et al., 2011).This motivated a couple of field investigations to collect information on subsurface characteristics using soil sampling, lab and in situ measurements, and combined sprinkling and tracer experiments.The findings relevant to this modelling study are summarised in the following; some of these have been published in further detail before (Wienhöfer et al., 2009a, b).

Soils and their hydraulic properties
Soils are siltic and vertic Cambisols in the mid-slope, and stagnic and gleyic Cambisols and Gleysols at the hillslope toe.Porosities in the topsoil (0-10 cm) are between 0.48 and 0.73, with a median of 0.58.Bulk densities are low and range from 0.5 to 1.1 g cm −3 , with a median of 0.63 g cm −3 .Soil texture is sandy loam.Below a depth of 10 cm soil texture is significantly finer, classified as silt loam and silty clay loam.Soil depths were measured with a manual auger at 63 locations, and were found to vary between 0.12 m to > 1.10 m (median value 0.70 m; at 8 locations bedrock was not reached at 1.10 m depth).There was no clear trend in variation of soil depth with measurement position along the slope line.
The bulk hydraulic conductivity of the soil was measured in situ under field-saturated conditions with a compact constant head permeameter, and was found to decrease from values around 2 × 10 −5 m s −1 at 12.5-25 cm depth to 10 −6 to 10 −7 m s −1 at 30-100 cm depth.The measured values include the hydraulic effect of macropores, which becomes apparent by the fact that the device's maximum measurable outflow rate of ca. 1 × 10 −4 m s −1 (Sobieraj et al., 2004) was exceeded at one-fifth of the measurement locations (n = 41).The bulk values measured in the field were corroborated by laboratory constant head permeability tests on two large (0.3 × 0.3 × 0.8 m 3 ) undisturbed soil columns from the hillslope.Additionally, multistep outflow experiments were performed on an undisturbed soil column (0.3 m diameter, 0.72 m height) under unsaturated conditions.These allowed determination of soil hydraulic parameters of the soil matrix; for example, average saturated hydraulic conductivity of the soil matrix between 0.36 to 0.72 m depth was determined to be 1.77 × 10 −7 m s −1 (Germer, personal communication, 2011).
The relevance of vertical and lateral preferential flow via macropores at this hillslope was further supported by plotscale dye staining experiments.Dye infiltration was spatially uniform in the upper (0-15 cm) organic-rich soil layer; flows converged vertically into desiccation cracks with apertures up to 1.5 cm and root pipes with diameters of up to 4.8 cm in the lower horizons.Besides vertical percolation, also lateral flow of the infiltrated dye in cracks, horizontal root pipes, and along the bedrock surface was observed.On average, 18.75 % of the plot area were stained at depths below 5 cm (experiment BB1; Wienhöfer et al., 2009a).Prominent pipes were also observed (visually) during excavation of the soil columns in the study area and at the location "cut-bank" described below.

Tracer and rainfall simulation experiments
Tracer experiments at the site showed that these distinct structures form a preferential flow network which generated fast subsurface transport at the hillslope scale.These experiments, which are described in detail by Wienhöfer et al. (2009a), involved rainfall simulation with sprinklers at four plots located along the slope line, tracer application at these plots, and measurements of tracer concentrations and discharge at the hillslope toe (Fig. 1).We focus in this paper on the hillslope tracer experiment conducted in 2007 and measurements taken at the location "cut-bank" (approximately 1.50 m high, 1.20 m wide), where a hiking trail cuts the hillslope and water was observed seeping out from soil pipe outlets located 0.3 to 0.5 m above the bedrock.The seepage from the cut-bank was funnelled into a temporarily installed V-notch weir equipped with a pressure gauge.
Two periods of sprinkling with 12 mm h −1 on a total area of 106 m 2 (the area of the plots in upslope direction were 15.1, 28.4,33.0, and 29.2 m 2 , respectively; Fig. 1) produced a nearly steady-state discharge of 0.08-0.10L s −1 ; no surface runoff was observed.The duration of the first rainfall simulation period was 32.58 h (total input 41.33 m 3 ), and the duration of the second rainfall simulation was 24.75 h (total input 31.39 m 3 ); the time between the two simulations was 12 h.Natural rainfall with a total of 94 mm occurred during the 48 h after the rainfall simulations (total input 115.81 m 3 ) and produced a much higher discharge (Fig. 3).From the available tracer data documented in Wienhöfer et al. (2009a), we have chosen the experiment "Uranine 1" for this modelling study because it was the first experiment under steady discharge conditions and resulted in a smooth breakthrough curve.In this experiment, the fluorescent dye tracer uranine was applied 7.33 h after the beginning of the first sprinkling period at the third sprinkler plot 28.7 m uphill from the cut-bank (Fig. 1).Tracer breakthrough was fast, as in all of the tracer tests at this site; in this case, breakthrough and peak velocities were 1.04 × 10 −2 m s −1 and 3.95 × 10 −3 m s −1 , respectively.
The breakthrough curve was analysed by the method of moments and fitting to a one-dimensional convectiondispersion model.The parameters obtained from the moments of the travel time probability density function resulted in low Peclét numbers (3.3 for experiment "Uranine 1").This illustrates that flow and transport after a distance of almost 30 m was still in the "near field" and far from being well mixed.Yet, this approach did not allow further conclusions on the underlying structures and processes.This analysis indicated, however, that all measurements at the location cutbank under steady-state conditions sampled the same flow field, and that the uranine was not retarded, e.g.due to reversible adsorption, compared to a conservative salt tracer (sodium chloride), which was applied at different plots in shorter distances during later experiments (Wienhöfer et al., 2009a).Another important aspect was the low recovery of the tracer; only 2.93 % of the total applied tracer mass was recovered in the breakthrough curve of the experiment "Uranine 1" considered in this paper.The recovery of uranine in a soil column experiment with an undisturbed soil block (surface area 0.25 m × 0.25 m, depth 0.35 m) was only 22 % after eight days leaching (Wienhöfer et al., 2009a).This apparent loss of tracer was attributed to irreversible sorption in the topsoil.Consequently, at maximum 22 % of the input mass should be expected to be mobile and able to be recovered at the hillslope scale, provided that all flow paths were completely sampled.Even when correcting for this immobile fraction of tracer, the experiment "Uranine 1" would have yielded a recovery rate of only 13.32 %, which could have been due to additional irreversible sorption along the much longer transport distance, or incomplete sampling of flow paths.

Conceptual model of discrete preferential flow paths at the study site
The field observations summarized above point out that prominent macropores in form of pipes and cracks constitute a connected network of vertical and lateral preferential flow paths within the fine-textured soils at the hillslope.A conductive top soil layer of low density and the soil-bedrock interface were additionally observed to influence infiltration and subsurface flow.The perception of the presence and characteristics of these different features stemmed from direct observations at separate spots.We did not have direct information on their spatial configuration over the extent of the hillslope.Nevertheless, with the observed fast breakthrough of the tracers in the experiments it is straightforward to hypothesize that the different structures form a connected preferential flow network which spans the entire hillslope.This implies that vertical and lateral pathways are present at many, if not all segments of the hillslope, and that these structures directly or indirectly connect with each other.This is similar to the conceptual models of subsurface flow paths in a forested slope segment presented by Noguchi et al. (1999) and Sidle et al. (2001).

Modelling approach
The objective of the paper is to test whether an explicit consideration of discrete preferential flow structures allows successful reproduction and prediction of water flow and solute transport at the hillslope described above.To this end we employ the numerical modelling software CATFLOW, and test different spatial model set-ups that are consistent with the available field observations.

Numerical model
CATFLOW is a physically based model for simulation of water and solute transport at the hillslope and catchment scale (Maurer, 1997;Zehe et al., 2001), which has been applied successfully in a number of studies (Zehe and Blöschl, 2004;Zehe et al., 2005Zehe et al., , 2010b, a;, a;Gräff et al., 2009;Klaus andZehe, 2010, 2011).The basic modelling domain in CAT-FLOW is a hillslope, which is represented in the model as 2-D cross section along the line of steepest descent.The third dimension perpendicular to the slope line is only represented by the width of the slope for each node; otherwise, uniformity is assumed.The 2-D profile is discretized by curvilinear orthogonal coordinates, and soil water dynamic is described by the Richards equation in its potential form, which is solved numerically by an implicit mass conservative Picard iteration scheme (Celia et al., 1990).The simulation time step is dynamically adjusted to achieve optimal convergence of the iteration scheme.Soil hydraulic functions can be described using several parameterizations; in the simulations presented here the parameterization after van Genuchten (1980) and Mualem (1976) has been employed.
Rainfall is partitioned into throughfall and interception storage, from which water may evaporate.Evaporation and transpiration are simulated based on the Penman-Monteith equation, taking into account the annual cycles of plant morphological and plant physiological parameters, albedo as a function of soil moisture and the impact of local topography on wind speed and radiation.In the case of infiltration excess or saturation excess, surface runoff is routed along the slope line using the diffusion wave approximation of the 1-D Saint-Venant equation, which is solved numerically with an explicit upstream finite difference scheme.Solute transport is simulated in CATFLOW with a particle tracking scheme based on a random walk approach.The deterministic part of a particle step is determined by the current seepage velocities in each principal direction of the curvilinear grid using a backward two level Runge-Kutta scheme (Roth and Hammel, 1996).The random part of the particle step involves the time step, a dispersion coefficient and a uniformly distributed random number in the interval (−1, 1).In the original version of the CATFLOW code, the seepage velocities acting on a particle at its current position are interpolated from the surrounding simulation nodes, ensuring a continuous velocity field.To model flow and transport in distinct structures, which are represented by individual simulation nodes, we sought to preserve the sharp contrasts in seepage velocities between adjacent nodes, for example macropores and soil matrix.We have therefore slightly modified the CAT-FLOW code; in the version used in this study, the seepage velocities for the particle step are not interpolated between nodes, but the seepage velocity of the actual simulation cell is used.To account for the fact that the distinctive structures in the 2-D profile did not extend over the entire width of the hillslope in the third dimension, the factor f a is introduced, representing the fraction of the macroporous cross section.The factor f a is used as a scaling factor in the determination of the seepage velocities for advective transport as the ratio of the filter velocity and the active cross section.Solute transport via surface runoff is not implemented in either version of CATFLOW.

Model set-up and structure generation
For the implementation of the perceptual model in CAT-FLOW, we adopted the concept of Klaus and Zehe (2010) of representing preferential flow paths explicitly as an artificial porous medium with low hydraulic resistivity (i.e.high hydraulic conductivity and low retention properties).This approach has also been followed by other studies (Nieber and Warner, 1991;Castiglione et al., 2003;Lamy et al., 2009;Nieber and Sidle, 2010).In contrast to Klaus and Zehe (2010), who used the approach for modelling a tiledrained field site with explicit representation of earthworm burrows, our study site is a steep forested hillslope, which is characterised by a shallow soil cover and macropores in form of pipes and desiccation cracks (Wienhöfer et al., 2009a).These structures were conceptualised being less tortuous and distributed more regularly over the hillslope compared to earthworm burrows, and we used a variable spatial resolution of the model with 0.05 m for the preferential flow paths and their surroundings, which appears more realistic compared to the constant resolution of 0.3 m used by Klaus and Zehe (2010).Further differences are that we assigned a uniform set of soil hydraulic parameters to each of the different types of material, without any random components, and that we did not apply any scaling of the width of the model domain to match the peak heights of the hydrographs.
Because the configuration of preferential flow paths was not known, we tested different conceptual representations of flow path structures.In doing so, we sought to vary the configuration of "structures", while keeping fixed what we assumed to be known, e.g.surface topography, soil parameters, or rainfall input.Five different types of structural features that had been observed to facilitate preferential flow at the study site were considered in the simulations: -a loose and litter-rich top soil layer, hereafter referred to as "litter layer", -vertical preferential pathways, hereafter referred to as "vertical structures", -a lateral preferential pathway, hereafter referred to as "lateral pathway", -a preferential pathway at the soil-bedrock interface, and -bedrock topography with two different representations (constant and variable) of measured soil depths.
The starting point for setting up the models was a simulation grid in fine spatial resolution, which is necessary to specify the preferential flow structures explicitly.We chose a grid size of 0.05 m × 0.05 m for the initial discretization of a cross section with a horizontal length of 65.0 m and a thickness of 1.8 m.The surface topography, and thus the geometry of the upper boundary, was taken from a laser-scan digital elevation model with 1 m resolution.The geometry of the lower boundary was defined by shifting the upper boundary by the thickness of 1.8 m perpendicular to the start of the slope line (Fig. 2).This model geometry was then combined with implementations of structures at the respective grid nodes.The litter layer was assigned at the topmost row of the simulation grid.As the outer nodes of the simulation grid are considered with only half of the discretization distance, the thickness of this layer was 0.025 m.The cracks and pipes within the soil matrix were conceptualised as vertical and lateral pathways in the two-dimensional cross section.As the exact configuration of these structures remained unknown, they were generated using random components.We used a Poisson process to allocate the starting points of the vertical structures sequentially along the soil surface, and specified three different minimum distances (1 m, 2 m, and 4 m) between two neighbouring starting points.The variation of this minimum distance effectively determined the density and thus the total number of vertical structures.While extending the structures stepwise into depth, a lateral step was allowed with a probability of 10 % in order to make the structures slightly tortuous.The final depth of the structures was drawn from a normal distribution with a mean depth of 0.9 m and a standard deviation of 0.05 m in order to generate vertical structures that majorly extend down to the mean bedrock depth of 0.85 m, while allowing for some small variation that also produced some structures ending in the soil matrix, especially with the variable bedrock topography.In either case, the vertical structures were cut off when crossing other structures (lateral pathway, soil-bedrock interface, bedrock).A lateral pathway within the soil matrix was generated in a similar manner, starting at the right boundary at a depth of 0.45 m, which corresponds to one-quarter of the total thickness of the modelled profile.This structure was extended stepwise towards the left boundary, allowing for upward and downward steps with a probability of 3 % each while keeping a minimum separation of 2 m between two bendings.To ensure comparability between different model set-ups, a constant random seed was used for generating the stochastic components.To determine the grid nodes with bedrock material, the measured soil depths were interpolated using ordinary kriging.In the simulations we used either a variable bedrock topography obtained by mapping the line of steepest descent of the interpolated bedrock topography onto the 2-D cross section, or a constant soil depth of 0.85 m, which was the mean soil depth of the variable topography (the standard deviation of the variable soil depth was 0.11 m).A soilbedrock interface was implemented as a continuous layer framing the resulting bedrock topography (Fig. 2).
These basic variants for the five types of structures (Table 1) resulted in 64 possible combinations, which formed the "basic model variants" together with the respective parameter values given in Table 2.The 64 basic set-ups (32 with variable and constant soil depth, respectively) formed the core for the subsequent analysis, and were complemented with several modifications.These additional model variants were made in a directed way in order to investigate the effect of an "additional model variant" (Table 1) in comparison with the "basic model variants", but not all possible combinations of all tested modifications were tested.We tested the effect of widening the laterally oriented structures (lateral pathway, soil-bedrock interface, litter layer) from a thickness of one node in the basic set-up to two or four nodes, and we tested different combinations of these wider structures with the basic variants for the other structures and a variable soil depth, resulting in 36 additional set-ups.Six additional setups were obtained by combining very densely arranged vertical flow paths having an average spacing of 0.5 m in all pos-sible ways with the basic variants and a variable soil depth.
In four set-ups, vertical structures were limited to the upper half of the hillslope and combined with a litter layer and/or a lateral pathway, together with a variable soil depth and a soil-bedrock interface.Higher hydraulic conductivities for the structures (Table 1) were tested with seven selected setups.Finally, five "homogeneous" set-ups completed the set of tested configurations.In four set-ups, a homogeneous soil mantle with the parameters of the litter layer and the preferential structures (Table 2) was combined with the constant and variable soil depth, respectively.The last set-up was a "null set-up" containing none of the basic structural elements (i.e. the soil matrix parameters (Table 2) were used for the entire model domain).In total, 122 different model set-ups (64 basic model variants and 58 additional model variants) were simulated, plus a number of preliminary test runs that were performed beforehand and helped in defining the final modelling procedure.
After combining the various structure realisations with the base geometry, the initial discretization was thinned out in model regions without preferential pathways in order to reduce the total number of nodes and thus the computational cost of the simulations.The fine grid size of 0.05 m was retained in the horizontal dimension for the vertical structures including the adjacent matrix nodes; in the vertical dimension it was kept at 0.05 m for the topmost three rows, for the lateral pathway and the soil-bedrock interface and the rows directly adjacent to these structures as well as for the endings of the vertical structures.For all other nodes, the spacing was widened up to a maximum of 0.5 m in the horizontal and 0.15 m in the vertical dimension.All pre-and postprocessing steps were carried out with help of the R environment (R Development Core Team, 2011).

Parameterization of soil and structures
The hydraulic properties of the different materials were modelled with a van Genuchten-Mualem parameterization.For parameterization of the soil matrix we used a parameter set that had been determined by multistep-outflow experiments on large (0.108 m 2 ) undisturbed soil columns from the centre of the hillslope (K.Germer, University of Stuttgart, unpublished data).These parameters had been determined under unsaturated conditions to exclude hydraulic effects of macropores to the greatest possible extent.The parameters for the macroporous structures were chosen to represent a material with low flow resistivity and water retention following Castiglione (2003) and Klaus and Zehe (2010).The litter layer was likewise parameterized as a highly conductive medium with high porosity, whereas a low hydraulic conductivity and a low porosity were assigned to the bedrock material (Table 2).
The transport parameters were chosen to model an ideal and nonreactive tracer.The isotropic effective dispersion coefficients of the different materials were chosen to include the effect of molecular diffusion and hydromechanical microdispersion due to sub-scale structures.Hydromechanical macrodispersion was not considered in the dispersion parameters, as macrostructures were modelled explicitly in this study, and thus rather low dispersion coefficients were selected.The highest dispersion coefficient was chosen for the soil matrix, whereas the value for the bedrock was only twice as high as the diffusion coefficient in water, which for uranine is of the order of 5 ×10 −10 m 2 s −1 (Casalini et al., 2011).All parameter values are given in Table 2.
A scaling factor f a = 0.1875 was specified as initial value for the fraction of the macroporous cross section in the solute transport simulations.To check this value, f a was varied from 0.025 to 0.25 in steps of 0.025.This yielded a set of 10 different factorsf a , which were tested with the set-ups found acceptable for the water flow simulations.

Sequence of simulations and boundary conditions
The various model set-ups were subjected to a succession of simulations, namely two one-week spin-up runs, the simulation of the sprinkling phase of the experiment, during which there was only input at the four experimental plots, and simulation of the natural rainfall phase, which occurred after the sprinkling experiment and during which the rainfall forcing comprised the entire hillslope.
A constant width of 1.75 m, corresponding to the width of the tracer application plot, was assigned for the spin-up runs and the simulation run of the sprinkling phase ("simulation area I", Fig. 1).The natural rainfall phase was simulated with variable widths along the slope line, representing the shape of the subcatchment with a total surface area of 1231.58 m 2 ("simulation area II", Fig. 1).
The final states of the preceding run served as initial condition for the following run.The first spin-up run was started from field-saturated conditions, and was then rerun starting from the simulated final conditions to produce the initial conditions for the simulation run of the sprinkling phase.To determine the initial conditions for the total area run from the final state of this plot-scale run, we calculated a weighted average of the water contents and solute concentrations of the areas affected and not affected by sprinkling, respectively.This was done for each soil type individually.
The boundary conditions at the surface were determined using meteorological data from the climate station at Heumöser, which is located approximately 250 m to the north-west, the known sprinkling rate during the experiment and rainfall data from a tipping bucket rain gauge located next to the hillslope.Plant transpiration was simulated assuming a uniform root distribution over the soil profile and a parameterization for coniferous forest provided by Lindenmaier (2008).A free outflow boundary condition and a gravitational flow boundary condition were prescribed at the toe of the slope (right boundary) and the lower boundary, respectively.

Model evaluation
To evaluate the simulation results, the observations were compared with total simulated runoff, calculated as the sum of surface runoff and water fluxes across the right boundary from the consecutive simulation runs of the sprinkling phase and the natural rainfall phase.Runs were deemed acceptable when they showed a Nash-Sutcliffe efficiency (NSE) greater than 0.75, and matched the observed water balance by 10 %.As no significant amount of surface runoff had been observed during the sprinkling experiments, model set-ups with a surface runoff ratio greater than 10 % of the total runoff during the sprinkling phase were discarded.Solute breakthrough curves were taken from the simulated solute transport over the right boundary of the model domain.For comparison of simulated and observed solute breakthrough, we calculated The identifier and NSE are given in the upper left corner of each panel.The total outflow is the sum of surface and subsurface runoff.In case of the acceptable model set-ups (a-e) total outflow mainly was subsurface flow, while surface runoff was all dominant for the set-up without any structures (f).The rainfall input (g, h) was either provided by the sprinkling experiments at the experimental plots or by natural rainfall; the input rate is related to the total hillslope area.Time is given as hours since beginning of the discharge measurement at the cut-bank.
the times to first breakthrough and peak concentration, and the correlation of the breakthrough curves (Pearson's r).

Simulated and observed hillslope runoff
The model evaluation criteria of the 64 basic model set-ups are summarised in Table 3.Five runs fulfilled all three criteria (NSE, water balance error, surface runoff ratio).Acceptable matches of simulated and observed hydrographs with a NSE higher than 0.75 (maximum NSE = 0.86) were achieved by 22 of the basic runs, and 38 model set-ups matched the observed water balance within an error of ± 10 % (minimum error 1 %).A surface runoff ratio of less than 10 % of total runoff during the sprinkling phase was found for 24 model set-ups, while in 27 simulations surface runoff constituted more than 90 % of total outflow during the sprinkling phase.The hydrographs of the five successful basic model set-ups are displayed in Fig. 3a-e.The details of these model set-ups are summarised in Table 4.It is noticeable that all of these five set-ups involved the presence of vertical structures and constant depth to bedrock.In each set-up at least one lateral pathway was present, either the lateral structure or the soilbedrock interface, or both.A litter layer was not present in two of the five runs, but these set-ups had a higher number of vertical structures.
The importance of structures became apparent in comparison with the homogeneous set-ups of the additional model variants.In the null set-up without any structures ("h-01"), the entire hillslope outflow occurred as surface runoff.This set-up is thus to be rejected, because no significant amount of surface runoff had been observed in the field.Additionally, the NSE for this set-up is rather low (NSE = 0.31), as the simulated response during the rainfall phase is rising and falling much more abruptly than observed, although the hydrograph during the sprinkling phase is matched well (Fig. 3f).The four set-ups with a homogeneous, conductive soil mantle above the bedrock yielded only right boundary flux and not any surface runoff, but the hydrographs resulting from these uniform parameterizations were strongly damped and delayed compared to the observations (NSE between −134.8 and −71.4).
The additional model variants derived by modifications of the basic set-ups, which included widening of the lateral structures (litter layer, lateral pathway, soil-bedrock interface), increasing the density of vertical structures and increasing the hydraulic conductivity, did not improve the re-sults of the water flow simulation in terms of the selection criteria.Only one additional set-up (v-45) fulfilled all three criteria.Besides a lateral pathway and a litter layer, this setup comprised a widened interface above the variable bedrock topography.The widened interface brought about a significant reduction of the water balance error, while the NSE was inferior to the corresponding base case set-up (Table 4).We Table 5. Summary of solute transport characteristics of successful simulations and observation: Times to breakthrough and peak, recovery of tracer, correlation of observed and simulated breakthrough curves (Pearson's r), and the scaling factor f a (percentage of preferential flow paths across the slope width) that produced the maximum correlation in the solute transport simulations.The observed percentage of preferential flow paths is a mean value found with a plot-scale dye-staining experiment (Wienhöfer et al., 2009a) for depths greater 0.05 m. therefore focus in the following on the five acceptable basic set-ups.

Observed and simulated solute dynamics
Transport of solute through the subsurface to the hillslope toe was generally found in 51 of the 64 basic set-ups, and 24 of these also fulfilled the surface runoff criterion.In these cases, the bulk of the simulated solute transport occurred via the implemented preferential pathways.No solute transport to the hillslope toe was simulated with set-ups that either contained no structures, vertical structures without any lateral structures, or the soil-bedrock interface and/or the lateral pathway without any vertical structures.The solute transport simulations of the additional model set-ups essentially produced similar results.The first breakthrough of tracer in the experiments at 28.7 m distance along the slope line was recorded after only 0.77 h, and the peak concentration was reached after 2.00 h (Wienhöfer et al., 2009a).The initial solute transport simulations (f a = 18.75 %) of the five acceptable basic set-ups yielded breakthrough times between 1.42 and 1.67 h, and peak times between 2.08 and 3.42 h, respectively.Reduction of the scaling factor f a to 12.5-17.5% further improved the match of the solute breakthrough curves (Table 5); breakthrough times were between 1.17 and 1.58 h, and peak times between 1.92 and 3.25 h, respectively.Two set-ups matched the observed shape of the tracer breakthrough best ("c-07", "c09"; Fig. 4); in both cases, the correlation of simulated and observed was high (r = 0.94) with f a = 12.5 % (Table 5).
In the experiment, only a small fraction of 13 % of the mobile tracer mass was recovered by the end of the first experimental stage considered in this paper (see Sect. 2.1.2).Recovery was much higher for the majority of the simulations  5).Consequently, the height of the simulated and observed solute fluxes differed considerably (Fig. 4).
The observed maximum concentration in the outflow was 30.55 µg L −1 , which corresponds to a maximum transport rate of 3.12 µg s −1 , while the highest maximum transport rate in the five acceptable basic set-ups was 288.43 µg s −1 ("Run c-19", Fig. 4).Only set-ups that contained a soil-bedrock interface and no additional lateral pathway had a recovery rate of less than 20 % due to increased storage of solute within the soil matrix (e.g."Run c-20"; Table 5 and Fig. 4).

Discussion
The main objective of this study was to explore the modelling approach of representing preferential flow paths as distinct, connected elements in a 2-D numerical model.It expands on earlier studies (Klaus andZehe, 2010, 2011) by refining the approach and testing it on a different setting.Similar to the Klaus and Zehe studies, the modelling approach allowed successful simulations of water flows and solute transport at the hillslope scale.In the following we elaborate and expand on the specific experiences made with the simulation of water flows and solute transport in this work, and evaluate the advantages and limitations of the modelling approach with reference to the literature.

Simulation of preferential flow and hydrodynamic hillslope response
The modelling approach of representing preferential flow paths as distinct, connected elements of low flow resistivity was successful in several aspects.The approach allowed modelling the dominant processes of preferential infiltration into vertically oriented flow paths and subsequent preferential flow in laterally oriented structures, and the outflow hydrograph of the hillslope was matched satisfactorily.Especially well fitted were the height and the onset of the outflow in response to sprinkling, and the magnitude and timing of the major peaks in response to natural rainfall, although the models were not calibrated on peak heights or in any other way.It is particularly remarkable that set-ups which matched the observed response to the steady-state rainfall simulation on parts of the hillslope also matched the observed response to natural rainfall on the entire surface area (Fig. 1), because only the hillslope width was used for scaling the input between the two phases of the simulations.This suggests a certain predictive capability of these set-ups in conjunction with the modelling approach.
Of course, the match of simulated and observed hillslope outflow was not absolutely perfect.The simulations differed from the observations during the recession phases, the peak heights of the three major peaks during natural rainfall were not matched equally well, and the small peak after the first major peak was not modelled by any of the simulations (Fig. 3).Although these deviations from the observed hillslope outflow may appear rather small in light of the fact that we are modelling the complex system of natural hillslope with a high degree of heterogeneity and a perfect fit of the model would never had been expected, it is illustrative to discuss this topic in further detail.One possible explanation would be incorrect observations, which cannot be fully excluded during field experiments.To simplify matters, we assume that the observations reported by Wienhöfer et al. (2009a) depict the hillslope hydrology correctly within typical ranges of uncertainty, and that these are reflected in the chosen acceptance criteria.Other possible reasons for the mismatch of simulated and observed hillslope outflow could be due to the modelling approach in general or due to the specific implementation of the approach in this particular study.Limitations of the modelling approach are related to the conceptualisation of preferential flow paths as highly porous media and the process representation using the Darcy-Richards equation, as well as to the reduction of the three-dimensional hillslope to a two-dimensional cross section.These aspects are discussed in further detail below.But even if conceptualisation and process description were perfect, the imperfect knowledge about the system itself would still lead to considerable uncertainty in setting up and parameterizing a spatially explicit process model.The general lack of complete information on the internal build-up of a hydrological system basically makes it a "black box" with many degrees of freedom for the modeller.This black box might be illuminated at selected "grey spots" where field observations for constraining the model set-up are available, and has to be described by assumptions otherwise.
In our study we have tried to pursue this approach by keeping parameters fixed for which we had some data (e.g.soil matrix parameters, topographic gradient, and soil depth along the slope line), while we used other field evidence (e.g. from dye-staining experiments at the plot scale) to guide our conceptual model of preferential flow paths at the hillslope scale.As we did not have information on the exact arrangement of subsurface flow paths, we chose to generate different realisations corresponding to the conceptual model.The flow paths were modelled with a random component, but in a rather regular basic arrangement for better comparability.It cannot be ruled out that other and perhaps more irregular patterns would deliver comparable or better results.The set of basic model variants, however, seems to have been adequately reflecting the internal architecture of the hillslope; at least the tested modifications, for example limiting vertical flow paths to the upper half of the hillslope, did not improve the results.Other aspects for which hard information was lacking included, for example, the spatial variability of soil parameters, which was not accounted for in the model set-ups; heterogeneity was represented solely by the different types of structures.Measured soil hydraulic parameters were only at hand for the fine-grained soil matrix from a single location, whereas soil hydraulic parameters for preferential flow structures, bedrock and litter layer were chosen arbitrarily or from the literature.Variations of the hydraulic conductivity values in the additional model variants, however, did not yield better-fitting parameter combinations.Likewise, the assumption of spatially uniform rainfall input and canopy interception could be revisited, even if the influence of spatially variable throughfall on subsurface flow processes should only be secondary (Hopp and McDonnell, 2011;Bachmair and Weiler, 2012).All these assumptions could be replaced if site-specific data was available, which possibly, but not necessarily, might further improve the simulated hillslope response.
In our model, water flows are simulated using Richards' equation, and preferential flow pathways are conceptualised as an artificial porous medium with low flow resistivity and low retention capability.With application of this concept, we accept the trade-off between the possibility to incorporate preferential flow in distinctive structures into an existing numerical model, and possible errors resulting from the use of Richards' equation for water flows in these structure, which would rigorously have to be deemed inappropriate for describing flow and frictional losses in macropores (Beven and Germann, 1982).Despite this inconsistency, the concept has been proposed for representing macropore flow in single-domain (Nieber and Warner, 1991) and dual-domain (Gerke and van Genuchten, 1993) soil hydrological models.Especially with spatially explicit single-domain models, the implementation of distinctive structures is straightforward by choosing a respective parameterization for corresponding model regions as done in the present study, and this approach was successfully applied in modelling controlled experiments at the lab-scale (Castiglione et al., 2003;Lamy et al., 2009) and the plot-scale (Vogel et al., 2006;Nieber and Sidle, 2010).From this research it was concluded that exact flow rules are not the only concern (Lamy et al., 2009), and successful modelling of preferential flow is possible even with an approximate flow law such as the Richards equation if at the same time an approximate representation of structures is taken into account (Vogel et al., 2006).
In our application of the approach at the hillslope scale, we used a hypothetical network of vertical and lateral flow paths of limited spatial extent to conceptualise the structural heterogeneity of the hillslope observed at the plot-scale.These flow paths are not supposed to represent single structures spanning the entire hillslope, or structures of a single origin, but we rather hypothesise a network of connected flow paths constituted by several individual macropores, such as root holes, desiccation cracks and animal burrows, which are either connected directly to each other or via zones of higher porosity sustained by biological and/or hydrological processes.When these pathways are modelled as a highly porous medium, we implicitly include the surrounding matrix that might as well contribute to preferential flow (Lamy et al., 2009).Functional connectivity of individual macropores controlled by saturation state is implicitly modelled as well, as dry portions of the flow network will act as flow barriers in the simulation.
Another possible limitation is the simplification of the hillslope as a vertical 2-D cross section.The reduction to 2-D tends to underestimate connectivity compared to a 3-D realisation when treating heterogeneous porous media as a random field (Fiori and Jankovic, 2012).In our study this was far less a problem, as preferential flow paths were modelled explicitly and connectivity was hence prescribed a priori.The topology of the flow network in the 2-D model was represented sufficiently for the simulation of hydraulic response, not least because the study hillslope was much longer than its width, and the line of steepest descent in potential energy serves as symmetry axis.The simplification to 2-D, however, restricted the explicit representation of distinctive structures in the third dimension, which was unproblematic for the water flow simulations, but possibly caused inconsistencies in the solute transport as discussed in the following section.

Simulation of preferential flow and solute transport
The modelling approach allowed simulating solute transport via the preferential flow paths, and timing and shape of the observed breakthrough curve was matched well by several simulations.This corroborates that the flow velocity distribution was modelled acceptably, in which case adequate modelling of macroscopic dispersion is a by-product of the explicit consideration of distinctive structures (Vogel et al., 2006).Because the tracer was modelled as a conservative solute with low molecular diffusivity, solute transport was closely related to the water flow simulation.A similar reasoning as for the water flow simulation hence would apply for explaining minor discrepancies in the solute transport simulations, which could for instance result from the specific parameterization or the arrangement of structures.
In the context of the tracer simulations, also the observations and the implementation of the approach warrant a critical assessment.One peculiarity of the tracer observations was the low recovery rate, which could only partly be explained with irreversible sorption in the top soil that was observed in a column experiment (Wienhöfer et al., 2009a).Most of the solute transport simulations resulted in high recovery rates, while a low recovery was only found with some set-ups which featured vertical structures and a soil-bedrock interface, but no lateral pathway (e.g."c-20"; Tables 4 and  5).If timing and shape of the simulated breakthrough curves were closer to the observations, the lower recovery rate could possibly be considered an additional criterion for selecting these set-ups, provided that the underlying reasons for the apparent loss of tracer were also better understood.If a major part of the tracer was transported or stored in different flow paths, this would have to be reflected in different model set-ups with a more diverse arrangement of flow paths (which still would have to fulfil the criteria for water transport).
On the other hand, if the low recovery was due to sorption or decay processes, this should be accounted for in the model.While modelling decay and sorption principally is possible with CATFLOW, it could be required in this context to implement a more realistic splitting of the amount of solute that enters the preferential flow structures and the soil matrix.In the present 2-D approach the amount of tracer entering the structures is possibly overestimated, as the structures extend over the entire width of the hillslope in the third dimension.This is much more complicated to handle in the model than the incorporation of a scaling factor for the calculation of seepage velocities, and probably a dual-domain approach with a detailed treatment of the exchange processes between soil matrix and preferential flow paths or a 3-D approach would be needed.
In order to generally allow simulation of solute transport in the spatially explicit flow paths, a slight modification of the numerical tool CATFLOW was required.We found during preliminary tests that the internal interpolation of flow velocities for the random walk particle tracking led to a trapping of particles adjacent to the implemented structures.This phenomenon has been described earlier (LaBolle et al., 1996), and is particular serious in the case of finely resolved materials with highly contrasting properties as in our study.Turning off the interpolation of local flow velocities for the particle step kept the velocities contrasts between matrix and structures, and minimised unintentional overshoot of particles out of the structures and trapping of solute in the soil matrix.The change in the code only affected the simulation of solute transport and not the water flow calculation.

Equifinality and hydrological relevance of structures at the hillslope scale
Several different combinations of structural features were successful in simulating the hillslope hydrograph.This equifinality in structural set-ups was also reported in earlier studies (Weiler and McDonnell, 2007;Klaus andZehe, 2010, 2011).The occurrence of structural equifinality, however, is not necessarily a drawback of the concept, since it implies that the exact configuration of subsurface flow paths does not need to be known explicitly to simulate lateral preferential subsurface flow at the hillslope scale (Weiler and McDonnell, 2007).Equifinal model set-ups could also be used in an ensemble approach to assess the possible range of system behaviour in the light of uncertain model set-ups.
Testing different realisations will be required as long as our information on the subsurface flow network (or other structures at the hillslope scale) is incomplete, which will probably always be the case, even if all available evidence was used to set up a model and reduce the degrees of freedom.If complementary information is available that has not been used during model set-up, the equifinality can be further reduced, and hence the picture of the investigated system will become finer.Similar to Klaus and Zehe (2011), we were able to reduce the number of acceptable set-ups by evaluating the capability to simulate tracer transport.Another possibility to reduce equifinality a posteriori would be to test the different configurations for their long-term behaviour, if long-term data are available for comparison.
The number of equifinal set-ups will increase with the number of "similar" set-ups that are tested.Assessing the similarity of equifinal set-ups in contrast to unsuccessful setups can help in learning about possible configurations of the system under investigation.In the case of the study hillslope, only set-ups with a combination of interconnected lateral and vertical structures were successful.This corroborates our conceptual model based on plot-scale observations of preferential flow paths at the field site, and is in line with findings from earlier modelling studies that considered flow in distinctive structures in their models (Sidle et al., 2001;Jones and Connelly, 2002;Weiler and McDonnell, 2007).Of course, we cannot exclude that some configuration we have not tested would give similar or even better results, although we tried to cover quite a range of configurations constrained by the available field observations.The observed hillslope hydrograph, however, was not reproduced well by homogeneous set-ups without preferential pathways, which included different set-ups with a (conductive or less conductive) soil layer on top of bedrock, or with a (thin or thick) litter layer on top of soil matrix, which resemble some kind of layered soil profile.
In contrast to our expectation, the litter layer and the variable bedrock topography were not as important in the model.We would have expected a litter layer necessary to facilitate inflow into the vertical structures and thus the network of preferential pathways.Although this was the case, setups without a litter layer produced also acceptable results (Table 3).As the importance of bedrock topography has been highlighted by a number of studies, for example those made at the Panola research site (Freer et al., 2002;Trompvan Meerveld and McDonnell, 2006a;Hopp and McDonnell, 2009), we considered it necessary to measure the variable soil depth at our site and include it in our model.This allowed for comparison with a constant soil depth (i.e. a bedrock topography that resembles the surface topography).The results show that there is a difference in modelled hillslope response between these two bedrock configurations, but the set-ups with constant soil depth performed better.Of course, this finding strongly depends on the specific variable bedrock topography that was used in the models.A comprehensive analysis focussing on the effect of variable soil depth would have to investigate the influence of different interpolations and projections, or consider using a 3-D model, which is beyond the scope of this paper.We implemented one representation of the interpolated measurements that, however, captures the observed bedrock depressions (assuming that flows lateral to the slope line would take the route to the maximum depth), and which we consider adequate and representative for our purposes.
Set-ups where the variable bedrock topography exerted a strong influence did not match the observed response, for example the "homogeneous" set-ups with a variable soil depth and a conductive soil mantle, or set-ups with vertical flow paths and a conductive structure (soil-bedrock interface) along the variable bedrock topography.These set-ups were principally suitable to give rise to the 'fill-and-spill' mechanism (Tromp-van Meerveld and McDonnell, 2006a), when infiltrating water percolates down to the bedrock and builds a water table at the interface to the less permeable bedrock.The gradient of this water table, which in turn is determined by the interplay of percolating water and bedrock topography, then drives subsurface flow processes.This mechanism might be of greater importance at sites like Panola than at our study site.The Panola hillslope is less steep (13 • ; Tromp-van Meerveld and McDonnell, 2006b) than our hillslope (18 • to 54 • ), and the soil matrix is much more coarsetextured and permeable (saturated hydraulic conductivity K s = 1.8 × 10 −4 m s −1 ; Hopp and McDonnell, 2009) than at the Heumöser study site (K s = 1.77 × 10 −7 m s −1 ).This means that at Panola it is much more likely that a water table develops freely above the less permeable bedrock and "fills" bedrock depressions, before it "spills" downslope.
At the studied hillslope at Heumöser, infiltration and percolation are strongly bound to macroporous structures.If flow is directed into a network of pipes before it percolates to the bedrock, the geometry of the preferential flow network will determine the driving gradient much more than bedrock topography.This indicates that with the presence of structures above the bedrock, the role of the bedrock topography becomes secondary at our study hillslope (please see the Supplement for an example of the development of relative saturation in a simulation with vertical and lateral structures, litter layer, a soil-bedrock interface and variable bedrock topography).Generally, it thus is the geometry (topography) of the dominating structure that determines the water table gradient and in turn the flow response of a hillslope.This could be bedrock, but also a preferential flow network within the soil cover.
Variable bedrock topography can appear to have greater importance in the context of different modelling paradigms without spatially explicit consideration of structures.For example, Stadler et al. (2012) employed a 2-D dualpermeability model at the same site at Heumöser, and successfully simulated the hydraulic response to the first phase of the sprinkling experiment.They also found that flow predominantly occurred in the macropore domain, but as this was modelled spatially constant and isotropic within the soil mantle, the variable bedrock surface -the same as used in this study -exerted much more control on lateral preferential flows than in our model, in which lateral pathways were present.This is an illustrative example for the basic fact that possible model outcomes are generally bound to the inherent assumptions of the underlying perceptual and conceptual models.Modelling studies that investigated the role of variable bedrock topography and the 'fill-and-spill' mechanism have not considered preferential flow paths in their spatially explicit models used for virtual experiments at the hillslope scale (Weiler and McDonnell, 2004;James et al., 2010); Hopp and McDonnell (2009) even excluded pipe flow and pipe-flow observations from their modelling study of the Panola research site, although during the event that was chosen for calibration, pipe flow was contributing nearly half (45 %) of the observed hillslope outflow (Freer et al., 2002).Including the effects of distinctive preferential flow paths in future modelling studies that systematically analyse different controls on subsurface stormflow will help to continue learning from these virtual experiments using hillslope hydrological models.In our opinion, the results of the present study encourage further investigating the explicit consideration of distinctive structures in hydrological process models towards this aim.

Conclusions
This paper implements and evaluates the concept of incorporating preferential flow paths explicitly as distinctive structures in a process-based model.The approach was applied within the 2-D numerical model CATFLOW for modelling water flows and solute transport observed during a field experiment at a steep forested hillslope.The model successfully represented hillslope hydrological response by depicting the sharp contrast in flux density between structures and matrix, and the solute transport simulations matched timing and shape of the observed breakthrough curve well, indicating that macrodispersion induced by preferential flow was captured well by the topology of the preferential flow network.Furthermore, the tracer simulations could further reduce the small number of equifinal model set-ups from the hydrograph simulation.
The configurations of successful model set-ups suggest that preferential flow bound to structures is a first-order control on the hydrology of the study hillslope, whereas spatial variability of soil depth is secondary in presence of a preferential flow network.We employed an established numerical model as a virtual reality in the sense that we modelled the unknown instead of searching for a suitable process description in a completely controlled system.The results of this study show that not only the flow equations and the numerical implementation have to fit the processes to be modelled, and that this has to be checked critically, but that also perception and conceptualisation of the system play a decisive role in the modelling process.Possible model outcomes are always bound to the assumptions inherent in the underlying perceptual and conceptual models and limitations of the numerical tool, and it is therefore mandatory to use a model with adequate complexity to include a wide range of possible processes.In order to discriminate their role against other possible controls on hillslope hydrology, distinctive structures

Fig. 1 .
Fig. 1.Map of the study hillslope showing locations of field observations and experimental plots as well as simulation areas considered for model set-up.The inset shows the location of the study area within Europe.

Fig. 2 .
Fig. 2. Set-up of hillslope model in profile view; the insets show enlarged detail with different realisations of explicit structures: (a) setup with litter layer, vertical flow paths (2 m spacing), lateral flow path, and soil-bedrock interface layer on interpolated bedrock topography; (b) set-up with widely spaced vertical flow paths (4 m spacing) and soil-bedrock interface layer at constant soil depth over the entire profile; the inset is at the same scale as (a).

Fig. 3 .
Fig.3.Observed hydrographs and simulated total outflow of the five acceptable basic model set-ups (a-e) and of a set-up without any structures (f).The identifier and NSE are given in the upper left corner of each panel.The total outflow is the sum of surface and subsurface runoff.In case of the acceptable model set-ups (a-e) total outflow mainly was subsurface flow, while surface runoff was all dominant for the set-up without any structures (f).The rainfall input (g, h) was either provided by the sprinkling experiments at the experimental plots or by natural rainfall; the input rate is related to the total hillslope area.Time is given as hours since beginning of the discharge measurement at the cut-bank.
Fig. 4. (a-e) Solute breakthrough curves of the five basic set-ups with acceptable water flow simulation in comparison with the observation.Time is given as hours since tracer application; please note the different scales for simulated and observed solute fluxes.

Table 1 .
Overview on model set-up variants: five different types of structures have been considered in different realisations for generating the structural set-ups; all 64 combinations of the "basic model variants", plus selected modifications of these combinations given as "additional model variants", were tested for simulation.

Table 2 .
Hydraulic and transport parameter values used for different materials in the basic model variants.

Table 3 .
Results of the water transport simulations of the 64 basic set-up variants in relation to the three model evaluation criteria ("N": NSE greater than 0.75; "W": water balance error less than 10 %; "S": surface runoff during the sprinkling phase less than 10% of total outflow); letters indicate which criteria were fulfilled, and text in bold, bold+italic, and underlined indicate how many criteria were fulfilled.

Table 4 .
Structural features and evaluation criteria of selected model set-ups ("+" indicates presence of the respective structural feature): the runs c-07, c-09, c-11, c-19, and c-20 from the basic model set-ups and the run v-45 from the additional model set-ups were acceptable for the hydrograph simulations with a NSE greater than 0.75, a water balance error less than 10 %, and surface runoff less than 10 % of total outflow during the sprinkling phase (surface runoff ratio); the run v-21 from the basic model set-ups is given for comparison with the modified set-up v-45 with a widened soil-bedrock interface.