A novel explicit approach to model bromide and pesticide transport in soils containing macropores

Introduction Conclusions References


Introduction
Since the mid-nineties it has become evident that preferential flow is the rule rather than the exception in structured soils (Flury, 1996).Preferential flow processes are highly relevant in runoff generation at the hillslope (Lindenmaier et al., 2005;Weiler and McDonnell, 2007;Wienhöfer et al., 2009;Zehe and Sivapalan, 2009) and headwater scales (Zehe and Blöschl, 2004;Zehe et al., 2005;2006) and consequently are a prime cause for the spatial variability in soil water content at hillslope scale (De Lannoy et al., 2006;Zehe et al., 2010b).Originally, the term "preferential flow" was coined after realising that water flow and transport in soils containing noncapillary structures -often worm burrows, root channels or soil cracks -was much faster than could be expected from classical theory of flow and transport in porous media (Beven and Germann, 1982;Zehe and Flühler, 2001).Rapid flow in morphologically connected preferential flow paths controls transport and residence times of pesticides (e.g.Flury et al., 1995;Elliot et al., 2000;Šimůnek et al., 2003), anions (Flury et al., 1995;Villholth et al., 1998;Lennartz et al., 1999;Köhne et al., 2006;Stone and Wilson, 2006) and of strongly adsorbing phosphorus (Stamm et al., 1998;2002).Vertical flow distances through the unsaturated zone are in most cases too small for the central limit theorem to apply (Blöschl and Zehe, 2005), thus we cannot assume fully mixed conditions.Therefore it is not possible to describe the imperfect mixing by a well defined macro dispersion coefficient.
Consequently a variety of approaches has been proposed to effectively represent preferential flow in soil physical or soil hydrological models.Šimůnek et al. (2003), Gerke (2006) and Köhne et al. (2009a) published very nice and helpful reviews of different approaches, like composite hydraulic function models (Zurmühl and Durner, 1996), dual porosity and dual permeability models.Recent studies suggested explicit representation of preferential flow paths as connected structures in the model domain either at the pedon scale (Allaire et al., 2002a, b;Vogel et al., 2006;Sander and Gerke, 2009) or at hillslope scale (Klaus and Zehe, 2010;Zehe et al., 2010a).These studies showed promising results to reproduce water flow or solute transport.The number of studies dealing with modelling preferential transport of solutes at the hillslope and field scale is so large that an exhaustive overview is beyond the scope of this paper.Hendriks et al. (1999) simulated bromide and nitrate transport in a soil with apparent shrinkage cracks using a modified version of the FLOCR/ANIMO model.Model parameters were estimated by using data of groundwater levels, soil moisture, and concentrations of bromide in both the groundwater body and the soil profile.While the model results for nitrate were poor, bromide concentrations in soil could be adequately reproduced, and predicted bromide concentrations in groundwater matched the mean observed values.Gerke and Köhne (2004) modelled water flow and bromide transport at a tile drain field site using a dual porosity approach (based on the DUAL model of Gerke and van Genuchten (1993)) and a single domain approach.While the single domain approach failed to reproduce observed transport, the DUAL model was able to reproduce the long term dynamics of water and bromide transport well.Nevertheless, the authors report that the DUAL model failed in exactly reproducing observed bromide concentrations.Köhne et al. (2006) used seasonal discharge and bromide data observed at one plot to calibrate different models and then predicted seasonal transport behaviour at two nearby field plots.They compared a single porosity and a mobile-immobile model where both were set up in one and two dimensions (referred to as 1-D and 2-D hereafter).Calibration and prediction of bromide plateau concentrations were only successful when using the 2-D mobile-immobile model approach.Haws et al. (2005) compared a single-porosity and a dual-porosity approach based on HYDRUS-2-D ( Šimůnek et al., 1999) to simulate tile drain discharge and the transport of applied chloride and bromide.Although both model approaches were able to reproduce the shape of measured cumulative outflow well, simulated transport deviated considerably from the observations.The authors stated that no approach was able to reproduce rapid transport of solutes to the subsurface tile drain in an acceptable manner and concluded that a well matched hydrograph does not automatically imply that the corresponding model parameters are physical meaningful (see also McGuire et al., 2007).Köhne and Gerke (2005) compared a model based on 2-D Richards Equation and 2-D Convection-Dispersion-Equation (CDE) with a model based on the mobile-immobile approach when simulating tile drain discharge and bromide transport observed at an experimental field site in northern Germany.As only the mobile-immobile approach yielded a good match of the observed bromide peak, the authors concluded that this peak was caused by preferential transport.Gerke et al. (2007) achieved a clearly better match of tracer concentrations observed at the same site when using a 2-D-dual-permeability model and injecting the bromide tracer exclusively into the soil matrix domain of the model.Pang et al. (2000) studied bromide and pesticide transport at a tile drained field site in New Zealand using HYDRUS-2-D.Simulated soil moisture was in good accordance with observation.The model failed, however, in reproducing observed bromide and pesticide concentrations.Pang et al. (2000) thus concluded that a successful match of preferential transport requires a refinement of the dualpermeability approach.
Leaching and long term fate of non-conservative solutes such as pesticides and pharmaceutics is even far more difficult to predict as sorption behaviour and degradation processes come into play (Köhne et al., 2009b).Gärdenäs et al. (2006) compared four different approaches to simulate long term soil water flows and pesticide transport at a tile drained field site in southern Sweden.Transport behaviour was observed by taking 13 water samples during a period of 6 weeks.One approach was based on a modification of soil hydraulic conductivity near saturation and the remaining were a mobile-immobile, a dual-porosity and a dualpermeability approach, respectively.The dual-permeability approach performed best, followed by the dual-porosity approach, while the other approaches failed both to reproduce water flows and pesticide transport.Boivin et al. (2006) also used HYDRUS-2-D to simulate tile drain discharge and pesticide concentrations observed during a 100-day period for three different soils.Water samples were taken as flow proportional daily mixing samples.Simulated discharge was, after the soil hydraulic functions were modified to account for preferential flow, in good accordance with observations.Boivin et al. (2006) stress that the nonequilibrium mobile-immobile approach allowed acceptable reproduction of the observed pesticide concentrations, while the advection-dispersion approach failed.
The listed studies corroborate that it is most difficult to assess a model parameterisation that allows an acceptable match of both water flows and transport behavior, even when using a dual permeability approach.There is furthermore no straight forward link between parameters that characterise the preferential flow domain in dual permeability models and those parameters we can measure to characterise the density and depth distribution of a macropore network in natural field soils.In the present study we thus test whether an alternative approach, which may be parameterised on field observables, is feasible to simulate both water flows and transport behavior in structured field soils.The essence is to represent vertical and lateral preferential flow paths explicitly as morphologically connected paths of low flow resistance in the spatially highly resolved model domain.Klaus and Zehe (2010) gave evidence that this approach allows successful reproduction of tile drain event discharge recorded during an irrigation experiment at a tile drained field site in the Weiherbach catchment (Zehe and Flühler, 2001) (Beven, 1993) in their spatial model setup.Several "hillslope architectures" that were all consistent with the available extensive data base allowed a good reproduction of tile drain flow response (compare Sect. 3.1).In this follow-up study we use the 13 best hillslope architectures (effective Nash-Sutcliff ≥ 0.9) from Klaus and Zehe (2010) to address the following questions: -Does the approach proposed by Klaus and Zehe (2010) additionally allow reproduction of the bromide and the Isoproturon (IPU) concentrations that have been observed during the irrigation experiment at this site?
-Can we reduce the set of the best model setups (Klaus and Zehe, 2010) by rejecting those that do not match observed bromide transport behaviour?
-Can we employ those behavioural model setups that allowed a good match of both the discharge and bromide loss at the tile drain for simulating transport behaviour of the pesticide at this site?
In Sect. 2 we present a short characterisation of (a) the study site and the underlying irrigation experiment (Zehe and Flühler, 2001), (b) the corresponding modelling study by Klaus and Zehe (2010) and our model approach.In Sect. 3 we compare simulated and observed breakthroughs of bromide and IPU.The paper closes with discussion and conclusions in Sect. 4.

Study area, irrigation experiment and available data
The Weiherbach valley is a rural catchment of 3.6 km 2 size in south-west Germany.About 95 % of the total catchment area is used for growing agricultural crops.The climate is semi-humid with an average annual precipitation of 750-800 mm yr −1 , average annual runoff of 150 mm yr −1 , and annual potential evapo-transpiration of 775 mm yr −1 .The geology consists of Keuper sandstone, marl and mudstone (lower and middle Triassic) covered with a Loess layer up to 15 m thick.
The focus of this study is to simulate solute breakthrough of an irrigation experiment at a tile-drained field site that was performed in April 1997 by Zehe and Flühler (2001).One day before irrigation of the 900 m 2 site the herbicide IPU was applied by the farmer with a mounted spray bar and the site was instrumented with 25 TDR-Probes of 30 cm length.Irrigation was performed with twelve evenly distributed sprinklers in three time blocks with a sum of 46 ± 5 mm.Ten min-  Zehe and Flühler (2001).The upper panel presents the bromide and Isoproturon (IPU) concentrations, while the lower panel presents the tile drain discharge and the irrigation rate.
utes after onset of the first irrigation a bromide and a brilliant blue pulse were added to the sprinkling water.Water samples were taken at the tile drain outlet, which was located approximately 1.2 m below the field level.Tile drain discharge was measured using a calibrated Venturi tube and a pressure sensor.Zehe and Flühler (2001) report a first flush of high bromide and IPU concentrations in the early phase via the tile drains into the brook, while the tile discharge showed only a very small reaction.In contrast, later on during the second and third irrigation blocks, peak concentrations and peaks in tile drain discharge correlated well (Fig. 1).Additional dye tracer experiments showed that connected burrows of anecic worms that linked the surface and the gravel layer above the tile drain likely acted as short-cuts for water and solutes to bypass the soil matrix during the irrigation experiment.

The CATFLOW model
CATFLOW is a physically based model for plot, hillslope and catchment scale modelling of water and solute transport (Maurer, 1997;Zehe et al., 2001;Zehe and Blöschl, 2004;Zehe et al., 2005).A catchment is treated as an ensemble of hillslopes that are interconnected via the drainage network.In CATFLOW a hillslope is represented as a twodimensional cross section along the steepest descent line that is discretised by two-dimensional curvilinear orthogonal coordinates.The hillslope is thus assumed to be uniformly perpendicular to the slope line.Soil water dynamic is described by the Richards equation in the potential form that is numerically solved by an implicit mass conservative Picard iteration (Celia et al., 1990).The simulation time step is dynamically adjusted to achieve optimal convergence of www.hydrol-earth-syst-sci.net/15/2127/2011/ Hydrol.Earth Syst.Sci., 15, 2127-2144, 2011 the iteration scheme and ranges from several seconds during rainfall conditions up to one hour, which is the allowed maximum.Soil hydraulic functions are described after van Genuchten (1980) and Mualem (1976).Evapo-transpiration is represented by an advanced Soil Vegetation Atmosphere Transfer (SVAT) approach based on the Penman-Monteith equation.Surface runoff down the hillslopes is based on the diffusion wave approximation to the 1-D Saint-Venant equations.It is numerically solved by an explicit upstream finite difference scheme.Solute transport in CATFLOW is calculated by a particle tracking scheme based on a Random Walk approach.The deterministic part of a particle step is determined based on the seepage velocities using a backward two level Runge Kutta scheme as suggested by Roth and Hammel (1996).The random part of the particle step is determined based on the square root of the dispersion/diffusion coefficient multiplied by the time step.Sorption can be described either by linear instantaneous adsorption, the Freundlich or the Langmuir isotherms.Additionally, the model may account for first order decay.

Representation of preferential flow paths as connected structures in the model domain
The key objective of the present study is work out whether an explicit and realistic representation of preferential flow paths as morphologically connected structures of low flow resistance is feasible for modelling water flow and transport behavior in structured soils.This idea has been mainly motivated by the studies of Zehe and Blöschl (2004) and Vogel et al. (2006).Vertical macropores at our study site are dominated by anecic earth worm burrows but also the tile drain can be seen as a very simple example of a lateral flow pipe.As the interaction of vertical macropore flow and lateral pipe flow is of key importance in hillslope hydrology (Weiler and McDonnell, 2007), we decided to represent worm burrows and the tile drain in the same way in our model.This offers the additional advantage to resolve flow processes in downslope direction as well as in direction of the dominating subsurface structures.The population of worm burrows is generated by an agent based approach as explained in the following.

Agent based generation of worm burrows
For each surface element we simulated the number of occurring worm burrows at the soil surface by assuming that they are Poisson distributed (Beven and Clarke, 1986).The Poisson parameter is the average worm burrow density per unit area multiplied by the area of a surface element.Next we simulated the length of each macropore located in the surface element (can be more than one) assuming a normal distribution.In the next step we used a stochastic agent, which extended the worm burrow into the domain by an anisotropic random walk.The probability for digging into the vertical was 0.9 while the probability for digging in lateral direction (upslope and downslope) was 0.05 each.This approach ensures that the burrow is continuous and thus connected and allows accounting for tortuosity.When used with a grid as fine as in this study (see Sect. 2.4), this procedure yields a worm burrow system that is consistent with our findings (Zehe et al., 2010a) and those reported by Shipitalo and Butt (1999).After the procedure was completed we assigned a flag or characteristic function to any grid cell in the domain that is one if the grid cell contains at least one macropore and zero otherwise (Fig. 2).
Our fundamental assumptions are (a) that vertical water flow in the grid elements with worm burrows is dominated by macropore flow and (b) this can be accounted for by assigning an effective medium to these model elements.To this end we multiplied the number of macropores in a grid element by the maximum volumetric water flow in the worm burrow and divided this value by the cross sectional area of the grid element to obtain its effective hydraulic conductivity.Klaus and Zehe (2010) compared values for maximum volumetric water flows in worm burrows from different experimental studies of Shipitalo and Butt (1999), Weiler (2001), and Zehe and Flühler (2001) for this purpose.Finally, we assigned an effective porous medium with very low water retention to the model elements containing macropores to minimise the effect of capillary forces on the simulation results.The study of Klaus and Zehe (2010) revealed that this approach is feasible at least for simulating soil water flows, as this explicit treatment of macropores as connected structures yielded successful predictions of observed tile drain discharge with a Nash-Sutcliffe (NS) of up to 0.95.The big drawback of this approach is the relatively high computational cost.

Infiltration into connected structures
Infiltration and runoff generation is treated in CATFLOW by means of a Cauchy boundary condition.The uppermost model layer is usually very fine (1-5 cm).Rainfall is treated as a flow boundary condition, as long as the uppermost grid Hydrol.Earth Syst.Sci., 15,[2127][2128][2129][2130][2131][2132][2133][2134][2135][2136][2137][2138][2139][2140][2141][2142][2143][2144]2011 www.hydrol-earth-syst-sci.net/15/2127/2011/ cell is unsaturated.In the case of saturation, the model switches to a Dirichlet boundary condition using the overland flow depth as a pressure boundary to drive infiltration.The infiltrating water flux is then calculated using Darcy's law.This approach works well to capture the onset of overland flow as shown, for instance, in Zehe et al. (2001) or to simulate infiltration at small scales (Zehe and Blöschl, 2004).Infiltration into the grid cells that contain the effective macropore medium works identically.In the case of overland flow, the effective drainage area of the macropore elements may increase due to runoff from upslope sectors.The boundary conditions at the remaining boundaries were no flow at the upslope boundary, free flow (water arriving the boundary is leaving the system) at the downslope boundary and gravity flow at the lower boundary.

Spatial model setup
Ideally the outlined approach is used at a very fine spatial grid size in downslope/lateral direction (2-3 cm).In this case each surface element of the 2-D domain contains either a single worm burrow or not (Zehe et al., 2010a).Such a fine lateral discretisation of the 30 m long hillslope, combined with a vertical grid size of 2 cm, resulted however in computation times of 200 h for a single model run.This was unacceptable as the purpose of the previous model study of Klaus and Zehe (2010) was to investigate whether the inherent uncertainty of key model parameters results in equifinality of the spatial model setup.Thus this study used a lateral grid size of 30 cm and a vertical grid size of 2 cm.Klaus and Zehe (2010) used available extensive data on soil hydraulic properties and initial conditions as well as statistical data on the worm burrow system (the density per unit area and average burrow length) to generate 432 setups of the tile drained field site.Tile drain flow response to the irrigation was simulated for each of these setups using CATFLOW.To keep computation time reasonable, Klaus and Zehe (2010) used a downslope/lateral grid size of 30 cm, which resulted in a computation time of up to 20 h per run.
As outlined above Klaus and Zehe (2010) found that several spatial setups allowed a very good match of the observed hydrograph at the study site.None of these possible setups could be rejected as all key parameters (the worm burrow density, the hydraulic conductivity of the worm burrows, etc.) were within ranges either reported in the literature (Shipitalo and Butt, 1999;Weiler, 2001;Zehe and Flühler, 2001) or of local observations.Also, simulated soil moisture patterns within the model domain appeared to be realistic.
Within the present study we employ the 13 best spatial model setups, simulate in addition transport of bromide and compare the results with the corresponding concentrations observed during the irrigation test.Table 1 characterises the soil profile at the field site and lists the corresponding soil water characteristics and characterises the macropore medium.Table 2 lists the most important characteristics of the 13 best model runs, the corresponding Nash-Sutcliffe efficiency and the differences between simulated and observed time to peak (denoted as timing error hereafter).

Representation of the tile drain and weighting factor for scaling length specific outflow
The tile drain and its gravel embedding is, as explained above, also represented as a connected flow path by assigning a band with the "worm burrow" medium to a depth of 95-120 cm that extended across the entire length of the domain.Klaus and Zehe (2010) compared three different hydraulic conductivities of the tile drain.Table 2 lists the corresponding value for each of the most likely model setups.
As CATFLOW is a 2-D model the output is length specific and must be scaled by multiplying simulated length specific outflow by the width of cross section that is drained by the tile drain.Consequently, Klaus and Zehe (2010) scaled the output such that simulated peak flows matched the observed peak flows.The study revealed that the "drained width" of the cross section was in most cases between 1 and 3 m (compare Table 2).Please note that within the present study these widths were not modified to match bromide concentrations and cumulated bromide loss from the tile drain.This is thus a generic test to find out whether the approach proposed by Klaus and Zehe (2010) allows a consistent prediction of flow and transport.For further details about the limitations of the approach the reader should refer to the section "Key assumptions underlying the 2-D approach" and the discussion section in Klaus and Zehe (2010).

Transport parameters
Bromide transport was simulated for the 13 model setups using the initial soil moisture state listed in Table 2. Simulated irrigation occurred as during the field experiment in three blocks of uniform intensity with a sum of 46 ± 5 mm.At the beginning of simulation 1500 g of bromide was evenly distributed within the upper layer (2 cm) of the model domain.
During all simulated cases we used a constant isotropic effective diffusion coefficient of D = 5 × 10 −7 m 2 s −1 for bromide (Zehe and Blöschl, 2004).This selection is corroborated by observations from nearby tracer profiles that were obtained under a similar irrigation conditions and similar a soil.Klaus (2011) determined a dispersion coefficient D of 5.4 × 10 −7 m 2 s −1 for a travel time of 24 h and a D of 7.5 × 10 −8 m 2 s −1 after 48h travel time.As we mainly intended to test the feasibility of the explicit representation of preferential flow paths we do not present a sensitivity study of the dispersion coefficient on the simulation results.IPU transport was simulated for the model setup named "run 1" that yielded the best discharge simulation and an acceptable match of the time series of cumulated bromide mass flow at the tile drain outlet (see Sect. 3.2).Due to the high computational cost we restricted the IPU transport modelling exercise solely to this spatial model setup.As initial state 270 g of IPU were spatially homogeneously distributed in the upper 2 cm of the soil corresponding to the amount of IPU that has been applied to the field one day before the experiment.This considers the redistribution of the IPU that was applied on the soil surface the day before irrigation and likely entered the top soil by means of diffusion.We used the same effective diffusion coefficient as for bromide transport.Non-linear adsorption behaviour of IPU was parameterised according to the Freundlich-isotherm: Where C a is the concentration in the adsorped phase [g g C is solute concentration in the water phase [gl −1 ] and β is the Freundlich exponent [-].
As a first step we assumed that adsorption was uniform in the entire domain and selected the parameter set that performed best.Next we assigned these parameters exclusively to the macropore medium and characterised adsorption behaviour in the soil matrix be means of average parameters β=0.8 and k f =10 l β g −1 g 1−β .

Simulated and observed hydrographs
Figures 3, 4, and 5 present simulated discharge plotted against observed tile drain discharge (panels in the left column).Simulated discharge is generally in very good accordance with the observations, except that the first discharge increase comes systematically too late during simulations.While Nash-Sutcliffe efficiencies are larger than 0.9 (compare Tables 2 and 3), some of the simulations (runs 7, 8, 11, and 13) have water mass balance errors (deviation between modeled and observed total discharge volume) that exceed 10 %.Peak time errors are in general of the order of 7 or 17 min, corresponding to 1 or 2 model output time steps.Please note that discharge measurements during the irrigation experiment stopped 5 h after irrigation, although tile drain discharge was still approximately two times above the pre-event level.For additional details on the entire Monte Carlo study, especially on the sensitivity of predicted discharge on key parameters such as average density of worm burrows or maximum flow rates in macropores, please refer to Klaus and Zehe (2010).

Simulated and observed bromide concentration and mass flows
Figures 3, 4, and 5 present the time series of simulated and observed bromide concentrations (in the middle column) and cumulated bromide loss from the 13 simulations plotted against the corresponding observation (right column).As can be seen from Fig. 1 or the panels in the middle column of Figs. 3, 4, and 5, observed bromide concentration peaked in the very first sample taken 36 min after irrigation onset.This first flush can be attributed to rapid preferential transport as no significant increase in discharge was detectable.Within the second irrigation block, 145 min after irrigation onset, bromide concentration levelled out at a concentration of 1.9 mg l −1 for almost 40 min, then declined to a pretty constant value of around 1.2 mg l −1 .
Simulated bromide concentrations are calculated from the cumulated bromide loss within a model output time step, which is 10 min, divided by the accumulated water volume that left the model domain within this time.Simulated bromide concentrations are clearly more "noisy" than the observed ones.Although the simulated bromide concentrations are different among model runs, the temporal pattern is similar for all runs.Simulated concentrations show an early peak in the order of 0.5 to 2 mg l −1 within the first 30 min, then decrease to zero or near zero values for all runs and rerise again to values of about 2 mg l −1 .The observed peak at minute 36 is often reproduced as a double peak by the model.During the second irrigation block some simulations reproduce a concentration peak instead of a plateau concentration, while some runs capture this plateau.Simulated concentrations are in acceptable accordance with the observed ones during the last irrigation block.Overall, we can state that simulated concentrations match the magnitude of the observed ones and also roughly the timing of the peaks.Part of the mismatch might be attributed to the fact that rather small errors in simulated discharge (magnitude and timing) propagate to considerable errors in concentrations (magnitude and timing).However, as the model was not tuned to reproduce the observed concentration we did not expect a really good match in absolute concentrations on a time scale of minutes.There is however a good accordance between the cumulated bromide losses derived from simulation when compared to the observed cumulated bromide loss (Figs. 3,4,and 5,right column).Despite the fact that all simulations underestimate the observed cumulated leaching within the early phase of the irrigation, simulation and observations are in very good linear accordance, as the coefficients of determination (not time corrected), R 2 , are always larger than 0.95.Within the model runs 1, 4, 9, and 10, the simulation matches the observed accumulated bromide loss to the tile drain with a relative error in the order of 5-12 %, which is within the error range of the observed data, 10 % for discharge and 15 % for bromide (compare Table 3).Table 3 lists furthermore the slope and intersect of a linear regression between observed and simulated cumulative bromide leaching; regression lines of model runs 1, 4, 9, and 10 are close to one.
Thus we may state that the spatial model setups of 1, 4, 9, and 10 reproduce the observed time series of cumulated bromide loss in an almost unbiased manner.This is remarkable as the model has not been attuned for this purpose.Klaus and Zehe (2010) showed that the soil moisture dynamics followed a realistic pattern.Figure 6 presents the hillslope scale bromide concentrations at five different time points.The simulated bromide concentration pattern is apparently controlled by the macropore system as patches of higher bromide concentration cluster along the border of the preferential flow paths.

Simulated and observed cumulative leaching of IPU
We selected the model setup one (run 1) for simulating IPU transport as it performed best with respect to discharge and was ranked fourth with respect to reproducing cumulated bromide loss.Simulated and observed cumulated leaching of IPU is shown in Fig. 7, each panel corresponds to a case were β was fixed at 0.2, 0.5, 0.8 or 1 and k f was varied between 0.26 and 20 l β g −1 g 1−β .To put it in short, none of the simulations based on these parameterisations came anywhere close to the observed leaching behaviour.Retardation was clearly far too strong.This finding is in accordance with observations during the irrigation experiment (Zehe and Flühler, 2001) who pointed out that the effective retardation coefficient of IPU against bromide was almost one.
We therefore reduced k f as explained in Sect.2.4.2 to stepwise reduce retardation coefficients to one and compared simulated and observed cumulated leaching of IPU (Fig. 8).The lowest k f -values of 1 × 10 −6 and 1 × 10 −7 l β g −1 g 1−β yield a good match of the early non-linear rise of the cumulative leaching curve regardless which value is used for β.However, in total they lead to a strong overestimation of the total IPU mass that leached into the tile drain.Figure 8 also shows that a combination of β = 0.2 and k f = 1 × 10 −5 l β g −1 g 1−β or β = 0.8 and k f = 1 × 10 −4 l β g −1 g 1−β allow an acceptable estimate of the total mass of IPU that leached into the tile drains.This suggests that there are several combinations of β and k f -values that allow either a good match of cumulated IPU leaching in the early stages or of total IPU mass at the end of the experiment.However there was no parameterisation that yielded a good match of both, the rapid increase at the beginning of the experiment and the total observed transport.We thus tested whether different adsorption characteristics in the soil matrix and in the macropore medium would improve the model performance.We selected to compare runs with a β-value of 0.8 for the soil matrix and the macropores and varied the k f -values in the macropores (1×10 −3 , 1×10 −4 , 1×10 −5 , 1×10 −6 l β g −1 g 1−β ), while the k f -values Hydrol.Earth Syst.Sci., 15,[2127][2128][2129][2130][2131][2132][2133][2134][2135][2136][2137][2138][2139][2140][2141][2142][2143][2144]2011 www.hydrol-earth-syst-sci.net/15/2127/2011/ in the soil matrix was kept constant (10 l β g −1 g 1−β ).Heterogeneous adsorption caused a slight "damping" of the cumulated leaching curves, but did not improve the shape of the modelled breakthrough curve (Fig. 9).We may thus state that without any prior knowledge about the retardation behaviour, simulated leaching is way off the observations.When such information is available, as in the present case, the proposed model approach can be attuned to reproduce the overall amount that leached to the tile drain, while the dynamics is hardly reproduced.Including non-equilibrium sorption into simulation could allow a better match of the observed leaching dynamics.

"Preferential structures" for preferential transport
The proposed approach to represent vertical macroporesin our case anecic earthworm burrow -and the tile drain as morphologically connected structures in a 2-D model domain (Klaus and Zehe, 2010;Zehe et al., 2010a) is feasible to: -Setup several spatial hillslope architectures that are consistent with detailed field observations, especially on the density and depth distribution of worm burrows and to reproduce observed tile drain discharge (see also Klaus and Zehe, 2010).This is consistent with findings at even smaller scales i.e. column experiments of Comegna et al. (2001) and Abasi et al. (2003).
-The 13 best hillslope architectures, which achieved a Nash-Sutcliffe efficiency larger than 0.9 in combination with a small peak time error, predicted dynamics of accumulated bromide leaching into the tile drain well (R 2 >0.95) without additional calibration.
Four of those hillslope architectures (run 1, 4, 9, and 10) matched the total accumulated water balance (compare Table 3) within the relative error range of 10 % according to Zehe and Flühler (2001).It is remarkable that the same hillslope architectures also matched the total amount of bromide that leached into the tile drain within the range of the observation errors (please note that due to Gaussian error propagation, the relative error of the bromide mass flow at the tile drain outlet equals the sum of the relative errors in discharge (10 %) and bromide concentrations (5 %)).We may thus state that only four of the 432 hillslope architectures tested in the study of Klaus and Zehe (2010) matched both flow dynamics and bromide leaching very well and in an unbiased manner.
However, the proposed model structures are, without additional calibration, unable to reproduce exactly the time series of bromide concentrations, but do capture the order of magnitude of the concentrations.Bromide concentrations and mass flows in the first phase of the irrigation experiment are systematically underestimated.A closer look at model runs 1, 4, 9, and 10 reveals that predicted concentrations match the observations in later stages (t >100 min) in an acceptable manner.Run 4 performs best with respect to the overall bromide mass balance, and run 1 captures even the magnitude of the maximum observed bromide concentration, although the timing is roughly 50 min too late.This is consistent with the fact that first discharge reactions within the model are also systematically 50 min too late (except for run 9 and 10).
Transport velocities and bromide concentration in the model are obviously both too small to match the first flush of bromide into the tile drain.One possible explanation for this is that the proposed approach is too slow to capture the obviously very fast initialisation of macropore flow in the real system.This might be true, but note that simulated first arrivals of bromide are around 10 min (compare Figs. 3, 4,  and 5), which is deemed to be rather fast.Another explanation for the mismatch of early concentrations is that the proposed approach does not sufficiently account for bromide exchange between the surrounding soil matrix and the connected structures (Weiler and Flühler, 2003;Coppola et al., 2009) -although exchange happens as indicated by Fig. 6.The proposed approach can certainly neither account for processes like solute exclusion nor for the possibly hydrophobic behavior of the worm burrow coating.An increased model complexity for instance by including mobile-immobile water would certainly add more flexibility to the model to closer reproduce the bromide concentrations.The drawback is however to include additional parameters that are not directly observable in field soils.This would thus also result in increased equifinality.
Furthermore, an explanation for the mismatch of early concentrations is that we simulate a spatially homogeneous irrigation, while Zehe and Flühler (2001) reported that their irrigation scheme was rather heterogeneous shortly after bromide application: a lot of irrigation water with a large amount of bromide concentrated along the middle axis of the 30 by 30 m 2 large irrigation site (compare their Fig.12).These locally increased irrigation and mass inputs could, according to Weiler and Naef (2003), cause fast preferential flow rates and enhance mass transport into the tile drain and result in the observed first flush.This effect cannot be reproduced by simulating (a) a homogeneous irrigation and (b) neglecting the microtopography that determines local scale lateral redistribution of ponded irrigation water.We have neither information about the exact spatial pattern of irrigation rates nor details on microtopography (which would be rather difficult to include into a 2-D model).We thus omit any trial to reproduce this first flush effect within a simulation, as the underlying assumption would be rather speculative.
As CATFLOW is a 2-D model, predicted discharge is length specific and model output has thus to be scaled to 3-D reality by specifying the width of the irrigated site that feeds tile drain discharge.As elaborated in Klaus and Zehe (2010) we selected this drained width such that the simulated peak discharge matched the observed peak discharge.As these values, which are in the order of 1-3 m (compare Table 2), were kept constant in the present study, we may state that the same scaling factor allows simulation of flow and transport in a consistent manner.This corroborates that this factor is more than a "quick fix" to reproduce water flows, but is related to the width of the irrigation site that contributed to the tile drain water flow and bromide mass flows and integrates different processes that lead to transport perpendicular to the direction of the tile drain.This does not mean that the scaling factor matches the contributing width of the hillslope in an exact manner.
We thus overall conclude that: -A realistic representation of dominating structures and their topology enables the prediction of preferential water and mass flows at tile drained hillslopes.This requires detailed knowledge about the density and topology of the preferential flow paths, in our case anecic worm burrows, that are structures that emerge at the hillslope scale.They are a fingerprint of population dynamics and behavior of earthworms as ecosystem engineers and obey ecological rather than physical principles.Thus the necessary realistic representation is needed to describe the structures and no upscaling from microscale physical principles is possible to obtain the necessary information.
-Those hillslope architectures that perform best with respect to the integral discharge response are also suitable for predicting dynamics of accumulated bromide leaching, partly without a bias.This finding is clearly opposite to what has been achieved by Haws et al. (2005) but in accordance with the finding of McGuire et al. (2007).
The first point implies that we had to be very open to learn from soil ecologists.The major outcome of this learning process is that the generated structures match geometry, tortuosity and topology of worm burrows well when used at a finer grid resolution (Zehe et al., 2010a).The necessary data for generation of these burrow systems furthermore can be directly observed in the field and on the long term hopefully inferred from species distributions models from soil ecology.

Clear deficiencies to reproduce observed pesticides leaching
Even if the model structure has been shown to reproduce both discharge and bromide transport behaviour, the reproduction of the observed accumulated leaching for IPU was rather challenging.A parameterisation of the Freundlich isotherm based on values reported in the Footprint data base (www.eufootprint.org)did not reveal any meaningful result.Based on the finding of Zehe and Flühler (2001) at this site, who report an effective retardation coefficient of one at the end to the experiment and similar findings of Kung et al. (2000) who found similar behaviour for bromide and the sorbing tracer rhodamin WT, we reduced the k f values and thus the retardation coefficient stepwise to one.With this, the proposed model approach could be attuned to reproduce the overall IPU mass that leached to the tile drain.Temporal dynamics of cumulated leaching was hardly reproduced.The strongly reduced retardation of reactive solutes is likely caused by the differences in the chemical and physical properties of the organic coating of the macropore walls compared to bulk soil matrix.While adsorption capacity of macropore walls can be stronger than in bulk soil material, the overall retardation may still be less, due to the smaller surface area per volume in macropores or rate limited adsorption (e.g.Jarvis, 2007).As travel times in macropores are likely small compared to the time scale of lateral diffusive mixing equilibrium approaches are deemed to be too simple to describe adsorption behavior.A more complex kinetic approach may lead to an improvement match of the temporal dynamics.A heterogeneous parameterisation of adsorption parameter between the soil matrix and the connected flow paths (e.g. as suggested by Ray et al., 2004) did not lead to a clear improvement.The retardation characteristics of the soil at the study site are probably more heterogeneous, both in the soil matrix and the soil macropores, when considering effects like macropore coating (e.g.Gerke, 2006).An alternative explanation for the differences between the modelled and observed IPU breakthrough is that the first flush of pesticides may be explained by facilitated transport, as suggested by Zehe and Flühler (2001) in their experiment.The IPU concentrations of their experiment included both facilitated and dissolved transport.The facilitated transport behaviour has been as documented in several studies (de Jonge et al., 1998(de Jonge et al., , 2004;;Villholt et al., 2000).Facilitated transport may still be active during the later stages of the experiment.This process cannot be captured with the proposed model.Overall, we may state that a simulation of pesticide transport, without having detailed local data on the adsorption characteristics in different subsurface compartments (e.g.Mallawatantri et al., 1996;Larsbo et al., 2009) and a much more complex description of the adsorption process, is largely an unsatisfying issue of trial-and-error, even when using data from the Footprint data base.

Equifinality of spatial hillslope structures
Several networks of worm burrows allow us to produce an integrated system response, i.e. dynamics of water flow and of accumulated bromide leaching, although the set of best model structures can be boiled down to four (runs 1, 4, 9, and 10) by taking the water and the bromide mass balances into account.Also, McGuire et al. (2007) were able to reduce model uncertainty and increased parameter identifiability by the additional use of tracer data in the hillslope model used (Hill-vi, Weiler and McDonnell, 2004), compared to a single calibration on discharge data.Figure 10 shows the relative difference between simulated and observed cumulated bromide leaching plotted against the Nash-Sutcliffe efficiency (left panel) and against the water balance error (right panel).As there is no clear relationship between these quality criteria, all of them have to be evaluated to reject non behavioural model setups.
To further elaborate the similarities and differences of the four acceptable model setups (see Table 2) we state that all those assumed wet initial conditions for the entire model domain (Klaus and Zehe, 2010).This is consistent with available soil moisture data from 25 TDR (Zehe and Flühler, 2001), which characterised top soil water content.This corroborates that detailed data on the initial soil moisture state reduces equifinality.All of the acceptable model setups are parameterised with the highest possible hydraulic conductivity for the tile drain (2.03×10 −1 m s −1 ) and an impermeable layer below the tile drain.The maximum water flow in macropores and the density of worm burrows per unit area show a clear interaction.Runs 1 and 4 are based on 30 worm burrows per m 2 and a lower maximum water flow of 7.9×10 −6 m 3 s −1 in macropores, while the runs 9 and 10 are parameterised with a reduced worm burrow density (10 per m 2 ) which is compensated by a higher maximum water flow of 1.33×10 −5 m 3 s −1 .The occurrence of saturated condi-tions below the tile drain showed no clear pattern between the different runs.
Prior to this study, we were not sure whether the proposed way of treating macropores and scaling the 2-D model domain is feasible at all.We thus did not select criteria for model rejection a priori, as it is recommended by Beven (2010) for hypothesis testing.However, based on the now available knowledge, we propose the following twostage procedure for rejecting model structures: -Select those hillslope architectures as behavioural that reproduced tile drain discharge with a Nash-Sutcliffe efficiency larger than 0.9 and match the water balance within the error range (10 %).Thus runs 7, 8, 11, and 13 would be rejected.
-Select from the behavioural hillslope architectures of the previous stage those that match dynamics of accumulated bromide leaching with a coefficient of determination larger than 0.9 and match the bromide mass balance within the error range (15 %).Thus runs 2, 3,5,6,8,11,12, and 13 would be rejected.
From this we conclude: -The number of hillslope architectures that are behavioural with respect to reproducing tile drain discharge dynamics can be reduced by using orthogonal data sources such as tracer data (McGuire et al., 2007;Klaus et al., 2008) and analysing the mass balances.
-Generic knowledge about the origin of dominating structures is crucial to reproduce these structures in a simplified yet a sufficiently realistic manner and thus reduce equifinality in the spatial model setup to a minimum amount.
Nevertheless, in our case four behavioural hillslope architecture structures remain.This is consistent with findings of Comegna et al. (2001), who showed that a good accordance of observed and modelled hydrographs does not lead to an appropriate unique description of the solute transport mechanisms.
Probably we have to accept that a possible cause of equifinality is that part of the indeterminacy may be essentially system inherent, in the sense that several types of architectures of subsurface flow pathways may yield the same integral response in discharge and mass flows.

Fig. 1 .
Fig. 1.The experimental results from the study ofZehe and Flühler (2001).The upper panel presents the bromide and Isoproturon (IPU) concentrations, while the lower panel presents the tile drain discharge and the irrigation rate.

Fig. 2 .
Fig. 2. One of the model setups.Hillslope with different soil layers and the effective macropore region as explicit structures (red).

Fig. 3 .
Fig. 3.The left column presents the observed discharge (solid line) and the modelled discharge (dotted line), the centre column presents the observed (solid line) and modelled (dotted line) bromide concentrations in the tile drain, while the right column presents the cumulated modelled bromide outflow versus the measured one and a linear regression line between them (b denotes the intersect and m the slope of the regression).The runs 1 to 5 are shown here.

Fig. 4 .
Fig. 4. The left column presents the observed discharge (solid line) and the modelled discharge (dotted line), the centre column presents the observed (solid line) and modelled (dotted line) bromide concentrations in the tile drain, while the right column presents the cumulated modelled bromide outflow versus the measured one and a linear regression line between them (b denotes the intersect and m the slope of the regression).The runs 6 to 10 are shown here.

Fig. 5 .
Fig. 5.The left column presents the observed discharge (solid line) and the modelled discharge (dotted line), the centre column presents the observed (solid line) and modelled (dotted line) bromide concentrations in the tile drain, while the right column presents the cumulated modelled bromide outflow versus the measured one and a linear regression line between them (b denotes the intersect and m the slope of the regression).The runs 11 to 13 are shown here.

Fig. 7 .
Fig. 7. Cumulated IPU loss plotted against simulation time for different parameterisations of the Freundlich isotherm.The solid line presents the observed IPU leaching, the dashed lines model runs with different parameterisation of the Freundlich parameters.Every plot includes eight model realisations with a fixed value for the Freundlich exponent β and varied values for the Freundlich coefficient k f between 0.26 and 20 l β g −1 g 1−β (lowest line) within 9 steps.

Fig. 8 .
Fig. 8. Cumulated IPU loss plotted against simulation time for different parameterisations of the Freundlich isotherm.The solid line presents the observed IPU leaching, the dashed lines model runs with different parameterisation of the Freundlich parameters.Every plot includes eight model realisations with a fixed value for the Freundlich exponent β and varied values for the Freundlich coefficient k f between 0.1 (lowest line) and 1×10 −7 l β g −1 g 1−β within 8 steps.

Fig. 9 .
Fig. 9. Modelling IPU transport with spatial homogenous (left) and spatial heterogeneous (right) sorption parameter.The solid black line denotes the sampled cumulative breakthrough curve, the dashed lines display the modelled breakthrough curves.

Fig. 10 .
Fig.10.Left panel: differences in bromide (observed to modelled) versus the Nash-Sutcliffe efficiency of the discharge modelling.Right panel: differences in total modelled bromide versus total observed bromide concentrations differences in total modelled to total observed discharge.

Table 2 .
Parameters and goodness of fit criteria to characterise the 13 model setups.NS: is the Nash-Sutcliffe efficiency, DT: differences between modelled and simulated time to peak, NP: number of macropores per unit area [m −2 ], MWF: maximum volumetric water flow in a macropore [m 3 s −1 ], IMC: initial moisture conditions (the wet case corresponds to a suction head of −50.66 hPa throughout the profile; the intermediate case had a suction head of −303.98 hPa in the upper 40 cm of the profile and −101.33 hPa below this layer; the dry setting corresponds to a head of −405.3 hPa in the upper 40 cm of the profile and −151.99 hPa in the remaining soil profile), k D : tile drain conductivity [cm s −1 ], GW: Groundwater occurrence, C-Sec: drained width of the cross-section in meters adjusted to match the peak.

Table 3 .
Number of model run, Nash-Sutcliffe efficiency (NS) of hydrograph modelling, relative water balance error in % (WB), relative error in bromide mass balance in % (BMB), the slope of the regression line fitted to the compared cumulative transport of measured and modelled values (m), the intercept of the regression line (b) and the corresponding coefficient of determination (R 2 ).