Journal of Contaminant Hydrology 104 (2009) 4–35

Contents lists available at ScienceDirect

Journal of Contaminant Hydrology
j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j c o n h yd

A review of model applications for structured soils: a) Water ﬂow andtracer transport
John Maximilian Köhne a,⁎, Sigrid Köhne b, Jirka Šimůnek c,1aHelmholtz Centre for Environmental Research — UFZ, Theodor-Lieser-Straße 4, D-06120 Halle (Saale), Germany
Martin-Luther-University Halle-Wittenberg, Department 6 Research, 06099 Halle (Saale), GermanycDepartment of Environmental Sciences, University of California Riverside, Riverside, CA 92521-0424, USAba r t i c l e

i n f oArticle history:
Received 10 October 2007
Revised 20 August 2008
Accepted 3 October 2008
Available online 17 October 2008
Keywords:
Preferential ﬂow
Non-reactive tracer transport
Non-equilibrium
Modelling
Structured soila b s t r a c t
Although it has many positive effects, soil structure may adversely affect the ﬁltering function of
the vadose zone that protects natural water resources from various sources of pollution.
Physically based models have been developed to analyze the impacts of preferential water ﬂow(PF) and physical non-equilibrium (PNE) solute transport on soil and water resources. This
review compiles results published over the past decade on the application of such models for
simulating PF and PNE non-reactive tracer transport for scales ranging from the soil column tothe catchment area. Recent progress has been made in characterizing the hydraulically relevantsoil structures, dynamic ﬂow conditions, inverse parameter and uncertainty estimations,
independent model parameterizations, stochastic descriptions of soil heterogeneity, and 2D or
3D extensions of PNE models. Two-region models are most widely used across all scales; as astand-alone approach to be used up to the ﬁeld scale, or as a component of distributed, larger
scale models. Studies at all scales suggest that inverse identiﬁcation of parameters related to PF
is generally not possible based on a hydrograph alone. Information on ﬂux-averaged andspatially distributed local resident concentrations is jointly required for quantifying PNEtransport. At the column and soil proﬁle scale, model predictions of PF are becoming
increasingly realistic through the implementation of the 3D soil structure as derived fromhydrogeophysical and tracer techniques. At the ﬁeld scale, integrating effects of the soil
structure and its spatial variability has been attempted by combining 1D PNE approaches withstochastic parameter sampling. At the catchment area scale, the scarcity of data makes
validation of PF related model components a task yet to be accomplished. The quest for easily
measurable proxy variables, as ‘the missing link’ between soil structure and model parameters,
continues in order to improve the practical predictive capability of PF–PNE models. A follow-up
paper complementing this manuscript reviews model applications involving non-equilibriumtransport of pesticides, as representatives of reactive solutes.
© 2008 Elsevier B.V. All rights reserved.

Abbreviations: 1D, one-dimensional; 2D, two-dimensional; 3D, three-dimensional; BC, boundary conditions; BTC, break through curve; CD, convection–
dispersion; CDE, convection–dispersion equation; CLT, convective lognormal transfer function; CT, computer tomography; CV, coefﬁcient of variation; DEM, digital
elevation model; DFM, discrete fracture model; DPM, dual-permeability model; DP-MIM, dual permeability model with mobile and immobile zones; GIS,
geographic information system; IC, initial conditions; KDW, kinematic dispersive wave; KW, kinematic wave; MIM, mobile–immobile model; MRI, magnetic
resonance tomography and imaging; PC, personal computer; PF, preferential ﬂow; PNE, physical non-equilibrium; SPM, single-permeability model; STM, stream
tube model.
⁎ Corresponding author. Tel.: +49 345 5585 5406.
E-mail addresses: max.koehne@ufz.de (J.M. Köhne), jiri.simunek@ucr.edu (J. Šimůnek).
1
Tel.: +1 951 827 7854.
0169-7722/$ – see front matter © 2008 Elsevier B.V. All rights reserved.
doi:10.1016/j.jconhyd.2008.10.002

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

1. Introduction
The phenomenon of non-equilibrium preferential ﬂow (PF)and transport in fractures, macropores or other local highpermeability zones, is encountered in various environmental
systems. PF has received increasing attention in different ﬁelds,such as hydrogeology, petroleum engineering, (soil) hydrology,
and waste disposal sites construction. For instance, aquifer
systems in heavily fractured rock, or in limestone with
dissolution conduits, are often vulnerable to contaminants
and need special protection when used as sources for drinking
water (Birk et al., 2006). The remediation of heterogeneous
sedimentary aquifer systems contaminated by organic chemicals at former military or industrial sites is an intricate, longlasting, and costly task if the contaminants are distributed
across layers with contrasting hydraulic and chemical properties. Prolonged tailing of concentrations has been encountered
due to the slow release of residual contaminants from ﬁnetextured and organic (e.g., lignite) layers to coarse textured
layers (Grathwohl et al., 2004). For similar reasons, a considerable effort is required to extract residual amounts of oil fromsubsurface petroleum reservoirs in fractured rock or heterogeneous sediment systems (Al-Huthali and Datta-Gupta, 2004).
PF processes in the vadose zone of fractured rocks have also
been intensely studied (Dasgupta et al., 2006a), e.g., in theframework of assessing the suitability of potential nuclear
waste repositories (Wu and Pruess, 2000). In the agriculturalenvironment, PF and agrochemical transport in structured soils
may signiﬁcantly affect the growth of crops (Kramers et al.,2005) and water resources (Sinkevich et al., 2005).
A review on PNE processes, including model applications,
was presented a decade ago by Jarvis (1998). Overviews of therelevant processes and conditions triggering PF in soils are
discussed in the classic paper by Beven and Germann (1982)
and in the recent review by Jarvis (2007) which covers the stateof-the-art in our understanding of PF and (agro-) chemicaltransport, and the effects of soil and crop management
practices. The latter is also treated in another recent review(Ma et al., 2007). Finally, the value of PF to ecosystem services
was assessed with a global perspective by Clothier et al. (2008).
Researchers have attempted to translate process information gained from experimental observations into conceptual
theories and mathematical models that are ‘physically-based’or ‘mechanistic’. The governing equations are implemented incomputer simulation tools, which can then be tested, modiﬁed,
and/or (in)validated against new experimental data. Increasingly sophisticated models have been developed to analyze PF
in various environmental systems, where there exist considerable differences in spatial and temporal scales, type of ﬂowing
liquids (e.g., water, brine, oil) and transported agents (e.g.,
dissolved chemicals, NAPLs, colloids, particles, microbes). Given
the limited scope of this review, the focus is on vadose zone
models developed to assess the impact of human activities onwater and solute balances in the agrosphere.
The objective of this paper is to give an overview of recent(past decade) developments in the application of physicallybased models describing preferential water ﬂow and solutetransport through structured agricultural soils. The PF model
concepts are only brieﬂy explained here. For details, the reader
is referred to recent review articles by Šimůnek et al. (2003),
Gerke (2006), Jarvis (2007), and Šimůnek and van Genuchten

5(2008). Here, the emphasis is on the applications of such models
in describing experimental observations. This paper reviews
model applications to conditions involving PF and nonequilibrium transport of conservative solutes at different spatial
scales. First (Section 2), an overview of model concepts andcorresponding computer models is given. Then, model applications concerning preferential water ﬂow (Section 3) and solutetransport (Section 4) are reported, followed by (Section 5) a
discussion of two-region approaches and (Section 6) methodsof independent parameter determination. Subsections within
Sections 2, 3, 4, and 6 are organized with regard to increasing
spatial scale.
Modelling of ﬁngering ﬂow (e.g., De Rooij, 2000; Wang et al.,2004), funneling ﬂow (e.g., Ju and Kung, 1997), and reactive,
colloid or particle-facilitated solute transport are excluded fromthis contribution. Model applications involving transport ofpesticides, as representatives for reactive solutes, will bereported in a companion review article.
2. Overview of model approaches
2.1. Conceptual models
Classical physically based model concepts for water ﬂow
and solute transport in structured soils can be broadlyclassiﬁed into continuum, bi- or multi-continuum, and network approaches (Fig. 1).
The traditional conceptualization of soil is that of a porous
medium with continuous properties (‘Continuum’, Fig. 1a).
Variably saturated water ﬂow through such a porous system isuniform and at local equilibrium, and is often described using
the Richards equation. Bimodal or multimodal hydraulic functions (e.g., Othmer et al., 1991; Wilson et al., 1992; Ross andSmettem, 1993; Durner, 1994; Mohanty et al., 1997) are an
indirect way of representing soil structure in such a uniform ﬂow
approach. A one-domain PNE concept of inﬁltration (Fig.1b) was
realized by assuming a rate-limited equilibration between thesoil water content and the pressure head in the Richards
equation (Ross and Smettem, 2000).
Other conceptualizations of PNE regard the soil as being
composed of two or more regions (domains) (Fig. 1c–f). The
stream tube model (STM), or stochastic convective model,
assumes displacement of tracer at (usually) steady-state ﬂow
through independent ﬂow tubes (Fig. 1c); it can be used witha convective lognormal transfer function (CLT) (Jury and Roth,
1990).
Two-region models of soil conceptualize structural pores as
the fast ﬂow region, whereas the soil matrix has a lower or zerosaturated hydraulic conductivity. Multi-region approaches
consider additional regions — up to three so far in model
applications. Domains are not geometrically separated, but areconceptually superimposed over the same soil volume. Common to these approaches are zero-dimensional, usually ﬁrstorder, transfer terms for lateral exchange of water and/or solute
between the regions.
The mobile–immobile model (MIM) (Fig. 1d) assumes thatthe soil micropore network is poorly connected and of such alow permeability that water is immobile in the verticaldirection. The MIM approach was originally developed for
simulation of ﬂuid displacement in fractured petroleum
reservoirs (Warren and Root, 1963; Coats and Smith, 1964).

6

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Fig. 1. Conceptual models of water ﬂow and solute transport in structured soils. Model dimensions may vary, e.g. 1D, 2D or 3D (a.–f.), 2D or 3D (g.–j.), 3D (k.).
Sources: g. modiﬁed after Graf and Therrien (2006), h. modiﬁed after Vogel (2002), i. modiﬁed after Deurer et al. (2003), j. Anderson et al. (2000), reprinted withpermission from the Soil Science Society of America, k. Vogel et al. (2006), reprinted with permission of the authors.

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

The approach was later adapted for structured soils (vanGenuchten and Wierenga, 1976; Gaudet et al., 1977).
The dual-permeability model (DPM) (Fig. 1e) assumes thatthe porous medium consists of two overlapping poredomains, with water ﬂowing relatively fast in one domain(often called the macropore, fracture, or interporosity
domain) and slow in the other domain (often referred to as
the micropore, matrix, or intra-porosity domain). The DPM
approach was proposed initially for studying fractured rockreservoirs in hydrogeology and petroleum engineering(Barenblatt et al., 1960; Duguid and Lee, 1977).
The multi-region model (Fig. 1f) comprises three ﬂow
regions (e.g., matrix, primary and secondary fractures) with
two transfer terms (Gwo et al., 1995; Hutson and Wagenet,1995). The same concepts shown in Fig. 1a–e were applied, orwere combined in different pairs, to describe water and solute
movement (see e.g. Šimůnek et al., 2003).
The above approaches can be further classiﬁed regarding
their scope. The ﬁrst group comprises continuum and bicontinuum models which consider ﬂow in the entire soil
as being controlled by both capillary and gravity forces(‘capillary PF models’). Bi-continuum models assume that theRichards equation can be used also for the network ofstructural pores. Flow in the PF system can be directed
upwards during periods of evapotranspiration.
The second group includes two-domain models whichassume that ﬂow in the PF domain is controlled by gravity
only and is always directed downwards (‘gravity-driven PF
models’). PF is often described by the Kinematic Wave (KW)
equation. Models of this group are tailored to soils withmacropores or heavy clay soils with distinct shrinkage cracks.
A third group comprises simpliﬁed descriptions of macropore ﬂow, such as capacity (“tipping bucket”) or routing
approaches (‘empirical PF models’). The capacity approachassumes zero downward ﬂow from each layer until ﬁlled toﬁeld capacity. The routing approach assumes instantaneous
bypass ﬂow.
While the above concepts were often implemented inagricultural system models, other concepts were more rarely
applied, and for research only (Fig. 1c, g–k).
The ‘network’ conceptual models (Fig. 1g–k) here comprise three entirely different approaches. Discrete fracture
models (DFM, Fig. 1g) were developed for describing ﬂow andtransport in discrete fractures and in the matrix of fractured
rock (e.g., Wang and Narasimhan, 1985). Exchange of waterand solute between the fracture and matrix domains is
calculated as advection–diffusion process without a need toresort to the ﬁrst-order approximation. However, a dense
spacing of nodes in the numerical grid is required near
fractures, which triggered the development of advanced
numerical techniques for adaptive grid reﬁnement (Slough
et al., 1999; Li et al., 2005) and time discretization (Park et al.,2008). A fundamentally different, since microscopic, approach is the pore network model (Fig. 1h), which calculates
Poiseulle ﬂow through the void soil pores arranged in a 2D or
3D geometric network. Pore network models have demonstrated potential for prediction of hydraulic properties ofstructured soil (e.g., Vogel and Roth, 2001, 2003). Solutetransport calculations are based on either particle tracking or
local advection and diffusion in the pore network. However,calculations are still limited to steady-state ﬂow conditions at

7

sample dimensions of 10− 1 m. The drainage network (Fig. 1i)concept (Deurer et al., 2003) identiﬁes ﬂow path networks byspatially linking the scaling factors of water retention curves.
Fractal (Fig. 1j) and structure-based (Fig. 1k) concepts are
brieﬂy described later with their applications (Sections 3.2.1
and 4.2.3).
For applications at larger scales, distributed models (i.e.,models that account for spatial variations in variables andparameters) were developed (Fig. 2). Those models can befurther subdivided into (here so-called) hydrological unit
models (Fig. 2a), which are discretized into sub-catchments,
or smaller units, based on similar physical, geomorphologicaland hydraulics characteristics. In these units, the vadose zone
is represented by 1D or 2D continuum or two-domain
submodels. Other distributed approaches considered here
are the 3D surface/vadose zone/groundwater models (Fig. 2b)with 3D subsurface discretization and a PF description for thevadose zone. Ideally, all submodels are fully coupled byappropriate boundary conditions. Examples of computermodels (Table 1) based on the above concepts are describedin the following section, with focus on the PF components.
2.2. Computer models
2.2.1. Capillary PF models
The PF representations in the models below are closely
related to the concept and equations of the DPM of Gerke andvan Genuchten (1993). Thus, the following model descriptions can focus on deviations and modiﬁcations.
DUAL (Gerke and van Genuchten, 1993) is a 1D-DPM basedon two Richards and convection dispersion equations formatrix and fracture pore systems, coupled by ﬁrst-orderterms for bi-directional exchange of water and solute. Watertransfer is driven by the pressure head difference between theﬂow domains, while solute transfer proceeds by advection
and by diffusion driven by the inter-domain concentrationdifference. The original research model was adapted for
application to ﬁeld conditions (Gerke and Köhne, 2004).
DPOR-1D (Lewandowska et al., 2004) employs a single
integrodifferential mass balance equation that includes anintegral term for water exchange between matrix and PF
subdomains. Despite an increased accuracy, fewer BCs andparameters need to be deﬁned than in DUAL. In MURF and
MURT (Gwo et al., 1995, 1996), (vertical) MUlti-Region Flowand Transport models, three overlapping pore regions withdifferent soil hydraulic properties are considered. First-orderwater and solute exchange occurs between regions arranged‘in series’. TRANSMIT (Hutson and Wagenet, 1995) is anothersimilar three-region 1D model. HYDRUS-1D (Šimůnek et al.,1998, 2003, 2005, 2008a,b,c; Šimůnek and van Genuchten,2008) describes water, heat, and solute movement in thevadose zone. The latest version includes a variety of SPM,
MIM, and DPM approaches; the more frequently used MIM
and DPM are again based on DUAL (Gerke and van Genuchten,1993). Inverse parameter optimization is implemented using
the Marquardt–Levenberg method which can accommodate
different types of data in the objective function.
DPOR-2D (Szymkiewicz et al., 2008) is the 2D equivalentof DPOR-1D. 2DCHEM (Nieber and Misra, 1995) is a 2D DPM
for application in drained soils. S_2D_DUAL (Vogel et al.,2000) may additionally consider spatially distributed

8

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Fig. 2. Distributed models. a) Hydrological unit approach, b) 3D surface/vadose zone/groundwater models.

parameters in the 2D grid. HYDRUS-2D (Šimůnek et al., 1999),
and its latest update HYDRUS (2D/3D) (Šimůnek et al., 2006),
are the two- and two/three-dimensional equivalents of
HYDRUS-1D, respectively. The origin of HYDRUS models canbe traced back to the early 1970s (Neuman, 1972; see Šimůnek
et al., 2008b). Additional features in HYDRUS-2D concern
boundary conditions (e.g., tile drainage) and tools forcalculating spatially distributed model parameters. The
TOUGH model family is mentioned for completeness. The
integral ﬁnite differences method allows TOUGH2 to be used
as a DPM with multi-phase (e.g., water/air) ﬂow equations in
one to three dimensions. TOUGH codes have been extensively
used for problems of ﬂow and chemical transport in fractured
rock (Pruess, 2004).
2.2.2. Gravity-driven PF models
A kinematic wave (KW) one-region model (Beven andGermann, 1982; Germann, 1985) is based on the boundary
layer ﬂow theory and was used for describing PF (e.g.,
Germann and Di Pietro, 1996, 1999). The KW model assumes
that the wetting front proceeds by convective ﬁlm ﬂow in themobile region and does not exchange water with theimmobile region. The kinematic dispersive wave (KDW)
model (Di Pietro et al., 2003) was developed by applying asecond-order correction to the KW model. MACRO (Jarvis,
1994; Larsbo and Jarvis, 2003; Larsbo et al., 2005) is a 1D DPM
that combines a KW description of water ﬂow and soluteconvection for the macropore region with Richards waterﬂow and solute convection dispersion in the matrix. Watertransfer into the matrix is treated as a ﬁrst-order approximation to the water diffusion equation and is proportional to thedifference between actual and saturated matrix water contents. The latter is deﬁned at a particular boundary pressure
head that partitions the total porosity into matrix pores andmacropores. Reverse transfer of water in excess of matrix
saturation is instantaneously routed from the matrix into
macropores. In an analogue fashion, solute advection anddiffusion are calculated as ﬁrst-order rate processes. MACRO
considers swell–shrink dynamics of clay soils. The KW modeland MACRO theoretically calculate identical hydrographs for
saturated matrix conditions without lateral water transfer.
The version MACRO 5.0 (Larsbo and Jarvis, 2003; Larsbo et al.,2005) is Windows-based and provides options for MonteCarlo uncertainty estimation and for parameter estimation
based on the SUFI program (Abbaspour et al., 1997).
SIMULAT (Diekkrüger, 1996) is a 1D DPM that considers
gravitational ﬂow and solute convection in the macropores,
advective solute exchange with the matrix, and vertical
Richards ﬂow and solute transport by convection and dispersion in the matrix. The Root Zone Water Quality Model, RZWQM(Ahuja et al., 2000) utilizes a DP-MIM description of 1D verticalsoil water ﬂow and chemical movement. Three transport
regions are assumed to exist in the soil: cylindrical macropores,
the mobile soil matrix, and the immobile soil matrix. In themacropores, water ﬂow is calculated using the Poiseulle
equation and solutes are displaced by convection. In the mobile
matrix region, water ﬂow is described using the Green-Ampt
approach during inﬁltration and the Richards equation during
redistribution, while solute moves by convection. Watertransfer is assumed to proceed uni-laterally from the macropores into the matrix by (Green-Ampt) imbibition, which can be
restricted by a sorptivity (impedance) factor. Solute exchange

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

9

Table 1
Characteristics of common 1D models for preferential ﬂow and tracer transport a
MACRO

HYDRUS-1D b

CRACK-NP

SIMULAT

Documentation (Ahuja et al., 2000;
Malone et al., 2004)

Larsbo and Jarvis, (2003, 2005)

Armstrong et al.
(2000)

Model type(s)Water ﬂow,
matrixSPM, DPM
Richards(Šimůnek et al., 2005; 2008a,Šimůnek and van
Genuchten, 2008)
SPM, MIM, DPM, DP-MIM
Richards

Aden and
Diekkrüger(2000)
DPM
Richards

RZWQM

Water ﬂow,
PF paths
Potential ET

Real ET
Solutetransport,
matrixSolutetransport,
PF paths
Water transfer

DPM, DP-MIM
Green-Ampt(downward), Richards(redistribution)
Poiseulle (gravity ﬂow
in cylindricalmacropores)
Penman Monteith orprescribed input

MIM
Inﬁltration: Philips(but immobile waterin peds within soil)
Poiseuille (gravity
Kinematic wave
ﬂow in cracks)(gravity ﬂow inmacropores)
Prescribed input
Penmandata
Monteith orprescribed
input
No
Reductionfunctions for ET
Immobile
Convection–
dispersion

Kinematic wave
(gravity ﬂow in macropores)

Richards

Penman Monteith orprescribed input

Penman Monteith, Hargreaves,or prescribed input

Reduction functions
for ET
Convection

Reduction functions for ET

Reduction functions for ET

Convection–dispersion

Convection–dispersion

Convection

Convection

Convection–dispersion

Convection

Convection

First-orderFirst-order

Philip inﬁltration

Darcy

First-order (convection–diffusion)

First-order (convection–diffusion)

Diffusion (Ficks law)

Convection

Yes
Yes
Water, solute(Yes) — lumped dissipation rate
(Yes)
No
Tillage

Yes
Yes
Water, solute
Yes
(Yes)
No
No

No
Yes
Water
No
(Yes)
No
No

Yes
Yes
No
No
No
No
No

Yes
SUFI

Yes
Levenberg–Marquardt

(Yes)
No

Yes
No

Yes

Yes

No

No

Green-Ampt, one-way(macro-matrix)
Solute transfer matrix-macro: instantmixing; mo-im:
ﬁrst-order (diffusion)Heat ﬂux
Yes
Root growth
Yes
Root uptake
Water, solute
Volatilization
Yes
Surface runoff
Yes
Solute in runoff Yes
Management
Crop rotations, tillage,
options
fertilizer, manure,pestic. application
mode
Tile drainage
Yes
InverseNo
methodGraphical User Yes
Interfacecrop yield, nutrient
Selectedcycle, equilibrium
furtherchemistryfeatures

Particle transport, CO2 transport,
Swelling–shrinking, consolidation,
Swelling–shrinking,surface sealing, canopy intersection
canopy interception
chemistry of carbonate system,
and washoff, Monte Carlo uncertainty major ion chemistry, metabolite fateestimation, Metabolite transport(ET — evapotranspiration, T — temperature, mo — mobile region, im — immobile region, macro — macropores, PF — preferential ﬂow, θ — soil moisture content,(Yes) means: is considered in a simpliﬁed way, SPM — single porosity model, MIM — mobile–immobile water model, DMP — dual-permeability model, DP-MIM —dual-permeability model with immobile water).
aModel features relevant to transport of reactive solutes such as pesticides and management are shown in the companion paper.
b
Similar features are included in HYDRUS-2D and HYDRUS (2D/3D).

proceeds by instant mixing between the macropore solution
and a thin matrix layer (e.g., 0.6 cm), and by ﬁrst-order diffusion
between mobile and immobile matrix regions (Malone et al.,2001).
VIMAC. (Greco, 2002) is a three ﬂow domain model which
considers ﬂow in swelling and shrinking clay soils. Flow in thematrix is modelled by means of Darcy equation. Water ﬁlm
ﬂow in shrinkage cracks of variable aperture depending on
the matrix saturation is described using the KW equation, as is
ﬂow in permanent macropores, a fraction of which may
be deﬁned as dead-end pores with internal catchment.
Horizontal water exchange between the three regions is

described by ﬁrst-order equations taking into account partial
wetting of macropore and cracks, and other effects (Greco,
2002).
CRACK-NP. (Armstrong et al., 2000) is a MIM that describes 1D
vertical water ﬂow and solute transport in clay soils withplanar shrinkage cracks and peds (soil aggregates) of cubic
shape. The crack and matrix porosities are dynamic functionsof the soil moisture. Water ﬂow in the cracks is calculatedusing Poiseuille's law modiﬁed to account for the cracks'
tortuosity and connectivity (Germann and Beven, 1981).
Lateral water uptake of peds is calculated using the Philips
equation modiﬁed for the fractional wetted aggregate surface

10

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

area. No vertical ﬂow is considered in peds. Solute transport inthe cracks proceeds by convection, with lateral advective anddiffusive losses of solute to peds, and potential reverse
diffusive transfer from peds into cracks.
M-2D was developed by Kohler et al. (2001) by combining
the 2D Richards-CDE based SWMS-2D model (Šimůnek et al.,1994) with a 1D DPM adopted from MACRO (Jarvis, 1994).
First-order water transfer is described similarly as in Gerkeand van Genuchten (1993), whereas solute transfer is byadvection only.
2.2.3. Empirical PF models
Physically-based two-region approaches with parsimonious parameter demand were developed. Two types can be
distinguished with respect to treatment of PF: instantaneousrouting from the surface to the bottom (PF routing approach),
and a layer-wise cascade of macropore water in excess of acertain saturation threshold (capacity approach).
2.2.3.1. PF routing approaches.
The FRACTURE submodel ofthe HYDRUS-1D (Novák et al., 2000) formally divides inﬁltration into the soil matrix into two components: verticalinﬁltration through the soil surface (1D-Richards equation),and lateral inﬁltration via soil cracks (Green-Ampt approach).
Excess water that does not inﬁltrate the soil surface is routed
into soil cracks. The Soil Water Atmosphere and Plant (SWAP)
model (van Dam et al., 1997; van Dam, 2000) comprises 1D
Richards and CD equations for the matrix and a routine forrouting water and solute in dynamic shrinkage cracks, coupledby ﬁrst-order water and diffusive solute transfer. Lateralinﬁltration from the cracks into the soil matrix is calculated asa ﬁrst-order process and depends on the active crack area. An
express fraction of solute may be routed directly to the drains.
2.2.3.2. Capacity approaches.
The COUP model (Jansson andKarlberg, 2004) describes matrix ﬂow using the Darcy and thecontinuum equations, while inﬁltration excess of the soil
matrix is cascaded as bypass ﬂow to lower soil compartments.
Water inﬁltration into frozen soil and heat transport are also
considered. COUP combines the SOIL (Jansson, 1991) and
SOILN (Eckersten et al., 1998) models. The LEACHW model
based on the Richards equation (Hutson and Wagenet, 1992)
was extended with a description of water moving downward
through fractures in a tipping bucket fashion (Booltink, 1994).
The Solute Leaching Intermediate Model, SLIM, is a MIM of
capacity type. SLIM has two main parameters: the amount ofimmobile water and the permeability (Addiscott et al., 1986;
Addiscott and Whitmore, 1991). The 1D-DPM IN3M (Weiler,
2005) calculates inﬁltration into the matrix and horizontal
water transfer using Green-Ampt equations, and assumes theDarcian ﬂow for redistribution. Inﬁltration exceeding the soil
matrix capacity triggers surface initiation of PF. The excesswater is distributed into the macropores according to thestochastic Weibull frequency distributions of observed
macropore drainage areas. Additionally, subsurface PF initiation may occur at layers boundaries. Macropore ﬂow is routed
downwards by layer-wise mass balance (Weiler, 2005). Basedon concepts of Bouma (1990), Kohler et al. (2003) developed
two simple capacity DPM of water ﬂow: the surface trigger
model assumes that gravitational PF is generated at the soil
surface, while the subsoil trigger model assumes that PF iscreated at the topsoil-subsoil boundary when the water
storage capacity of the topsoil is exceeded.
2.2.4. Catchment scale models
2.2.4.1. 3D surface/vadose zone/groundwater models. MODHMS
is a 3D water ﬂow model based on the 3D Richards equation(Panday and Huyakorn, 2004). Areal overland ﬂow and ﬂow
through a network of streams and channels is described bydiffusion wave equations. A bypass routing approach accounts
for PF to tile drains that depends on weekly rainfall, irrigation,
and an empirical bypass fraction parameter. FRAC3DVS is a 3D
DFM (3D matrix, 2D parallel fractures) for variably saturated(‘VS’) ﬂow and solute transport in fractured porous formations(Therrien and Sudicky, 1996). Vertical, horizontal and inclined
2D fractures are incorporated into the 3D grid by superimposing
2D parallel plate elements onto 3D matrix block elements.
Fracture faces and matrix blocks share common nodes along thefracture walls, such that continuity of hydraulic head andconcentration is ensured in the calculation of mass exchanges.
Water ﬂow is calculated using Richards equation in the matrix
and an adapted form of Richards equation for the fractures.
Transport is by advection and dispersion in both the fractures
and the matrix. FRAC3DVS can be regarded as an extension of
the 2D FRACTRAN groundwater model (Sudicky and McLaren,1992). A recent successor of FRAC3DVS is the 3D fully integrated
subsurface and overland ﬂow and solute transport model
HydroGeoSphere (Therrien et al., 2005). The Integrated
Hydrology Model, InHM (VanderKwaak, 1999; VanderKwaak
and Loague, 2001), also gives a fully coupled description of thesurface and subsurface water ﬂow and solute transport over the
2D land surface and in the 3D dual continuum subsurface undervariably saturated conditions. Overland ﬂow is calculated using
the 2D nonlinear diffusion-wave equation together with
Manning's equation. In the subsurface 3D Richards equation isassumed; optionally in continuum or dual-continuum form.
2.2.4.2. Hydrological unit models.
WEC-C (Croton and Barry,
2001) is a distributed water and solute displacement model ina GIS-like environment. Vertical PF can be alternatively
considered as a 1D DPM (Gerke and van Genuchten, 1993)
or 3D DPM. Groundwater ﬂow is modelled as a 2D Darcy
system and overland ﬂow using 2D Manning and kinematicwave equations. Boundary conditions account for interception, canopy storage, transpiration, evaporation, stream
routing and groundwater discharge. WEC-C was validated
against synthetic data examples (Croton and Barry, 2001).
CATFLOW (Zehe et al., 2001), based on the original model
HILLFLOW (Bronstert, 1994), subdivides the catchment area
into a number of hillslopes connected to a drainage network.
Each hillslope is conceptualized as 2D vertical cross-sections
through its main slope. Subsurface ﬂow and transport in eachhillslope is described by 2D Richards-CDE equations. The
effect of macropores on water ﬂow is represented by a linear
increase of the bimodal hydraulic conductivity function above
a saturation threshold of 0.8. Surface runoff along hillslope
elements, subsurface drainage, and runoff in the channel
network is calculated based on the advection diffusion
approximation of the 1D Saint Venant equations. The
Distributed Hydrology Soil Vegetation Model, DHSVM (Wigmosta et al., 1994), was extended with PF descriptions

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

(Beckers and Alila, 2004). DHSVM uses the 1D Green-Amptequation to calculate vertical unsaturated matrix water
movement. Macropore ﬂow is initiated by surface ponding
and is instantaneously routed to the bedrock surface. Lateralgroundwater ﬂow in the matrix is calculated using a quasi 3D
Darcian approach, while lateral PF in the groundwater iscomputed using a Darcy-based routing approach. Preferential
and matrix ﬂows are disconnected. Overland ﬂow is calculated using a linear reservoir (capacity) scheme. The 1D
Distributed Subsurface-Surface Flow Dynamics Model,
DSFDM (Mulungu et al., 2005) adds a PF component to the
Famiglietti and Wood (1994) catchment scale model. Threesoil layers are considered in DSFDM: in the top 1-m soil layer,
macropore ﬂow and lateral interﬂow may occur; the middle‘transmission’ zone assumes vertical ﬂow in the soil matrix,
and the bottom layer with Darcy (lateral) ﬂow of groundwaterin the soil matrix. Preferential inﬁltration through macropores and interﬂow in pipes, as well as overland ﬂow, aredescribed using KW equations.
3. Modelling of preferential ﬂow systems
A primary reason for the interest in PF is its impact on
contaminant transport. However, a number of studies have
utilized water ﬂow observations only, without solute observa-

11

tions, to evaluate PF models and to study PF processes atdifferent scales. The review of such model applications (Table 2)should indicate the model complexity and data required todescribe PF under the given conditions.
3.1. Scale I (soil column, lysimeter)
3.1.1. Constructed soil columns
The advantage of using soil columns (length approx. 0.1 m
to 1 m) for a basic ﬂow process analysis is that IC and BC canbe optimally controlled, and the heterogeneity of soil properties is negligible. Water ﬂow experiments were conductedwith dual-porosity media with known geometry to eliminate
unknowns that might otherwise obscure the PF process understudy. The water transfer term and domain speciﬁc hydraulicfunctions in capillary PF models (DPM) were evaluated.
Continuum models were also tested regarding their capabilityof simulating PF. Examples are discussed below.
3.1.1.1. Water transfer term.
Ahuja et al. (1995) used RZWQM
to describe the effect of macropore ﬂow under rainfall in a soilcolumn with a 3-mm artiﬁcial macropore. In order to matchthe macropore bottom outﬂow, the lateral absorption of waterfrom macropore to soil matrix had to be adjusted forcompaction of the macropore wall (Ahuja et al., 1995).

Table 2Model applications for preferential WATER ﬂow
Scale Model type; nameDimension, equation name or description

Authors of model application

Scale IAC
SPM, MIM, DPM; HYDRUS-1D

1D; Richards (2× for DPM), f.o. (or s.o.) water transfer

AC
AC a
UC
UC
UC
L

1D; single integro-differential equation
1D; Convective dispersive kinematic wave, no water transfer
1D; Richards, kinetic θ–h relation
1D; Richards (matrix), KW (water ﬁlm in cracks), KW (macropores)
1D; Richards (2× for DPM), f.o. water transfer
1D; Richards (matrix), KW (macropore), f.o. water transfer(Castiglione et al., 2003; Köhne and Mohanty,
2005; Köhne et al., 2006a)
Lewandowska et al. (2005)
Di Pietro et al. (2003)
Ross and Smettem (2000)
Greco (2002)Šimůnek et al. (2001)(Jabro et al., 1998 (L); Ludwig et al., 1999 (F);Alaoui et al., 2003 (UC); Akhand et al. 2006 (P))

DPM; DPOR_1D
MIM; KDW model
SPM; n.a.
TPM; VIMAC
SPM, MIM, DPM; HYDRUS-1D
SPM, DPM; MACRO

Scale IIS
DPM/S_1D_Dual
S
DPM (stochastic)/IN3M
P
SPM/MIN3P + Rosetta
P
MIM/Fractal active region
modelP

1D; Richards (2×)
1D; Green-Ampt, Darcy
1D; Richards, modiﬁed K(h)1D; Richards, θ-h is a function of active (i.e., mobile) region,
no water transfer
0D-distributed, initiation of macropore ﬂow at surface or subsurface
according to geostatistic microtopography and macropore distribution

Zumr et al. (2006)
Weiler (2005)
Gérard et al. (2004)
Liu et al. (2005)
Weiler and Naef (2003)

Scale IIIH
DPM/Hillslope PF modelW, H DPM/DHSVM
W, S SPM (distributed)/CATFLOW

Tsutsumi et al. (2005)
Beckers and Alila (2004)(Zehe and Blöschl, 2004; Zehe et al. 2005)

W
W
R

Vrugt et al. (2004)
Bormann et al. (2005)
Mulungu et al. (2005)

3D-Richards (matrix), 1D-Manning (pipe ﬂow)
(Pseudo) 3D; Green Ampt, routing (no transfer), Darcy
2D-vertical Richards (DEM and GIS based) with bimodal K(h) for
hillslope elements, and 1D simpliﬁed Saint Venant equations
for surface runoff and drainage
DPM/MODHMS + PF
3D-Richards (matrix), 1D-bypass routing (macropore)
SPM, DPM/TOPLATS&SIMULAT 3D (quasi); 2D-(TOPLATS)& 1D (SIMULAT)DPM/DSFDM
1D-distributed, KW

Scale I: artiﬁcial or repacked soil column (AC), undisturbed column (UC), lysimeter (L); scale II: soil proﬁle (S), plot (P), hillslope (H), ﬁeld (F); scale III: watershed(W), river basin (R); PF — preferential ﬂow, DEM — digital elevation model, SPM — single porosity model, DPM — dual permeability model, TPM — triple
permeability model, MIM — mobile immobile model, KW — kinematic wave, KDW — kinematic dispersive wave, f — fast ﬂow region, m — matrix, f.o. — ﬁrst-order,
s.o. — second-order, SCM — stochastic convective model, STM — stream tube model, DFMD — discrete fracture matrix diffusion model, SOTS — second-order twosite model, MRM — multireaction model, BIM — Lagrangian transport model with bimodal lognormal probability density functions (pdfs), TFM — lognormaltransfer function model.
aRepacked soil columns inoculated with earthworms.

12

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Köhne and Mohanty (2005) tested the effect of having aﬁrst- or second-order water transfer term (Köhne et al.,2004a) in the DPM in HYDRUS-1D (Šimůnek et al., 2003).
DPM simulation results were compared with those obtainedusing the 2D Richards equation based HYDRUS-2D (Šimůnek
et al., 1999) as a reference ﬂow model. A column setup (80 cm
height, 24 cm diam.) provided outﬂows separately from thematrix and the central cylindrical PF region, as well aspressure heads and water contents in the PF and matrix
regions at various positions to enable PF model inversion. The
DPM could describe the observed hydraulic nonequilibrium
between the matrix and the PF domain. Region-speciﬁc
results came closer to those of the 2D-reference model
when using the second-order mass transfer term (Köhne andMohanty, 2005).
In a comparable experimental study, the same DPM withﬁrst-order term could model water transfer accurately only
when used with a scaling factor that varied with themacropore hydraulic conductivity (Castiglione et al., 2003).
3.1.1.2. Domain speciﬁc hydraulic functions.
Lewandowska
et al. (2005) conducted inﬁltration experiments with adouble-porosity medium composed of a mixture of sand andsintered clayey spheres arranged in a periodic manner. The
resulting inﬂow and outﬂow observations were used to verify
the DPM approach ‘DPOR-1D’ (Lewandowska et al., 2004).
With domain speciﬁc hydraulic parameters (van Genuchten,1980) identiﬁed independently for the sand material andclayey spheres using HYDRUS-1D (Šimůnek et al., 1998), themacroscopic prediction satisfactorily approximated observed
inﬂow and outﬂow for the composite dual-permeability
medium (Lewandowska et al., 2005). A similar procedure forparameterizing a DPM resulted in a satisfactory simulation ofwater drainage from a repacked soil column with cylindrical
PF region (Köhne et al., 2006a). The 2D counterpart of DPOR1D, named DPOR-2D, was successfully tested using dataobtained from an axi-symmetrical experimental inﬁltration
setup with a similar mixture of sand and sintered clay spheres(Szymkiewicz et al., 2008).
3.1.1.3. Continuum models.
Different (single-) continuum
models were also evaluated for their performance in describing preferential water drainage. For example, the KDW model(Di Pietro et al., 2003) was tested using water ﬂow observations obtained in laboratory experiments conducted with
repacked soil columns (28 cm height, 15 cm diam.) inoculatedwith earthworms. With model parameters estimated fromthe observed drainage hydrograph, the KDW model described
the outﬂow hydrographs accurately (Di Pietro et al., 2003). In
another example, the Richards model with bimodal functions(Durner, 1994) performed similarly well as the DPM in
describing water outﬂow from the bulk soil (Köhne andMohanty, 2005).
The success of the continuum models prompted thequestion of whether water ﬂow measurements alone are
sufﬁcient to obtain reliable values of hydraulic DPM parameters by inverse identiﬁcation. Using the experimental
setup of Köhne and Mohanty (2005), Köhne et al. (2006a)
found that the outﬂows from either the matrix or the PF
region could only be described when the DPM was ﬁtted toregion-speciﬁc outﬂow data. It was concluded that for naturalsoils where experimental data do not come in ‘separated
form’ for fast and slow ﬂow regions, domain-related hydraulicparameters of the DPM may be difﬁcult to identify from waterﬂow observations alone (Köhne et al., 2006a).
3.1.2. Intact soil columns and lysimeters
Undisturbed, naturally structured soil cores and lysimeters
more closely represent natural ﬁeld conditions than repacked
soil columns with artiﬁcial macropores. Both capillary andgravity-driven PF models (see Section 2.1) were evaluatedbased on experimental data from such conditions. The
following examples compare equilibrium and physical nonequilibrium models. Different types of data were used in theanalyses, including outﬂow, water contents, and pressure
heads.
3.1.2.1. Capillary PF models.
Ross and Smettem (2000)
applied their one-region nonequilibrium model (Fig. 1b) and acorresponding uniform ﬂow model to describe TDR measured
water contents and outﬂow responses from large (0.67 m
height, 0.3 m diam.) undisturbed soil columns subject toconstant rainfall. The surface 0.2 m was a ﬁne sandy loam andthe subsurface a kaolinitic clay. The presence of PF wasinferred from commencement of outﬂow prior to complete
wetting of the soil. The uniform ﬂow model failed to describe
this behaviour. By contrast, the PF model could describe both
cumulative outﬂow and water contents for ﬁve out of six
other soil cores, when only the macropore hydraulic conductivity was individually calibrated (Ross and Smettem,
2000).
Šimůnek et al. (2001) implemented the Marquardt–
Levenberg parameter optimization in HYDRUS-1D for inverse
DPM, MIM and SPM simulations of upward inﬁltration
experiments performed on undisturbed sandy loam soilcores (0.10 m height and diam.), without macropores,
sampled from a tilled soil horizon. The multi-objectivefunction was deﬁned in terms of the cumulative upwardinﬁltration rate, pressure head readings, and the ﬁnal total
water volume in the soil core. The observations could be
similarly well characterized using MIM and DPM, as opposed
to using the Richards equation based SPM, but the data wereinsufﬁcient to achieve uniqueness for many (up to 9)
estimated parameters required for two-domain models(Šimůnek et al., 2001).
3.1.2.2. Gravity-driven PF models.
Alaoui et al. (2003)
compared MACRO (Jarvis, 1994) and a corresponding onedomain KW approach to assess the extent of PF in a structured
soil column (0.45 m height, 0.39 m diam.). The soil column
contained aggregates, roots, and earthworm channels andtexture changed from loamy sand (top) to sand (bottom).
Measured water contents using TDR and outﬂow were studied
during inﬁltration into unsaturated soil. Outﬂow did not startearlier than the TDR response. The model results suggested
intense exchange between macropores and the permeable
sandy matrix. For these hydraulic equilibrium conditions, with
only one calibrated parameter, KW results came fairly close to
those using MACRO (with six calibrated parameters). However, the KW simulation showed step-increases of watercontents and outﬂow that were not apparent in the data
(Alaoui et al., 2003). This convective ﬂow effect was removed

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

by second-order correction to the KW model (Di Pietro et al.,2003).
VIMAC was applied to model inﬁltration experiments
carried out on a large undisturbed clay soil column (55 cm
high, 30 cm diam.) and lasting from 6 h to 8 days. Hydraulicmodel parameters were derived from independent measurements of hydraulic conductivity and water retention curves,
shrinkage characteristics, aggregate water diffusivity, and soil
morphology. The permanent macropores' cross-section areas
were calibrated. VIMAC was capable of reproducing outﬂow
and TDR-measured water content proﬁles. Results improved
when the internal catchment process in dead-end macropores was included (Greco, 2002).
Jabro et al. (1998) compared the gravity-driven PF model
MACRO (Jarvis, 1994) with the empirical two-domain PF
models SLIM (Addiscott et al., 1986; Addiscott and Whitmore,1991) and SOIL (Jansson, 1991), and with continuum models
LEACHW (Hutson and Wagenet, 1992) and NCSWAP (Molina
and Richards, 1984), to evaluate the signiﬁcance of considering PF for simulating water drainage. The silt loam soil
contained worm holes, decayed root channels and cracks. The
models were compared with regard to their accuracy inpredicting drainage ﬂuxes below the root zone of corn asmeasured with 18 pan lysimeters (area 0.465 m2, depth 1.2 m)
during 5 years. The ﬁrst of these years was taken for model
calibration. Simulated drainage ﬂuxes for all ﬁve models were
not signiﬁcantly different, and were all within 95% conﬁdence
intervals of the measured values (Jabro et al., 1998).
3.1.3. Summary for preferential water ﬂow modelling at scale I(column, lysimeter)
Most PF models applied at scale I are one-dimensional.
The typical time scale of the experiments is between hours
and days. Both continuum models and non-equilibrium
approaches with one or two domains were used. For artiﬁcial
dual-porosity media, PF processes as represented by region
speciﬁc observations of ﬂow, water contents, and pressure
heads could be successfully described with various DPM. A
second-order water transfer term could further enhance theaccuracy for certain geometries of the ﬂow system.
If the simulation objective was to describe the hydrographonly, simpler one-domain equilibrium ﬂow models often
performed equally well. Thus, structured soils may often be
regarded as a continuum for the purpose of modelling water
ﬂow. Consequently, inverse identiﬁcation of two-domain
model parameters generally yielded poor or uncertain results
when using data of water outﬂow or inﬁltration. Hence, ahydrograph must be considered insufﬁcient to make meaningful statements about PF. Conversely, combined data on
ﬂow and soil hydraulic state variables include PF informationthat can only be described by a non-equilibrium model.
3.2. Scale II (plot, ﬁeld)
Compared to scale I, additional uncertainties for modellingat scale II (plot, ﬁeld) arise due to the spatial heterogeneity
and temporal dynamics of the system with regard to soils,
crops, meteorological variables, boundaries, and the lateralsurface or subsurface runoff.
The plot and ﬁeld scales are not exactly deﬁned terms.
From an agronomic perspective, the ‘ﬁeld scale’ is the smallest

13

management unit that would represent the entire agriculturalﬁeld and its crop. In line with more common terminology insoil hydrology, the ﬁeld scale is deﬁned here as covering aconsiderable part of a ﬁeld (0.1 to 100 ha) under the prevailingboundary conditions, while the plot scale refers to a small partof a ﬁeld (1 m2 [soil proﬁle] to 1000 m2) often subject toexperimentally controlled boundary conditions (irrigation).
Depth is in the meter range. Several studies involved tiledrained ﬁelds, because of the straightforward delineation of
domain boundaries, the integral function of the drain outﬂow,
and the potential environmental effect of drains on the diffuse
pollution by agrochemicals.
The main conceptualizations of soil (texture and structure)
heterogeneity for water ﬂow studies involved (i) deterministic
application of a 1D PF-model with a single set of parameters,and (ii) stochastic application of a 1D PF-model with several setsof parameters sampled from their respective frequency
distributions. Approach (i) was frequently adopted due to its
relative simplicity. However, calibrated ‘effective’ parametervalues are lumped with heterogeneities and other scale effects,
and thus may lose their physical meaning. Approach (ii) allows
for uncertainty analysis of simulation results. Uncertainty may
have various sources, including uncertainty of parameters,
model structure, initial and boundary conditions, and user
subjectivity. An alternative to approach (ii) is the distributedmodel.
3.2.1. Plot scale
At the plot scale, typically the 1D deterministic approach(i) was applied. The same problems of identifying PF
parameters were encountered as on the column scale, so
only brief examples will be given. Ludwig et al. (1999)
modiﬁed and calibrated MACRO (Jarvis, 1994) for simulationof groundwater levels and matric potentials observed during
2 weeks in tile-drained aggregated clay and sandy loam soilplots (15 m by 100 m). DPM simulation results improved only
slightly compared to those using a one-domain (Richards)
model, and the water transfer rate parameters in the DPM
could not be identiﬁed (Ludwig et al., 1999).
Zumr et al. (2006) applied the 1D counterpart of the DPM
S_2D_Dual (Vogel et al., 2000) for the Levenberg–Marquardt
estimation of soil hydraulic parameters using measured
pressure heads below a plot with pseudogley soil containing
desiccated cracks, biopores and interaggregate fractures.
Problems with reproducing observations were attributed tothe uncertainty of allocating tensiometer values to the matrix
or macropores (Zumr et al., 2006).
Gérard et al. (2004) tried to predict 4 years of daily
measurements of soil water content in a forest soil plot withPF in a capillary pore network. They utilized pedotransfer
functions (ROSETTA, Schaap et al., 2001) and a bimodal
hydraulic conductivity function (Mohanty et al., 1997) with
the Richards' ﬂow submodel of MIN3P (Mayer et al., 2002).
However, soil water contents were under-predicted in the wet
range. The match improved after calibration, but calibrated PF
parameters appeared to be of limited physical meaning(Gérard et al., 2004).
The seasonal dynamics of structured soil hydraulic properties (Messing and Jarvis, 1990, 1993; Dasgupta et al., 2006b)
pose an additional challenge. Akhand et al. (2006) applied
MACRO to simulate 60 days of hourly readings of water

14

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

contents, matric potentials, and tile drain discharge from asilt–loam soil plot (7 m by 65 m) with cracks and worm
burrows. MACRO was calibrated for relatively wet autumn
conditions. In the drier validation period of spring andsummer when soil cracks developed, the model underestimated bypass ﬂow and tile discharge (Akhand et al.,2006).
There is increasing evidence that PF patterns in soil can becharacterized by fractal or multi-fractal geometries (Rieu andSposito, 1991; Hatano and Booltink, 1992; Baveye et al., 1998;
Boufadel et al., 2000; Olsson et al., 2002, 2007). Liu et al.
(2005) incorporated fractal PF patterns into a MIM. Fractal
information was included in terms of a constitutive relationship for the fraction of mobile water as a function of thesaturation. A preliminary validation of the model was
performed based on ﬁeld observations (van Dam et al.,1990) of the variation of mobile water content with soil
depth (Liu et al., 2005).

macropore density primarily controlled the total macropore
ﬂow, while surface micro-topography controlled the probability distribution; only a few macropores contributed
signiﬁcantly to the total macropore ﬂow. Qualitative agreement was found between results from model simulations andfrom several published laboratory experiments (Weiler andNaef, 2003).
Tiemeyer et al. (2007) developed and tested the spatiallydistributed model MHYDAS-DRAIN to assess the effects of tile
drainage and PF on the hydrology of lowland catchments. The
PF component was modelled using the transfer function
approach. The model was applied to two ﬁelds (4.2 and4.7 ha) in the Pleistocene lowland of North-Eastern Germany.
Modelled discharge rates and groundwater levels agreed
reasonably well with measured data, while parameter values
were found to depend on the calibration criteria as well as on
the spatial and temporal resolution of the modelling domain(Tiemeyer et al., 2007).

3.2.2. Field scale
Examples of the 1D stochastic approach (ii) and hydrological unit approach are presented as ﬁeld scale applications.
PF models traditionally do not consider a population ofmacropores, but assume a single set of properties representing some ‘equivalent’ PF path at the ﬁeld scale. Booltink(1994) was among the ﬁrst to represent heterogeneity in
bypass ﬂow modelling stochastically. The LEACHW model(Hutson and Wagenet, 1992) was extended with a tippingbucket description of fracture ﬂow in a clay soil at a tile
drained site of 100 m × 300 m. The Monte Carlo analysis was
performed for soil hydraulic parameters related to PF. Theirfrequency distributions were partly derived from measurements including the contact areas of soil aggregates adjacentto PF paths, stained macropores, water contents, andhydraulic conductivities. Monte Carlo parameters from theﬁnal set of acceptable model simulations were used as input
for model validation. The model reproduced the sharp
groundwater level rise and an instant drain discharge, afterstorm events and due to PF through shrinkage cracks. An
accurate characterization of bypass ﬂow after short-duration,
high-intensity rain events required a short time interval(b15 min) for rainfall measurement (Booltink, 1994).
In a more recent example, the effect of rain intensity on
the fraction of different macropores involved in PF was
evaluated in a stochastic application of the DPM IN3M (Weiler,
2005). A sensitivity analysis showed that ﬂow variability inmacropores at higher rainfall intensities markedly reduces
water transfer and increases bypass ﬂow. This might give anadditional, ﬁeld-scale explanation for the frequently observed
phenomenon of the small transfer rate coefﬁcient calibrated
in DPM to simulate non-equilibrium PF and transport (seediscussion in Jarvis, 2007). IN3M predicted water balance,water content variation, and dye coverage with depth as
observed in irrigation and dye tracer experiments at threeﬁeld sites reasonably well (Weiler, 2005).
Weiler and Naef (2003) simulated macropore ﬂow initiation at the surface or at a saturated subsurface layer.
Simulations were based on spatial distributions of earthworm
burrows and surface micro-topographies surveyed at four
grassland sites. A simple routing model of PF treated as matrix
saturation excess was applied. The results showed that

3.2.3. Summary for scale II (plot, ﬁeld)
Applied PF models simulating PF at the plot and ﬁeld scale
do not differ substantially from those used at the column
scale; they are still 1D and usually based on one- or twodomain concepts. The application time scale usually ranges
from weeks to months, whereas time intervals required for
accurate calculation of bypass storm ﬂow are still in the samerange as for scale I. Non-uniqueness problems were limiting
inverse PF parameter identiﬁcation at this scale. Some ofthese problems may be caused by the uncertainty in relating
local scale (microscopic) measurements to an effective(macroscopic) region of a two-domain model. Other limitations and problems observed already at scale I (see Section
3.1.3) are also transmitted to scale II. Stochastic PF modelstudies have shown potential to integrate effects of heterogeneity, parameter uncertainty, and PF triggering processes
into ﬁeld scale modelling. Macropore ﬂow variability was
found to reduce water transfer, which may serve as a new
explanation for the small transfer rate coefﬁcients in DPM
required to simulate PF and transport (see Section 4).
3.3. Scale III (hillslope, catchment area)
A multitude of hydrological processes operating atdifferent spatial and temporal sub-scales control the dynamic
catchment response. Catchment hydrology is a vast anddynamically developing research area far beyond the scope of
this review. Importantly, large catchment hydrological models address problems at a scale increasingly relevant to global
climate models. Interestingly, due to uncertainties involved instudying catchment hydrology, a paradigm shift appears to beunderway where traditional approaches are challenged, newphilosophical perspectives emerge, and innovative measurement and modelling approaches are proposed. The relations
between form and function, patterns and processes, genesis
and dynamics of catchments are moving into focus. The
reader is referred to an excellent overview of the state-of-theart in catchment hydrology research as covered in threeinvited commentaries (Tetzlaff et al., 2008; Soulsby et al.,2008; Dunn et al., 2008). Macropore ﬂow in the vadose zone
is but one out of several rapid ﬂow processes that may affectthe overall hydraulic response and water quality of a

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

catchment. A yet small but developing research area is on
utilizing models to analyze rapid ﬂow processes at thehillslope and catchment scales.
3.3.1. Hillslope scale
Many recent experimental investigations of hillslope
hydrology have considered the role of PF (e.g., Noguchi
et al., 1999; Buttle and McDonald, 2002; Tromp-van Meerveld
and McDonnell, 2006a,b; Lin et al., 2006).
3.3.1.1. Lateral preferential ﬂow: a threshold process.
LateralPF seems to be a hillslope characteristic responsible for peak
rapid runoff generation in humid catchment areas (e.g., Lin
and Zhou, 2008). Lateral PF is a complex process that can takea variety of forms including ﬂow in ‘pipes’ (i.e., large
macropores oriented parallel to the soil surface) along thebase of the soil proﬁle or in the topsoil (e.g., Uchida et al.,2001), or ﬂow through a dynamic network of PF paths
embedded in the soil matrix, movement in a thin saturated
layer or along micro-channels above bedrock, and transport in
exfoliation fractures in bedrock (Buttle and McDonald, 2002).
Evidence from hillslope studies suggests that although
individual macropore segments are often smaller thanapproximately 0.5 m in length, they have a tendency to selforganize into larger PF systems as sites become wetter(Noguchi et al., 1999; Sidle et al., 2001). There is evidence
that a threshold of total storm event precipitation has tobe surpassed in order to activate the lateral PF system andtrigger subsurface pipe ﬂow (e.g., Tromp-van Meerveld and
McDonnell, 2006a,b).
3.3.1.2. Hillslope PF models.
Different conceptual models
were proposed, including interconnected, self-organizing PF
systems in forested hillslopes (Sidle et al., 2001), hillslope pipe
ﬂow (Uchida et al., 2001), and ﬂow through a subsoil
macropore network as one of four main ﬂow pathways (Linet al., 2006).
Computational models include a 3D DPM that describes
lateral PF in a hillslope containing dry, partially ﬁlled, or fully
ﬁlled pipes in various arrangements. Matrix and macropore
ﬂows were computed using the 3D Richards and 1D
Manning's equations, respectively, coupled with an integral
expression for water transfer ﬂux. The model resultscompared well with bench scale experiments of Sidle et al.
(1995) (Tsutsumi et al., 2005).
Lehmann et al. (2007) applied a 2D percolation model that
represents the hillslope as a grid of sites forming a latticeconnected by bonds. Water-occupied sites connected by abond form a cluster. A cluster may grow until, at a certain
occupation probability (called the ‘percolation threshold’), it
spans the whole system and becomes a percolation cluster.
Model validation was attempted based on published hillslope
data of 123 events with measured rainstorm amount, outﬂow
and initial water content at 70 cm (Tromp-van Meerveld and
McDonnell, 2006a). Percolation theory appeared to be a valid
approach to quantify threshold processes controlling thehillslope outﬂow after rainstorm events (Lehmann et al.,2007).
A threshold of rainfall depth halting the rapid water tablerise and initiating lateral PF in hillslopes is also assumed in
modiﬁed DHSVM (Beckers and Alila, 2004).

15

3.3.2. Catchment area scale
3.3.2.1. Catchment area ≈ 4 km2.
Zehe and Blöschl (2004)
used CATFLOW (Zehe et al., 2001) to study how uncertainty of
initial soil moisture values is propagated into model prediction uncertainty of the hydraulic responses in the agricultural
Weiherbach catchment (3.6 km2), Germany. The model
parameterization was based on laboratory and catchment
data. The results showed that the uncertainty in model
predictions of discharge hydrographs was highest for initialsoil moistures around the threshold value separating matrix
and macropore ﬂows. The predictive uncertainty decreased
with increasing rainfall intensity. Moreover, satisfactory
discharge predictions were only obtained when using
spatially distributed precipitation values, not just an averagevalue. Finally, calculations were repeated for a plot within thecatchment. Perhaps surprisingly, the model uncertainty was
reduced from plot (CV = 0.21) to catchment (CV = 0.08). This
net reduction was explained by averaging of small-scalevariability; this countervailing effect exceeded additionalsources of variability at the catchment scale (Zehe andBlöschl, 2004).
Zehe et al. (2005) found that the simulated Weiherbach
runoff pattern depended on average initial soil moisture. By
contrast, spatial variability of initial soil moisture had amoderate effect (approximately +/−10% difference in runoff
peaks) only for dry or intermediate saturation, where
different realizations caused the saturation threshold to be
surpassed at different catchment locations, thus triggering
the two key processes of PF and overland ﬂow. Spatial
patterns of precipitation and macroporosity represented
additional sources of runoff prediction uncertainty (Zeheet al., 2005).
3.3.2.2. Catchment area ≈ 10 km2.
Beckers and Alila (2004)
modiﬁed (see Section 2.2.4) and applied DHSVM (Wigmosta et al., 1994) to analyze the contributions ofpreferential hillslope runoff to peak ﬂow generation in atemperate rain forest catchment (10 km2) on Vancouver
Island, British Columbia, Canada. Hydrometeorological datafrom 1972 to 1990 were utilized. DHSVM was applied toevaluate hillslope contribution to peak ﬂow generation in
the 10-km2 catchment area. Downslope subsurface ﬂow
rates during PF network formation were prescribed basedon three tracer tests. Soil depths between 0.5 m and 1.25 m
were deﬁned assuming Gaussian random ﬁelds with acorrelation length of 50 m. Soil hydraulic properties wereﬁtted by comparing observed and simulated streamﬂows at
four weirs in the catchment. The model was able tosimulate peak ﬂows and groundwater responses for both
low and high intensity rains. While the differences between
the ‘matrix-ﬂow-only’ and the ‘matrix-PF’ model
approaches for simulating low and medium streamﬂow
were subtle, higher values of model efﬁciency wereachieved for the ‘matrix-PF’ model of storm runoff. The
model evaluation suggested that a unit area discharge
threshold of 2.8 mm/h separates matrix ﬂow dominatedrunoff from PF dominated runoff. Matrix ﬂow dominatedbasin streamﬂow 67% of the year, while storm event ﬂows
were almost entirely derived from PF (Beckers and Alila,
2004).

16

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

3.3.2.3. Catchment area ≈ 40 km2.
MODHMS (Panday andHuyakorn, 2004) was applied with global inverse parameter
optimization to simulate spatially distributed tile drainage fora 2-year period in the Broadview Water District (38.8 km2),
San Joaquin Valley, California (Vrugt et al., 2004). The domain
was horizontally discretized into 200 m × 200 m square cells,
which were vertically subdivided into 14 layers of 0.3 m (top)
to 0.9 m (bottom). Inverse estimation of 7 parameters,
including initial and saturated water contents, bypass fraction, crop coefﬁcient, drain depth and conductance, and aMiller–Miller scaling factor for hydraulic conductivity, wascarried out. The simulation results of MODHMS werecompared with those obtained using a simple storage-based
bucket model and a spatially averaged 1D ﬂow model. The
study concluded that the drain conductance and the bypass
coefﬁcient were determined well (while soil hydraulicproperties were not) and that these parameters dominatedthe hydrology of the Broadview Water District (Vrugt et al.,2004). However, the conclusions about PF remain uncertain,
since the simple bucket model gave similar results.
3.3.2.4. Catchment area ≈ 75 km2.
Jones et al. (2008)
presented an application of the InhM model (VanderKwaak
and Loague, 2001) to model hydrological processes in the
75 km2 Laurel Creek Watershed in Southern Ontario, Canada.
Subsurface ﬂow was approximated as a steady-state process.
Model calibration yielded good agreement between theobserved and computed river baseﬂow discharge and hydraulichead values in 50 wells. Subsequent transient overland ﬂow
simulations based on rainfall time series were utilized to predict
the observed discharge hydrographs. Moderate agreement was
obtained. The model was deemed to be feasible to analyze thespatially heterogeneous and temporally dynamic response of awatershed to rainfall events (Jones et al., 2008).
3.3.2.5. Catchment area ≈ 350 km2.
Modiﬁed DSFDM wasapplied to simulate river ﬂow in the temperate forested
Hikimi river basin (352 km2) in a mountain region in Japan(Mulungu et al., 2005). Ten model parameters were calibrated
against the river ﬂow in 1991 using hourly rainfall data. Of
these optimized parameters, the saturated hydraulic macropore conductivity was the most sensitive. By contrast, thewater transfer rate between the matrix and macropores was
not a sensitive parameter. After model calibration, river ﬂow
in three subbasins was predicted between 1992 and 1996.
Simulation results appeared satisfactory, except that peak
ﬂows were underestimated (Mulungu et al., 2005).
3.3.3. Summary for scale III (hillslope, catchment area)
Models operating at the catchment scale require aphysically sound integration of the surface, vadose zone,
and underlying aquifer systems. Preferable, although difﬁcult
to parameterize, is a fully coupled model scheme where thethree layers interact with each other. An additional requirement posed here was to include PF. Most of the above models
fulﬁl these requirements, including fully 3D model systems(MODHMS, InHM, HydroGeoSphere), and distributed models
with some empirical simpliﬁcations in their equations(DSFDM, DHSVM). CATFLOW was applied with free drainage
at the lower boundary of individual hillslope segments, andthus without groundwater system (Zehe and Blöschl, 2004).

This approach may still give valuable insights for relatively
small catchments where the focus is on fast runoff components, or for studying ﬁrst-order controls. Still other approaches, such as percolation network models, are rare.
Compared to the predominantly vertical PF at smaller scales,
overland ﬂow and lateral PF (e.g., pipe ﬂow) appear as important
rapid runoff components in hillslopes. A consensus is emerging
about a strongly nonlinear PF response to rainfall depth andintensity, with a threshold value triggering rapid runoff response
in a hillslope. The identiﬁcation of such thresholds for certain
types of hillslopes or catchments might be one way towards
simplifying the prediction of their hydrological behaviour.
Problems encountered in the identiﬁcation of bypass ortransfer rate parameters are not surprising, given that those
problems were already encountered at smaller scales. However, models could still be used to assess how uncertainty ofantecedent soil moisture and precipitation, and their spatial
variability, is propagated into model prediction uncertainty.
It seems to be an unresolved question if solving overland
or groundwater ﬂow equations for horizontal cell sizes of
100 m × 100 m, with average properties, can be considered aphysically-based approach. Likewise, the vertical discretization limit for solving Richards equation in the vadose zone canbe in the cm range (Vogel and Ippisch, 2008) and is also
difﬁcult to fulﬁl at catchment scale. Advanced numerical
techniques (Li et al., 2005; Park et al., 2008) may eventually
overcome these problems.
Nevertheless, considerable advances have been achieved.
Water ﬂow studies in highly instrumented hillslopes were
instructive for conceptual model development and evaluation.
Sophisticated catchment scale PF models were developed andﬁrst applications for the catchment scale have been reported. The
validation of rapid runoff model components appears challenging in catchment hydrology studies (i.e. studies that solely lookat water ﬂow not using solute or other additional information).
4. Modelling of tracer transport in preferential ﬂow systems
A major concern with PF is its effect on contaminant solutetransport. Studying conservative solute transport is often
regarded as a prerequisite for understanding reactive solute
displacement. Moreover, some agrochemicals, such as nitrate,
behave relatively conservatively under certain conditions.
Thus, preferential tracer transport was investigated innumerous model applications at different scales (Table 3).
The ﬂow mode, i.e., steady state vs. transient, may have a
distinct effect on PNE solute transport (Meyer-Windel et al.,
1999; Cote et al., 1999). Traditional column studies involve
transport under steady-state ﬂow conditions and are modelled using closed-form analytical solutions. However, waterﬂow in the vadose zone involves dynamic inﬁltration,
redistribution, and capillary rise processes. Thus, models ofsolute displacement in the vadose zone generally considertransient water ﬂow and require numerical solutions.
4.1. Scale I (soil column, lysimeter)
4.1.1. Constructed soil columns, steady state ﬂow
Similar to investigations of water ﬂow, studies of preferential solute transport involved repacked soil columns with
artiﬁcial PF paths. Due to the impact of lateral concentration

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

17

Table 3Model applications for preferential solute transport
Scale Model type; name
Scale IACMIM; n.a.
AC

MIM; n.a.

ACc
ACMIM; BIM (Bimodal advection
model)
SPM, DPM; VARSAT2D, 2DCHEM

AC

MIM, DPM, SPM; HYDRUS-1D

AC,
UC
UC

MIM; Single Spherical
Aggregate Model (SSA)
MIM; CXTFITUC

MIM; CXTFIT and MODFLOW/MT3D,
DFM; FRACTRAN

UC
UC
UC

SPM, MIM, STM; n.a.
DPM; n.a.
DP-MIM; 2-MIM

UC

TPM; CXTFIT + PFUCUCDFM (MIM); boundary layer (BL)
model
MIM, DPM (TPM); MURTUC

Temporal moments analysis (TMA)

UC

DPM, SPM; MACROUC

MIM, DPM, SPM; HYDRUS-1D

Scale IIS
SPM; n.a.
S

DPM, SPM; MACROS

Multiscale random cascade;
Universal Multifractal Model (UM)

S,P,F

MIM, DPM, SPM; HYDRUS-2D

P

Drainage network construction; n.a.

F

SPM, STM, MIM; n.a.

F

MIM, DPM, SPM; HYDRUS-1D

F

DPM; SWAP + PF component

F

DPM; (sub)surface trigger modelF

Convective lognormal transfer
function; n.a.
DPM; M-2D

F

Scale IIIH
SPM (distributed, DEM based); HillVi

Model dimension; ﬂow state; description

Authors of model application

1D; steady state, MIM for macropore with radial(Fickian or ﬁrst-order) diffusion into matrix
1D, steady state, MIM, ﬁrst-order(vs. Fickian diffusion) transfer term
1D, steady state, Langrangian mobile–immobile
advective transport with bimodal lognormal pdf
2D, Richards (2× for DPM), CDE (2× for DPM), f.-o.
transfer of water and solute
1D, Darcy (2× for DPM), CDE (2× for DPM), f.-o.
transfer of solute
1D, steady state, MIM, Fickian diffusion transfer,analytical solution
1D, steady state, MIM, f.-o. transfer, analytical solution.
Also: CDE, method of moments

Young and Ball (1997, 1998)

1D, steady-state, saturated, MIM (CXTFIT), 2D-Darcy,
CDE-matrix, cubic law, CDE-fracture (FRACTRAN);
3D-Darcy, CDE (MODFLOW/MT3D)
1D, steady-state, model comparison
1D, steady-state, Darcy (2×), CDE (2×)
1D, steady-state, CDE (2×)
mobile–mobile–immobile model, transfer
1D, steady-state, saturated, CDE (matrix), Poiseuille(macropore, laminar), Manning (macropore, turbulent)
1D, steady-state, kinematic wave advective transport inmacropore, diffusion equation for solute transfer
1D, CDE (2×), ﬁrst-order solute transfer with embedded
small-scale structure
1D, steady-state, statistical integral functions ofdimensionless BTC moments
1D, Richards, CDE (matrix); kinematic wave,
convection (macropores)
1D, Richards (1–2×), CDE (1–2×), f.-o. transfer ofwater and solute

3D, Richards, CDE, stochastic or deterministic descriptionof spatially distributed hydraulic properties
1D, Richards, CDE (matrix); kinematic wave,
convection (macropores)
1D, no ﬂow model, random cascade with multifractal
spatial distribution of hydraulic conductivity, e.g. toanalyze spatial dye distribution pattern
2D, Richards (2×), CDE (2×), f.-o. transfer of water andsolute2D, no ﬂow model, reconstruction of paths connecting
similar scaling factors of water retention curves.
3D, Richards or Darcy-Buckingham, CDE, stochastic or
deterministic spatially distributed properties
1D, Richards (1–2×), CDE (1–2×), ﬁrst-order transfer ofwater and solute
1D, Richards, CDE (matrix); empirical mass balance androuting of water and solute (cracks), ﬁrst-order waterand diffusive solute transfer; dynamic shrinkage cracks
0D, steady-state, two-domain capacity model for
preferential ﬂow and solute advection
1D and 2D, no ﬂow, convective lognormaltransfer function
2D Richards and CDE (matrix);
1D kinematic wave and convection (macropores)

Grifﬁoen (1998), Grifﬁoen et al. (1998)
Rosqvist et al. (2005)
Allaire et al. (2002a,b)
Zhang et al. (2004) (AC4)Cote et al. (1999)
Vanderborght et al. (1997) (UC2); Langner et al.
(1999), Kamra et al. (2001), Ersahin et al. (2002),
Seo and Lee (2005), Stagnitti et al. (2000, 2001)
Jørgensen et al. (2004)

Vanderborght et al. (2002, 2006)
Schwartz et al. (2000)
Kjaergaard et al. (2004)
Perret et al. (2000)
Wallach and Parlange (1998, 2000)
Gwo et al. (1998)
Stagnitti et al. (2000), Kamra et al. (2001)
Kätterer et al. (2001), Logsdon et al. (2002),
Larsbo et al. (2005), Larsbo and Jarvis (2006)b
Köhne et al. (2004a,b)

Vogel et al. 2006
Merdun and Quisenberry (2004)
Olsson et al. (2002, 2007)

Abbasi et al. (2003) (P), (2004) (F); Haws et al. (2004)(S2), (2005) (F); Köhne and Gerke (2005) (P); Köhne
et al. (2006a) (F)Deurer et al. (2003)
Vanderborght et al. (2006) (UC2)
Gerke and Köhne (2004), Köhne et al. (2006a,b)
Van Dam (2000)

Kohler et al. (2003)
Gaur et al. (2006)
Abbaspour et al. (2001)

Weiler and McDonnell (2006) (Hb), (2007); McGuirei) 1D Vertical unsaturated (preferential) ﬂow. Mixing
et al. (2007)
within cells, advection between unsaturated andsaturated zones and between cells; ii) quasi-3D, saturated
lateral ﬂow routing between grid cells(continued on next page)18

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Table 3 (continued)(continued)
Scale Model type; name

Model dimension; ﬂow state; description

Authors of model application

Scale III
W, P
W
W

W

DFM; FRAC3DVS

3D-Richards, CDE (matrix), modiﬁed 2D-Richards andCDE (parallel fractures)a
3D Darcy, 2D overland, 1D vertical PF: Richards (2×) +DPM (distributed); WEC-C
advection (2×)
SPM (distributed, DEM-and GIS2D Richards' and CDE; bimodal K(h) with linear
based); CATFLOW
increase above threshold represents macropore
fraction; 1D Saint Venant for surface runoff and drainage.
DPMa (distributed); Mike She/Daisy + 3D Darcy, 2D overland , 1D channels, 1D PF: Richards,
MACRO
CDE (matrix); kinematic wave, advection (macropores)

Van Der Hoven et al. (2002) (Wd)
Croton and Barry, 2001 (Wb)
Zehe et al. (2001)

Christiansen et al. (2004)

Scale I: artiﬁcial or repacked soil column (AC), undisturbed column (UC), lysimeter (L); scale II: soil proﬁle (S), plot (P), hillslope (H), ﬁeld (F); scale III: watershed(W), river basin (R); PF — preferential ﬂow, SPM (DPM, TPM) — single (dual, triple) permeability model, MIM — mobile immobile model, KW (KDW) — kinematic(dispersive) wave model, f — fast ﬂow region, m — matrix, f.-o. (s.-o.) — ﬁrst- (second-) order, SCM — stochastic convective model, STM — stream tube model, DFM— discrete fracture model, pdf — probability density function.
aIn the vadose zone part of the model, b“synthetic, virtual soil”; assumed (or statistically generated) model parameters, cmunicipal waste sample, dfoamed zeolite/
iron pellet material.

equilibration on solute transport in the vadose zone (Flühler
et al., 1996), a common theme was the effect of the solute
mass transfer term on PNE transport simulations. Steadystate ﬂow conditions are considered below. First, ﬁndings ofsome classical studies are reported.
4.1.1.1. Solute transfer term.
For well-deﬁned, uniformly
sized and shaped matrix spheres and rectangular blocks orhollow cylindrical macropores, exact solutions of solute diffusion based on Fick's second law can be employed as theirtransfer term in the MIM. For such geometries, equivalent shapefactors were derived for use in a ﬁrst-order rate model of masstransfer as a zero-dimensional approximation of the diffusion
equation (e.g., Rao et al., 1980a,b; >van Genuchten and Dalton,1986). The MIM with ﬁrst-order transfer term is less accurate
for transport in spherically aggregated systems than fortransport through cylindrical macropores (van Genuchten and
Dalton, 1986), and for conditions of severe PNE in sphericallyaggregated systems may be even less accurate than the CDEwith an equivalent dispersivity that lumps the effects of intraaggregate diffusion (Parker and Valocchi, 1986). Moreover,analytical solutions for the ﬁrst-order mass transfer term andFick's second law were compared to the cyclic solute transfer
into and out of an immobile region. Limitations of the ﬁrstorder concept were found for short time periods withinsufﬁcient equilibration between phases of inward and outward diffusion (Grifﬁoen, 1998).
Young and Ball (1997) compared the effects of input
boundaries and other experimental conditions for soil columns
with cylindrical macropores vs. spherical aggregates. The ﬁtted
ﬁrst-order rate solute mass transfer coefﬁcient, α, in the MIM
systematically varied with column conditions and asymptotically approached a constant value only for large dimensionlesstime. The concept of a non-dimensional residence time scale
was useful in characterizing these variations, particularly forcylindrical macropore geometry (Young and Ball, 1997).
In accordance with this result, the examination of
published MIM analyses of laboratory column data sets
revealed a dominant trend of linear variation of α withresidence time and solute advective velocity (Grifﬁoen et al.,
1998; Maraqa, 2001; Haggerty et al., 2004).
Alternatively to the ﬁrst-order model, exact diffusion
equations were utilized in MIM for studying lateral solute

diffusion. A spherical aggregate model based on analyticalsolutions of the radial form of Fick's second diffusion law wasused to estimate the effects of intra-aggregate concentration
gradients on solute transport through aggregated soil columns. The model could be used to explain the transfer
mechanisms that increased leaching by 20% in an intermittent(as compared to continuous) ﬂow regime (Cote et al., 1999). Ina cylindrical macropore system, a MIM with diffusion
equation could be used to determine the effective diffusion
coefﬁcient of soil surrounding a sand ﬁlled macropore (Youngand Ball, 1998).
4.1.2. Constructed soil columns, transient ﬂow
4.1.2.1. Solute transfer.
Ahuja et al. (1995) used RZWQM tostudy transport of a surface-applied Br− tracer under rainfall
in a soil column with a 3-mm artiﬁcial macropore. A good
match of the rapid Br− breakthrough resulted when themacropore water was assumed to mix with only 0.5 mm of
matrix solution, revealing restricted solute transfer throughthe compacted macropore walls (Ahuja et al., 1995).
Cey et al. (2006) applied 3D DFM HydroGeoSphere(Therrien et al., 2005) to simulate the effect of lateral masstransfer on inﬁltration and tracer transport into a hypotheticalsoil column containing a vertical fracture. Rapid water transfer
restricted PF and model results were more sensitive to thematrix hydraulic conductivity than to fracture properties (Ceyet al., 2006). The above examples show two extremes; highly
restricted transfer with rapid breakthrough (Ahuja et al., 1995)
vs. fast transfer resulting in restricted PF and high sensitivity ofmodel results to soil matrix properties (Cey et al., 2006).
Since lateral solute transfer is now governed by the interplayof advective (water transfer driven) and diffusive components,
scaling rules derived for steady-state conditions will not generally
apply for transient ﬂow conditions. Mass transfer under dynamic
conditions, with changing directions to and from matrix, mighteven compensate some of the bias of the ﬁrst-order term
introduced for steady-state conditions, although this potential
effect has not been analyzed yet. At least in view of otheruncertainties, the inherent inaccuracy of the ﬁrst-order term may
be tolerable from a practical point of view (Jarvis, 2007). However,improved shape factors (see below, Section 4.1.3) and secondorder terms (Köhne et al., 2004a) were also tested.

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

4.1.2.2. Macropore arrangement and its conceptualization.
The effect of varying the continuity and tortuosity of(artiﬁcial) macropores on measured tracer transport andtransient ﬂow through soil columns was analyzed by Allaire
et al. (2002a) using the Richards ﬂow model VARSAT2D(Nieber et al., 1994) and CDE based 2DCHEM (Nieber andMisra, 1995). Both models were coupled. The PF region in
2DCHEM was deactivated; instead, the known macropore
structure and matrix were directly implemented in the 2Dmodel grid. Simulations showed that surface connected
macropores favoured the rapid transport of solutes appliedin solution, but delayed solute transport if the solute wasinitially incorporated in the soil. The latter result was
conﬁrmed repeatedly also for ﬁeld conditions (e.g., Larsson
and Jarvis, 1999a,b; Jarvis, 2007). Furthermore, several
tortuous macropores could be lumped into a single straighteffective macropore for simulations of the efﬂuent tracer BTC.
The importance of macropore continuity and tortuosity
increased with solute retardation (Allaire et al., 2002a).
In a companion paper, Allaire et al. (2002b) further
compared the 2D model explicitly describing the macropore
geometry with models considering various simpliﬁcations,such as a single straight ‘effective’ macropore, an SPM withbimodal hydraulic functions, and a DPM with the ﬁrst-ordertransfer term. The DPM was the best of the three simpliﬁed
approaches for predicting breakthrough curves during
steady-state ﬂow. For transient ﬂow (inﬁltration), differences
between the BTCs predicted by the different approaches
strongly decreased due to physical equilibrium caused byadevctive lateral solute transfer (Allaire et al., 2002b).
The acceleration of solute transport by macropores
reported above and in many other studies is environmentally
critical only when surface applied agrochemicals are exposed
to a rain storm following shortly after their application to thesoil surface. For other conditions, PF bypassing the matrix
may even slow down vertical leaching of solutes residing (andbeing ‘protected’) in the matrix (Larsson et al., 1999b; Köhne
et al., 2006b; Jarvis, 2007).
4.1.3. Intact soil columns, steady-state ﬂow
Studies discussed in this section addressed the questions
how natural soil structure can be represented in a two-region
model, what are the effects of different saturations on PNE
tracer transport at steady-state ﬂow, what controls parameter
identiﬁability in two-region models, and which alternative(other than two-region) model approaches can be used.
4.1.3.1. Embedding soil structure into shape factors.
Asmentioned above, the ﬁrst-order transfer rate coefﬁcient usuallycontains a shape factor representing some idealized geometry(e.g. spheres, cylinders). However, natural soils have a more
complex structure. Gwo et al. (1998) embedded the aggregate
scale structure and pore scale hydrodynamics into the ﬁrstorder mass transfer term by integrating solute ﬂuxes overmixed-size parallel slab matrix blocks using closed form
solutions of diffusion equations. Their time-independent shapefactor was calculated by summing up solute ﬂuxes over matrix
blocks of different categories (shape, size, water content andeffective diffusion coefﬁcient) with different probabilities. The
transfer term was implemented into a DPM version of MURT
(Gwo et al., 1995). Forward simulations of tracer BTCs from an

19aggregated saprolite soil column were conducted, assuming
aggregates grouped into 6 size categories of 0.3 to 18 mm radii,
with different observed frequencies. The BTCs could be
predicted for different transport velocities. Intra-region dispersivities were an order of magnitude smaller than those obtainedby ﬁtting a DFM or MIM to the BTC (Gwo et al., 1998).
A multi-rate solute transfer model was developed byHaggerty and Gorelick (1995). A MIM with this transfer modelyielded signiﬁcantly improved simulations of measured BTCs
in unsaturated soil columns as compared to a MIM with ﬁrstorder transfer. The assumed lognormal distribution oftransfer rates represented different time scales of diffusion
and required one additional parameter (Haggerty and
Gorelick, 1998).
The surface-area-to-volume ratio of the matrix can beclosely linked to shape factors (Gerke and van Genuchten,
1996). Diffusion into the matrix of a dual-porosity soil was
described by an equivalent cylinder system which preserves
the surface-area-to-volume ratio as derived from image
analysis of soil thin sections (Rappoldt and Verhagen, 1999).
The software package CXTFIT (Toride et al., 1999) withanalytical solutions for MIM and inverse parameter estimation procedures was frequently used to describe tracer
miscible displacement experiments, as for example in thefollowing studies also described below (Vanderborght et al.,1997, 2000; Langner et al., 1999; Wallach and Parlange, 2000;
Kamra et al., 2001; Ersahin et al., 2002; Seo and Lee, 2005).
4.1.3.2. Parameter identiﬁability in two-region models.
The
MIM mathematically degenerates to the CDE when concentrations in mobile and immobile zones are at equilibrium (Parker
and Valocchi, 1986). For such conditions BTCs described by CDE
and MIM become identical (Langner et al., 1999). The success ofparameter identiﬁcation may depend on different factors. The
extent of PNE directly affects the identiﬁability of MIM
parameters, as Vanderborght et al. (1997) demonstrated byanalyzing asymptotic parameter conﬁdence intervals derived
from sensitivity matrices. The MIM parameter identiﬁcation
depends on the solute application mode (pulse input is better
than continuous input), soil structure type (macropores are
better than aggregates), and types of data errors (constantcoefﬁcient of variation is better than constant variance)(Vanderborght et al., 1997).
Similarly, the extent of PNE also affects the identiﬁability of
DPM parameters. Physically unrealistic ﬁtted values of thesolute transfer coefﬁcient and the fraction of sorption sites wereobtained for Br efﬂuent BTCs from structured clay soil columns.
The almost similar performance of ﬁtted CDE and DPM
suggested that the transport regime was close to equilibrium.
Response surface analysis was used to analyze the problems ininverse parameter identiﬁcation (Schwartz et al., 2000).
Comegna et al. (2001) conducted seventeen leachingexperiments with three intact soil columns under bothsaturated and unsaturated conditions. A CDE, MIM, and CLT(convective lognormal transfer function) models all gave
indistinguishable results when applied to ﬁt the efﬂuent BTCs,
even those with an asymmetrical shape. However, ﬁtted
values of the retardation factor R b 1 in the CDE might suggest
the presence of PNE, in particular since the MIM could be usedwith R = 1 without decrease in goodness-of-ﬁt (Comegna
et al., 2001).

20

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Moreover, non-local techniques such as TDR for measuringresident solute concentrations usually fail to reveal PNE.
Therefore, ﬂux concentrations of efﬂuent BTCs were mostly
used to estimate PNE related model parameters (Vanderborght
et al., 2000). On the other hand, the MIM parameters for thenitrate BTC measured with a vertical TDR probe were
purportedly similar to those obtained from the efﬂuent BTC(Seo and Lee, 2005). This could indicate that vertical TDR probes
may be a suitable tool to sense a mixture of resident and ﬂux
concentrations in vertical PF paths and matrix.
4.1.3.3. Saturation dependence of PNE.
Ersahin et al. (2002)
analyzed bromide efﬂuent BTCs obtained from columns
collected from different horizons of a silt loam soil. For
saturated ﬂow conditions, columns with larger macroporosity
had larger values of pore water velocity, dispersivity, anddiffusive mass exchange rate coefﬁcient, but smaller values ofthe mobile water fraction than more homogeneous columns.
These differences were reduced with decreasing matric heads
between 0 and −5 cm, with fairly constant values between−5 cm and −10 cm. Thus, macropore ﬂow and PNE, as reﬂectedby asymmetric BTCs, disappeared with decreasing saturation(Ersahin et al., 2002). Similarly, the PNE of structured soilcolumns disappeared below the −5 cm matric head, where CDE
and MIM simulations of tracer (tritiated water and pentaﬂuorobenzoic acid) BTCs became alike (Langner et al., 1999).
Similar observations were made also by Pot et al. (2005),
who needed increasingly complex models (from CDE, MIM, to
DP-MIM; all as part of HYDRUS-1D) to describe tracer
displacement column experiments for increasing ﬂuxes andsaturations, and by Schwartz et al. (2000), who could describe
tracer breakthrough at −10 cm pressure head using the CDE,
while at −5 and 0 cm the ﬁt improved when using a DPM.
Evidence from these and numerous other studies suggeststhat PNE in soils is generated at pressure heads larger than
− 10 to − 6 cm in macropores of ‘equivalent cylindrical
diameter’ larger than 0.3–0.5 mm (Jarvis, 2007).
4.1.3.4. Alternative model approaches.
Different otherapproaches, such as temporal moments, DFM, CLT (see above,
Comegna et al., 2001), and three-region approaches werecompared with MIM to assess their individual weaknesses andstrengths. Stagnitti et al. (2000) found the statistical temporal
moments method to be as suitable as MIM for describing tracer
BTCs, while in another study, the moment method failed tocapture preferential peak concentrations of BTCs (Kamra et al.,2001).
A DFM that uses a KW equation to describe advective
transport in the fracture, i.e., a ﬁrst-order process to describe
dissolved chemical transfer through a boundary layer and adiffusion equation to describe the solute diffusion into thematrix, was successfully ﬁtted to BTCs measured for fractured
clay–till soils. By contrast, a corresponding model with local
concentration equilibrium at the fracture–matrix interface(without boundary layer) underestimated chemical breakthrough during early stages (Wallach and Parlange, 2000).
Jørgensen et al. (2004) compared the suitability of SPM,
MIM and DFM approaches for assessing Br transport insaturated columns (0.5 m × 0.5 m) of clayey till with vertical
fractures and root channel biopores. Simulations using
FRACTRAN (Sudicky and McLaren, 1992) in the SPM modeand MIM (Toride et al., 1999) required ﬁtted parameter values
that varied several orders of magnitude for different ﬂow
rates. Only the FRACTRAN in DFM mode could be applied withconstant calibrated values of fracture spacing and aperture toapproximate the observed BTCs during alternating ﬂow rates(Jørgensen et al., 2004).
DFM can simulate PNE transport at steady-state ﬂow moreaccurately than MIM with ﬁrst-order transfer. For fractured
soils with aggregate coatings or organic linings, DFM withboundary layer further enhances the accuracy.
For some soils, assuming two regions is not sufﬁcientfor transport characterization. A MIM and a DP-MIM (Fig. 1)
were compared for analyzing BTCs of tritiated water andcolloids, measured in 54 soil columns with different clay
contents under steady-state ﬂow conditions with a 5-cm
applied suction (Kjaergaard et al., 2004). Only the threeregion model could ﬁt all BTCs and supported the hypothesisthat there were two mobile and a diffusion domain in theclayey soils (Kjaergaard et al., 2004). The physical reality of
postulated ﬂow regions could be supported by combining Xray CT scans of macropores with measured tracer breakthrough to parameterize a three-mobile-region model (Perret
et al., 2000).
4.1.4. Intact soil columns and lysimeters, transient ﬂow
Modelling preferential solute transport under transient
ﬂow conditions requires a large number of model parametersthat usually cannot be measured. Some model packages
including MACRO (e.g., Jarvis, 1994; Larsbo et al., 2005) andHYDRUS-1D (Šimůnek et al., 2003, 2008a,b) have increasingly
come to be used with inverse and uncertainty estimation
techniques. Applications of these two models are reported
below as examples of other related approaches.
MACRO was coupled with the inverse parameter estimation program ‘Sequential Uncertainty Fitting’, SUFI (Abbaspour et al., 1997), for analyzing Br and Cl BTCs observed on
intact soil columns. The SUFI procedure starts with userdeﬁned prior uncertainty domains for parameters to be
ﬁtted. Each uncertainty domain is divided into equidistantstrata and parameter values are deﬁned by the ﬁrst momentof each stratum (Roulier and Jarvis, 2003). MACRO-SUFI wasapplied to study the effects of the initial saturation andthe mode of tracer application (in soil vs. in irrigation water).
A single ﬁtted parameter set could be applied to simulate
almost all elution curves, irrespective of initial saturation
and application mode, reasonably well (Kätterer et al.,2001).
The ‘Generalized Likelihood Uncertainty Estimation’ (GLUE)(Beven and Binley, 1992; Beven and Freer, 2001) was also
coupled with MACRO (Larsbo et al., 2005). The motivating
theory behind GLUE assumes that all models and measurements are to some extent wrong, so that many parameter sets
may represent the observations equally well. The GLUE
procedure runs multiple MACRO simulations with sets ofparameters sampled from predeﬁned distributions within
speciﬁed uncertainty limits. Simulations with model efﬁciencies larger than the predeﬁned threshold value are considered
acceptable. For a Cl displacement experiment, water percolation and Cl concentrations in both the efﬂuent and ﬁnal resident
proﬁles of intact soil columns were accurately reproduced by
MACRO with GLUE. Model efﬁciencies for acceptable parameter

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

sets can be used to recognize, which parameters can beconditioned (or inversely identiﬁed) in a particular transport
experiment (Larsbo et al., 2005).
Yet another inverse approach, called ‘parameter identiﬁcation method using the localisation of information’(PIMLI) (Vrugt et al., 2001), was applied with MACRO to
test the effect of the column design on parameter identiﬁcation (Larsbo and Jarvis, 2006). PIMLI identiﬁes the subsets ofdata that contain the most information for the identiﬁcationof each parameter. These subsets are then used iteratively to
reduce parameter uncertainty. The synthetic data representing transport experiments on clay and loam soil columns
contained sufﬁcient information to reduce uncertainty inthe mass transfer coefﬁcient, saturated matrix hydraulicconductivity and dispersivity. The information content waslargest in the early stages of the experiments, such that hightime-resolution measurements were suggested for the ﬁrst
irrigation(s) following solute application (Larsbo and Jarvis,
2006).
HYDRUS-1D (Šimůnek et al., 2005) includes the Levenberg–Marquardt method for multi-objective inverse parameter estimation. Köhne et al. (2004b) used the MIM
in HYDRUS-1D for inverse simulation of pressure heads,
water contents, water outﬂow, and Br breakthrough in the
efﬂuent of six aggregated soil columns with different initialwater contents subject to intermittent irrigations. Inverse
identiﬁcation of the required MIM parameters was fairly
successful, except for the saturated water contents in mobile
and immobile regions, which were highly correlated (Köhne
et al., 2004b). Using DPM and DP-MIM in HYDRUS-1D could
further improve model performance (Köhne et al., 2006a).
4.1.5. Summary for scale I (column, lysimeter)
For soil columns with artiﬁcial PF paths, and for naturalundisturbed columns, the effects of experimental conditions,
including the tracer application mode, ﬂow mode (steadystate and transient), saturation, ﬂow velocity, and residencetime, on PNE, solute transfer, and parameter identiﬁability,
were investigated. Steady-state and transient ﬂow conditions
required analytical and numerical solutions of PNE models,
respectively. Steady-state ﬂow conditions were frequently
used for investigations at scale I (less so at larger scales).
Model approaches used for both steady-state and transientconditions involved SPM, MIM, DPM, three-domain models,
DFM, and statistical method of moments. DFM (withboundary layer) was the most accurate approach for steadystate conditions. SPM, even if provided with bimodal
hydraulic functions, generally failed to simulate preferential
solute transport. The success of the other approaches often
increased with the level of model complexity, such as adding
another ﬂow domain, but this approach was rarely independently conﬁrmed using CT. However, model success was also
dependent on the parameter identiﬁability, as related toexperimental conditions, types of data available, and theextent of PNE. Signiﬁcant PNE was found for pressure heads
above −10 cm (see also Jarvis, 2007). For transient ﬂow,
different procedures were applied for inverse estimation ofparameters and uncertainty (Levenberg–Marquardt, GLUE,SUFI, PIMLI). Fitted parameters that are assumed to reﬂect
constant soil properties in the conceptual model, such as
mobile region, dispersivity, and transfer rate coefﬁcient, may

21vary with ﬂow rate and saturation if the model does notaccurately describe involved processes. However, it is notclear to which extent these effects are relevant for transient
vadose zone conditions. Structure related shape factors weresuggested as a means of including complex natural soil
structure in the solute transfer term.
4.2. Scale II (plot, ﬁeld)
Deterministic 1D and 2D two-region models with single
parameter sets were evaluated (compare Section 3.2.1.) andcompared with corresponding one-region continuum
approaches at the plot and ﬁeld scales to assess the relevanceof considering PF. Similar as for soil water ﬂow studies at thisscale (Section 3.2.), tile-drained ﬁelds were often involved in
Sections 4.2.1–4.2.2. A few alternative model applications aredescribed in Section 4.2.3.
4.2.1. 1D two-region approaches
4.2.1.1. Capillary PF models
Gerke and Köhne (2004) compared the DPM of Gerkeand van Genuchten (1993) with a corresponding RichardsCDE based SPM (Vogel et al., 1996) for simulating tile draindischarge and Br concentrations in a 0.5-ha ﬁeld with astructured loamy soil. The SPM and DPM described thedischarge hydrograph similarly well, but only the DPM
reproduced the Br breakthrough pattern. Hence, the DPM
was deemed to be capable of capturing relevant PNE transportmechanisms during transient PF (Gerke and Köhne, 2004).
Köhne et al. (2006a) studied the feasibility of the inverse(Levenberg–Marquardt) identiﬁcation of DPM parametersfrom a drainage hydrograph. The DPM implemented in
HYDRUS-1D was used to ﬁt hydraulic and transport parameters, either sequentially or simultaneously, using observed
tile-drainage hydrographs and Br concentrations. Only thesimultaneous ﬁtting procedure was successful in describing
Br breakthrough. From these and lab-scale results it wasinferred that a hydrograph alone is insufﬁcient for the inverse
identiﬁcation of soil hydraulic DPM parameters (Köhne et al.,2006a).
4.2.1.2. Gravity-driven PF models
Merdun and Quisenberry (2004) applied MACRO for
simulating drainage ﬂow and Cl transport in structured loamy
sand soil proﬁles at 6 plots. MACRO outperformed a onedomain model that over-predicted water content and Cl
concentrations with depth and time (Merdun and Quisenberry, 2004).
Similarly, MACRO gave a satisfactory description of soil
water contents, drainﬂow, and resident and ﬂux concentrations of Br observed in the 0.4-ha Lanna ﬁeld site (Sweden)with structured clay soil, when all available data were usedin the inverse parameter estimation (Larsbo and Jarvis,
2005).
4.2.1.3. Empirical PF models
Kohler et al. (2003) successfully applied their Surface
and Subsurface Trigger DPM to explain tile discharge and Br
ﬂux BTCs observed in a 2-year tile-drained ﬁeld experimentnear Zurich, Switzerland. Results were used to assess matrix

22

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

and PF components. Results suggested that 73% of Br was
exported in PF paths within 2 years after the applications, andthat the storage volume of the top soil had to be ﬁlled ﬁrst to
trigger PF in the subsoil (Kohler et al., 2003).
The empirical DPM SWAP was applied to ﬁeld observationsof soil water contents and Br concentrations in cracked clay soil,
groundwater levels, and ﬂuxes of tile-drainage water and Br
during 572 days. A single domain model was not able to mimic
the ﬁeld-average water ﬂow and solute transport. Routing of PF
through cracks considerably improved the SWAP simulation ofwater contents and Br leaching to groundwater. However, theobserved retardation and high dispersion of Br could still not be
reproduced (van Dam, 2000).
4.2.2. 2D two-region approaches
While lateral water ﬂow, e.g., in tile drained soils, can be
captured as a sink term to a 1D model, lateral solute transport
is problematic to model in a 1D framework. Hence, 2D tworegion models were developed and evaluated to model ﬁeld
scale preferential solute transport.
Kohler et al. (2001) evaluated their combined 2D/1D DPM
named M-2D. For a test data set based on a column
experiment with preferential tracer transport, M-2D performed much better than SWMS-2D (its 2D component), andsimilarly to MACRO (its 1D component). A ﬁeld tracer test
was conducted in a tile-drained municipal waste site in
Switzerland with measurements of tile discharge andelectric conductivity. M-2D approximated experimental
observations during both calibration and validation time
periods (Kohler et al., 2001).
Modiﬁed HYDRUS-2D model was used in several studies
comparing 2D-approaches based on SPM, MIM, or DPM for
simulating PNE tracer transport in agricultural ﬁelds (e.g.,
Abbasi et al., 2003, 2004; Köhne and Gerke, 2005; Haws andRao, 2004; Haws et al., 2005; Köhne et al., 2006b).
Abbasi et al. (2003) compared the 2D MIM and SPM using
inverse simulations of water ﬂow and tracer movement at a
sandy loam ﬁeld plot with furrows (3 m × 3 m) in Phoenix,
Arizona. The topography of the furrow cross-sections wasconsidered explicitly in the numerical grid. Water contents,
inﬁltration rates, and resident solute concentrations wereused in the objective function. The similarity of the inverse 2D
MIM and SPM results led to the conclusion that equilibriumtransport prevailed at this ﬁeld site (Abbasi et al., 2003).
Regarding this conclusion one must caution that the spatial
resolution of concentration data sampled at discrete positions
is usually insufﬁcient to detect PF (e.g., Kasteel et al., 2005).
Abbasi et al. (2004) extended investigations to a ﬁeld
(115 m × 20 m) covering several sets of three furrows each.
Fitted hydraulic and transport parameters at the ﬁeld scale
were more variable, yet relatively similar on average, to theplot-scale study (Abbasi et al., 2004).
At the Bokhorst loamy soil site in Northern Germany, theﬁnal Br distribution derived from multiple soil cores in a soil
proﬁle was described almost similarly by the 2D MIM and 2D
SPM, and thus reﬂected matrix transport. Br concentration
peaks in tile efﬂuent were only reproduced by the 2D MIM
and thus revealed preferential displacement. Nonequilibrium
was highly dynamic; it was limited to periods of highintensity rainfall suggesting a threshold process (Köhne andGerke, 2005). Simulations were further reﬁned using a 2D

DPM (S_2D_DUAL, Vogel et al., 2000) and improved boundary
conditions (Gerke et al., 2007).
At the Infeld site in North-West Germany, the observed
rapid Br efﬂuent breakthrough at low concentrations could
also only be simulated using the 2D MIM approach. Simulation results suggested that over 60% of the surface applied Br
was immobilized by transfer into the stagnant soil water
region, and that the 2D ﬂow ﬁeld induced by tile drainsenhanced Br dispersion (Köhne et al., 2006b).
Haws et al. (2005) simulated water and solute ﬂuxes fromsubsurface drains of two macroporous silty clay loam plots(48.5 m × 60 m) in agricultural ﬁelds in West Lafayette, Indiana.
The 2D SPM and 2D MIM simulations matched the drainagehydrographs, but not solute breakthrough. They concluded thata hydrograph ﬁt cannot be used as evidence of the physicalmeaning of model parameters (Haws et al., 2005). Similar
conclusions were drawn in other studies (Gerke and Köhne,2004; Köhne et al., 2006a; McGuire et al., 2007).
Another advantage of a 2D (vs. 1D) model is that soil spatial
heterogeneity can be considered in addition to PF (Vogel et al.,2000).
4.2.3. Alternative 2D or 3D approaches: fractals, networks, scaleway
Olsson et al. (2002) used a multiscaling approach todescribe the spatial distribution pattern of PF as visualized bydye tracer tests. The Universal Multifractal Model was basedon the random cascade process with a multifractal spatial
distribution of hydraulic conductivity. The model reproduced
hydraulic conductivity functions and key features in theobserved dye distribution, and was suitable for describing PF
pathways (Olsson et al., 2002).
The drainage network approach by Deurer et al. (2003)
was conceived for PF in hydrophobic sandy soils, but might
apply to structured soils as well. In three proﬁles (10 m length,
1.4 m depth) at a gleyic sandy Podsol in Northern Germany,
ﬂow path networks were identiﬁed by spatially linking thesimilar linear scaling factors of water retention curves. The
network fraction decreased exponentially with depth, corresponding to a strengthening of PF. Assuming 1D piston-ﬂow
through the effective transport volume of the network, thearrival times of the tracer peak concentration observed in a Br
tracer experiment could be predicted (Deurer et al., 2003).
This concept of an exponential decrease of the macropore
fraction with depth was incorporated by Haws and Rao (2004)
in the MIM within HYDRUS-2D (Šimůnek et al., 1999) toinvestigate the effect on solute transport. A strong effect on the
BTC peak position in the soil proﬁle was found, but the efﬂuent
BTC could be approximated by a depth-averaged macropore
fraction (Haws and Rao, 2004). However, a more generalapproach would be to relate transport parameters in models to
the texture and observed morphological features of soil
horizons, instead of making a priori assumptions (see discussion in Jarvis, 2007).
Vogel and Roth (2003) proposed the ‘scaleway’ approach,which is a conceptual framework for the explicit considerationof the structure (heterogeneity) in modelling ﬂow andtransport at the scale of interest, where microscopic heterogeneities are averaged and replaced by effective descriptions.
The scaleway needs a representation of the structure, a processmodel at the scale of interest, and corresponding effective
material properties. In an application of the scaleway, theJ.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

solute BTC in an undisturbed soil (Orthic Luvisol) column could
be accurately predicted based on computer tomography (CT)and network modelling. This approach seems to be more
physically based than, e.g., statistical multi-fractal approaches.
Vogel et al. (2006) explored the use of three different
approaches to incorporating the relevant 3D-structure of agiven soil for modelling ﬂow and transport. They performed a
Brilliant Blue dye tracer experiment in a 2.5 m2 plot. A 3D
hydraulic structure model (Fig. 1k) was constructed by
combining (i) direct measurement of soil layers, (ii) correlated
anisotropic random ﬁelds to describe heterogeneity withinsoil layers, and (iii) a genetic earthworm burrow model basedon percolation theory. This hydraulic structure model was
implemented into a 3D Richards-CDE model. The direct
simulation of dye distribution in the soil compared reasonably well with observations. Hence, a rough representation of
the soil structure and hydraulic properties may be sufﬁcient
to approximate tracer transport (Vogel et al., 2006).
Gaur et al. (2006) applied 1D and 2D convective lognormaltransfer functions (CLT) to evaluate whether concentrations intile drain outﬂow can be predicted from spatially variable soil
surface transport properties determined by multiple (45) TDR
probes that recorded the electrical conductivity (EC) in a ﬁeld
plot (14 m × 14 m) above 1.1-m deep tile drains. Additionally,
the EC of the tile drainage ﬂow was measured. Only
predictions made using 2D (not 1D) CLT reproduced theobservations. It was concluded that the surface solutetransport properties determined by TDR can be combinedwith a 2D CLT transport model to make reasonable predictionsof tile ﬂux concentrations for the soil that is relatively uniformwith depth (Gaur et al., 2006).
4.2.4. Summary for scale II (plot, ﬁeld)
The predominant ﬂow mode at scale II was transient bothin the experimental and model analyses. Contrasting
approaches were applied to conceptualize the effects of soil
structure on tracer transport. Two-region (MIM, DPM) models
are the most widely used approaches at scale II for
characterizing preferential ﬂow and transport. Increasingly,
these PF models also consider transient variably-saturated
ﬂow in 2D, and are complimented with options for inverse
parameter identiﬁcation.
Model comparisons again demonstrated better performance of the two-region as compared to continuum models
for describing conservative solute transport, while relatively
similar results were obtained for the hydrograph. Furthermore,
two-region approaches were used for better understanding ofPF effects on solute transport. For example, preferential ﬂow
and transport can be triggered in the subsoil, and bypass ﬂow
may under certain conditions reduce tracer leaching. Alternative models included stochastic multifractal description,
drainage network reconstruction, a 3D hydraulic structure
model, drainage network reconstruction, and 1D or 2D
convective lognormal transfer function models (CLT).
Non-ideal transport may be caused by the soil structure (micro-heterogeneity), the spatial variability (macroheterogeneity) of soil properties, or a combination of both. A
decade ago Feyen et al. (1998) suggested that modelling
concepts for micro- and macro-heterogeneity should becoupled in one overall mathematical framework, since bothtypes of heterogeneity are relevant for agricultural and23

environmental management. A few recent approaches havepresented possible avenues in that direction (e.g., Vogel andRoth, 2003; Vogel et al., 2006).
4.3. Scale III (hillslope, catchment area)
Weiler and McDonnell (2004) proposed virtual (numerical) experiments to explore ﬁrst-order controls on complex
hillslope hydrology. They applied their quasi-3D, distributed
hillslope model HillVi. Results suggested that high drainable
porosity increased channeling of lateral ﬂow at the soil–
bedrock interface. Weiler and McDonnell (2007) additionallyprovided HillVi with lateral PF pipes. The model was applied
to predict runoff, and transport of Br applied as a line source
tracer, at a hillslope in the Maimai catchment area (New
Zealand). Pipe network conﬁguration had a much higher
inﬂuence on Br movement than on runoff. Although sometracer was rapidly displaced through the pipe system, thebulk of tracer remained in the soil. This explained the(frequently) observed dominance of pre-event water in runoff
events (Weiler and McDonnell, 2007).
Van Der Hoven et al. (2002) also investigated how toseparate storm hydrographs into pre-event water and stormevent water, but based on the δ18O natural isotope content.
While this is a problem relevant to catchment scale, theexperimental analysis of δ18O values in shallow groundwater
focussed on one multi-port well in fractured saprolite. The
DPM FRAC3DVS (Therrien and Sudicky, 1996) was applied toanalyze the data. Results suggested that the δ18O content was
not constant over the course of a storm event, and that mostof the matrix water was isolated from ﬂow in fractures. It was
concluded that quantitative hydrograph separation using
isotope techniques is not possible in most situations (Van
Der Hoven et al., 2002).
Christiansen et al. (2004) coupled the 1D vadose zone —3D groundwater hydrology model MIKE SHE operating at thecatchment scale (Refsgaard and Storm, 1995), the 1D agroecosystem model DAISY (Hansen et al., 1991; Abrahamsen and
Hansen, 2000), and a 1D DPM based on MACRO (Jarvis, 1994).
The resulting model was used to study the effects ofmacropore ﬂow on stream discharge and groundwater levels
in a small catchment (1.5 km2) in Zealand, Denmark.
Simulation results indicated an erratic and spatially variable
nature of macropore ﬂow depending in a complex manner onsoil characteristics and hydrological regime. Three years oftracer applications were not sufﬁcient to estimate averageleaching to groundwater, or to predict which rainfall events
generated the highest macropore ﬂows and leaching of tracer
or pesticide. Model validation for macropore ﬂow washampered by the lack of experimental data at this scale(Christiansen et al., 2004).
Zehe et al. (2001) assembled the CATFLOW model basedon experimental studies carried out in the Weiherbach
catchment (6.3 km2) in S-W Germany. CATFLOW integrates
three different spatial scales (plot, ﬁeld, catchment) and threedifferent time scales for BCs and parameters (constant,season, event). On different space-time scales, water andsolute dynamics are dominated by different processes: On theplot scale, ﬂow and transport are approximated as 1D verticalprocesses; inﬁltration (PF) and solute adsorption are relevanton the event time scale and evapotranspiration, interception,

24

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

and degradation dominate seasonal scale processes. On thehillslope (ﬁeld) scale, runoff generation and inﬁltration (PF)govern water balance on the event time scale, while on theseasonal scale, precipitation variation, land use changes andsoil fauna activity affect water and solute displacement. On
the catchment area scale, water dynamics depend on
topography, land use, and geometry of the drainage channels.
The CATFLOW simulation results were compared to observations of water contents, inﬁltration rates, and ﬂow patterns
from the Weiherbach catchment. Again, model validation was
limited by the non-representativeness of point data for larger
areas (Zehe et al., 2001).
According to Weiler and McDonnell (2007), conceptualization and parameterization of the effects of lateral PF on
hillslope hydrology currently represent the greatest challenge
in modelling PF processes at the catchment scale.
4.3.1. Summary for scale III (catchment area)
Model analysis of PF effects on solute transport is aformidable task at the catchment area scale, where PF is oneof the many processes governing solute delivery from thecatchment to the stream. Although several models have beendeveloped in the past decade, only a few of them were applied
to realistic situations, and consequently the effect of PF on
catchment hydrology and water quality is still not wellunderstood. A principal problem at the catchment area scale
is to assess and validate models based on experimental data:at the lower scale end, point measurements (e.g., soil watercontents) are measured at the model sub-grid scale and areusually scarce, while at the upper scale end, observations of
the integrated catchment response (e.g., river hydrograph)
are above the process scale represented in the model. The
validity of quantitative hydrograph separation using massbalance techniques for natural isotopes (e.g., δ18O) wasquestioned.
Among approaches that may be worth further exploration
are nested experimental tracer tests and models (plot-ﬁeldcatchment area), advanced measurement methods including
geophysical and remote sensing techniques, and stochastic
simulations and virtual experiments.
5. Discussion
5.1. Comparing strengths and weaknesses of two-region models
While several other promising PF approaches are being
developed and tested as described in previous sections, thetwo-region concept is now by far the most widely applied
approach to model PF in agricultural structured soils and thus
our discussion here will be focused on them.
Capillary, gravity driven, and empirical two-region
approaches differ in their complexity and model assumptions.
For model applications the question arises if any one of theseapproaches can be considered superior to others, if only for aparticular purpose. It would be desirable to base a model
ranking on an evaluation of reported model applications.
However,, this is difﬁcult for three reasons: (i) most studies,with few exceptions (e.g., Logsdon et al., 2002; Aden andDiekkrüger, 2000; van Dam, 2000), only report overall successof a model after calibration, (ii) calibration methods andcriteria for the reported success differ between the studies,and (iii) there are too few studies on model comparison basedon the same data set and furthermore, of those studies appear
to be of limited signiﬁcance due to the uncertainty introduced
by different model users. Thus in the following we discussmodels' strengths and weaknesses on the basis of underlying
model conceptualizations.
The following general features should be included in a DPM
to be applicable to a range of conditions: 1) physically based
process descriptions with parameters (at least in principle)
linked to measurable quantities, 2) options to consider soil
heterogeneity; at least soil layering, 3) transient ﬂow in both
domains at a high time resolution, 4) program-controlled,
domain-separated, dynamic upper BCs to conceptualize thedynamics of water and solute inﬁltration into matrix andmacropores at different rainfall rates, and other surface
processes such as surface runoff and evaporation, 5) lower
BCs realistically representing different situations including
variable groundwater table, gravity ﬂow, or drainage, 6) bidirectional, rate-limited transfer of water and solute, 7)capability of initiating PF at both surface and subsurface
layers, 8) inverse and/or uncertainty estimation procedures,and 9) a range of options required to describe agricultural
management options (e.g., regarding tillage, drainage, crops,
mechanical compaction, different modes of chemical application) and different agrochemicals.
Even recent empirical capacity PF models meet several ofthese requirements, e.g., the relatively comprehensive COUP
model (Jansson and Karlberg, 2004). However, capacity androuting approaches cannot describe PF for moderate rainfall
rates below the saturated soil hydraulic conductivity. Furthermore, they usually do not consider subsurface initiation of PF atlayer boundaries. Exceptions are the Subsurface Trigger model(Kohler et al., 2003) and the 1D stochastic IN3M model (Weiler,
2005). Lateral transfer between the macropores and matrix is
usually considered by layer wise mass balance without
physically based equation. Thus, empirical PF approaches are
suited for speciﬁc conditions only.
Currently, DPM approaches which fulﬁl most or all of theabove requirements, are either based on gravity-driven
macropore ﬂow, such as MACRO (Larsbo et al., 2005), or
assume that ﬂow through the soil structure network can be
described similarly as ﬂow in a porous-medium, which iscontrolled by capillary and gravity forces, such as HYDRUS-1D(Šimůnek et al., 2008a,b) (Section 2.).
The most obvious conceptual difference between PF
models concerns the question of an adequate description of
bypass ﬂow. A recent detailed discussion (Jarvis, 2007) of ﬂow
processes in macropores based on ﬁndings from various
experimental and modelling studies concludes that PNE insoils is usually generated at pressure heads larger than −10 to−6 cm in macropores of a minimum equivalent diameter of
0.3–0.5 mm, and that gravity is the main driving force for PF.
This suggests the choice of a pressure head of −10 cm
separating matrix from macropore regions. While DARCY's
law might also be accepted as a reasonable quasi-empiricalmodel of macropore ﬂow, the gravity-driven KW approachshould give a (quasi-) physically based description of PF instructured soils (Jarvis, 2007). However, the advantage ofusing a gravity-driven over a capillary approach might not beuniversally given for all structured soils, conditions, andscales, as discussed below.

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

5.1.1. Microscopic scale
Current conceptual ideas about pore scale fracture ﬂow
processes appear to be largely based on laboratory studies
performed with low permeability fracture surfaces (e.g.,Tokunaga et al., 1997; Or and Tuller, 2000; Ghezzehei, 2004;
Dragila and Weisbrod, 2004). These processes may not be
representative of ﬂow in many macroporous soils, where thelarger matrix hydraulic conductivity causes faster waterexchange and slows down macroporous ﬂow. Even conclusions based on laboratory studies of macropore ﬂow velocities
in soils might not tell the whole story. For technical reasons,these studies were mostly performed for macropore diameters of several mm, where ﬁlm ﬂow was observed forunsaturated conditions and full ﬂow for saturated conditions(e.g., Logsdon, 1995; Ghodrati et al., 1999), compliant with thegravity-driven KW approach. Indeed, Jarvis (2007) alerted
against the possible bias introduced by the high irrigation
intensities of larger than 10 mm h− 1 used in the majority(60%) of the reviewed studies. By contrast, the lower end of PF,
i.e., ﬂow in smaller macropores or at lower ﬂux rates inducedby natural rainfall rates, is less well known. Signiﬁcant
preferential movement of NO3 occurred in soil columns
with earthworm burrows and root channels even at ﬂuxes
as low as 1.5% and 5% of the saturated hydraulic conductivityof the column, respectively (Li and Ghodrati, 1995, 1997). A CT
based analysis of ﬂow in macropores in three different soils
yielded, although uncertain, pressure head values of −40, −14,and −10 cm for the transition between matrix and macropore
ﬂow (Mori et al., 1999). Measurements of diffusion into soil
aggregates, and of BTCs in soil columns at different pressure
heads (0 to − 15 cm), suggested that with decreasing
saturation, transport through smaller macropores continues
at reduced degree of nonequilibrium (Haws et al., 2004). In
conclusion, limited experimental evidence does not seem toprove that capillary forces cannot affect PF in soils.
5.1.2. Macroscopic scale
Two-region models are not exactly based on pore scale
physics, but are macroscopic approaches (e.g., Jarvis, 2007).
Thus, the question is not how accurately the model represents
pore-scale processes, but how the model, supplied with
‘efﬁcient’ parameters, can describe some average processes
and phenomena. This question has been evaluated for STM(stream tube models) and CDE models (Vanderborght et al.,2006), but not for DPM. At least some qualitative considerations are attempted here. A characteristic of large soilcolumns, and in particular of an entire ﬁeld, is the heterogeneity of soil texture and structure, where the latter includes
fractures, root channels, and earthworm burrows of differingsize, directions and continuity, causing different transport
patterns in different soil horizons (e.g., Kulli et al., 2003). The
effect of structure variability on the hydraulic conductivity
function of the fast ﬂow region can be represented in thegravity-driven KW approach by the kinematic exponent, andin the capillary PF approach, e.g., by the tortuosity factor in
van Genuchten's (1980) model. Additionally, solute transport
within such a heterogeneous PF network will be highly
dispersive due to different local ﬂow velocities (Vogel et al.,2006). Dispersion within the structure network can be
represented in the capillary PF approach by using the CDE,
while in the gravity driven approach entirely by mass

25

exchange with the matrix. Moreover, a 2D capillary PF
approach can be used to simulate lateral PF.
5.1.3. Conclusion
In summary, the capillary PF approach appears to be auseful concept for the purpose of describing bypass ﬂow and
PNE transport in soils with hierarchical structure, e.g.,
structured loamy soils. For clay soils with open cracks, orsoils with large continuous macropores at the surface, thegravity-driven approach is better suited. Both gravity-driven
and capillary PF approaches are considered to be largely
physically-based for applications at the soil column scale, but
semi-physical when used at the ﬁeld scale, where thestructure of many soils is too complex to be classiﬁed into
exclusively non-capillary or capillary pores.
2D approaches were also used for both, the gravity-driven
approach (M-2D; Kohler et al., 2001) and the capillary PF
approach (e.g., HYDRUS (2D/3D) (Šimůnek et al., 2006b).
Adding a description of spatial variability of soil properties to
1D (Booltink, 1994; Weiler, 2005), or 2D (Vogel et al., 2000)
approaches would strengthen the physical basis for ﬁeld scale
applications. However, model applications involve striking abalance between parameterization effort and the required
accuracy. In that practical sense, both approaches have been
demonstrated to be useful in applications for different soils at
scales from the column to the ﬁeld.
5.2. Beneﬁts and limitations of inverse methods
In the past decade, the use of inverse techniques has
become a standard procedure to facilitate DPM applications(e.g., Section 4.1.4). Roulier and Jarvis (2003) argued that
inverse estimation procedures provide objective, reproducible, and unambiguous results, providing the problem is well
posed. This statement is certainly accurate, but, unfortunately,in DPM ﬁeld applications the problem is often not well posed.
This was shown in simulations of preferential water ﬂow
based on scarce data (Section 3). When modelling solutetransport, several parameters related to ﬂow, transport, andlateral transfer have to be estimated. To better constrain the
inverse solution, many studies considered different types of
ﬂow and transport data simultaneously in the objectivefunction. Two problems arise in this approach. (i) Dataobtained at a local scale, such as tensiometer or TDR readings
and concentration proﬁles, can only be represented in a
spatially resolved (2D or 3D) model. Additionally, in a tworegion approach which data type corresponds to which
domain must be deﬁned. (ii) Automated weighting procedures
cannot fully compensate the effects associated with different
data sets (with divergent mean and variance, and variable data‘noise’) in the objective function. Some manual ﬁne tuning of
weights or other factors is usually needed for optimal results(e.g., Scorza jr. et al. 2007, Köhne et al., 2004b). The uncertainty
involved in inverse procedures has given incentive to develop
uncertainty estimation procedures such as GLUE and SUFI
(Section 4.1.4.) and multi-objective optimization algorithms(Wöhling et al., 2008). Although inverse techniques havebecome very helpful tools, their further improvements areneeded for their objective use. This may include focussing on
subsets of data which are most important (Larsbo and Jarvis,
2006), or recognizing signiﬁcant patterns.

26

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

6. Independent parameter determination
6.1. Scale I (soil column, lysimeter)
6.1.1. CT, MRI, and tracer techniques
The X-ray CT scan technique was used to analyze the large
pore structure, while MRI was used to assess both the 3D porespace arrangement and ﬂow velocity. Progress was achievedby combining techniques including CT, MRI and solute or dye
tracer breakthrough to describe relations between porosity
distribution, ﬂow paths, and soil hydraulic and transport
properties (Table 4).
Several studies evaluated the usefulness of hydrogeophysical methods for independent model parameterization. For
example, Císlerová and Votrubová (2002) studied CT identiﬁcation of local high porosity zones. The CT scans identiﬁed
macropores in the undisturbed sandy loam. However, CT wasinsufﬁcient to separate fast and slow ﬂow domains in the ﬁnesand columns, where this supplementary information could
be provided by MRI and dye tracer distribution.
Perret et al. (1999) used CT for establishing the 3D
arrangement of the geometry and topology of macropores.
Perret et al. (2000) further applied CT for real-time examination of Potassium iodide movement through the macroporesand matrix of a large undisturbed soil column. Based on the
CT results, distribution density functions of macropore sizeand hydraulic radius were assigned to the two (laminar andturbulent ﬂow) macropore subdomains of a three-mobileregion model. The breakthrough curves measured in thematrix domain were ﬁtted using the CDE (Toride et al., 1999).
Measured potassium iodide breakthrough in the macropores
could then be predicted using the three-region model (Perret
et al., 2000).
Vanderborght et al. (2002) explained concentration
patterns of two adsorbing ﬂuorescent dye tracers in a soilcolumn with the large 3D pore structure derived from CT
imaging. Additionally, breakthrough of Cl and dyes in acolumn, and dye adsorption in batch tests, were measured.
The data were used together to explain the success of 1D MIM,
and the failure of SPM and STM, in describing observed
breakthrough of Cl and dyes. In MIM, the fraction of sorption
sites in the mobile region could only be inversely identiﬁed
when the mobile region fraction was independently determined from dye concentration maps. It was concluded thatthe combination of efﬂuent breakthrough data and patternsof local concentrations in the soil provided the information to
identify governing transport processes and model parameters(Vanderborght et al., 2002).
Kasteel et al. (2000) used the CT derived X-ray absorption
coefﬁcients, which are related to local bulk densities, as aproxy for hydraulic properties for two different ﬂow regions.
The relative hydraulic conductivity, water retention, andconnectivity functions of the fast ﬂow region were derived
using a network model based on 3D CT images of local
densities in a small undisturbed soil block (resolution
0.04 mm) (Vogel and Roth, 1998, 2001). These relative local
hydraulic functions were scaled using measured absolutevalues, and were then implemented as spatially distributed
functions in the 3D Richards-CDE model SWMS_3D(Šimůnek et al. 1995) for direct simulation of observed tracer
breakthrough. Tracer breakthrough was predicted nearly as

accurately as when MIM was ﬁtted to the data (Kasteel et al.,2002). Javaux et al. (2006) also considered the CT-derived 3D
microscale spatial distribution of hydraulic parameters in
SWMS_3D to reproduce tracer BTCs in a sandy monolith.
The MRI technique was applied for assessing the 3D porespace arrangement in a repacked soil column with sand andgravel layers (Baumann et al., 2000). Flow velocities within PF
paths and the diffusion coefﬁcient were measured using MRI
at a spatial resolution of 1.3 × 1.3 × 5 mm3 and at time intervalsof 10 s for different applied ﬂow velocities between 0.9 and40 cm/min (Baumann et al., 2000).
A 3D hydraulic tomography method, i.e., inversion of
travel times of transient pressure propagation measured inhydraulic and pneumatic tests in an unsaturated fractured
sandstone core, was used to estimate the local hydraulicproperty distribution and diffusivity (Brauchler et al., 2003).
Other hydrogeophysical techniques include self potential,
ground penetrating radar, induced polarization, electricalresistivity tomography, electromagnetic induction, microwave radiometry, and time domain reﬂectometry (TDR)(Vereecken et al., 2004).
6.1.2. Imbibition and diffusion in soil aggregates to estimate
mass transfer
Gerke and Köhne (2002) and Köhne et al. (2002a)
investigated the effect of aggregate coatings as an approximate estimator for water and solute mass transfer parameters. Aggregate coatings were found to reduce water
imbibition (Gerke and Köhne, 2002) and solute diffusion(Köhne et al., 2002a) into soil aggregates by approximately anorder of magnitude. Results were only a rough approximationof transfer mechanisms in the ﬁeld, since the effects of PF innatural soil position (partial wetting, hydraulic inactivity)
could not be considered.
6.1.3. Soil hydraulic measurements to derive PF-related soil
hydraulic functions
Hydraulic measurements on soil cores with PF paths wereﬁtted with bimodal (Coppola, 2000), dual-permeability(Köhne et al., 2002b), and multi-modal (Priesack and Durner,
2006) hydraulic functions. Physically-based dual-domain
hydraulic functions were developed by upscaling modelledliquid conﬁgurations in partially saturated fractured porous
media, and were ﬁtted to hydraulic data (Tuller and Or, 2002;
Or and Tuller, 2003). However, uncertainty will be introduced
for all such models by purely hydraulic ﬁtting (e.g., Köhne
et al., 2006a).
6.1.4. Column experiments and pedotransfer functions
Methods for estimating parameters of the steady-state
MIM involve parameter ﬁtting under optimized experimental
conditions. For example, the immobile water content, thesolute transfer coefﬁcient, and the dispersion coefﬁcient were
estimated using a vertical TDR method that was used toanalyze resident tracer breakthrough (Lee et al., 2000, 2001,
2002). A ﬂow interruption technique for intermittent soluteleaching through soil columns was applied to estimate themobile and immobile water contents, and the mass transfer
coefﬁcient (Ilsemann et al., 2002). Several other studies havereported similar analyses for the MIM (Shaw et al., 2000;
Ersahin et al., 2002).

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

27

Table 4
Independent parameter determination for preferential ﬂow and transport models
Scale Parameters (model)
Scale I
Ped
Diffusion coefﬁcients of ‘interface PFmatrix’
Ped
Hydraulic conductivity of ‘interface PFmatrix’AC
3D pore space arrangement (of repackedsoil), diffusion coefﬁcient, ﬂow velocity inpreferential ﬂow paths
AC

UC

UC

UC

UC

UC

UC

UC

UC
UC
UC

Determination approach

Diffusion measurement with intact and surface-removed
aggregates
Flow measurement with intact and surface-removed
aggregates
3D estimation of pore space (0.5 × 0.5 × 0.5 mm3 resolution),ﬂow velocity (1.3 × 1.3 × 5 mm3, time resolution 10 s), anddiffusion coefﬁcient in soil column using magnetic
resonance tomography (MRT).
Genetic-based computer search algorithm to invert scale3D equivalent preferential ﬂow paths in
representative equivalent fracture networks using BTCs andfractured porous media (ﬂow-transport
structural information
DFM)
3D pore space arrangement and its relation 3D visualization of pore space (resolution 0.5–1 mm) of dry
to ﬂow paths
undisturbed soil samples using X-ray computer tomography(CT) images
3D pore network and water retention and 3D scan (resolution 0.04 mm) of a small undisturbed soil
permeability functions of fractured porous block using X-ray CT images, reconstruction of 3D pore
network connectivity deﬁned by the 3D Euler number.
media
3D structure of local-scale (1-mm)
3D scan (resolution 0.5–1 mm) of an undisturbed soil
hydraulic properties
column using X-ray computer tomography (CT) images. Xray absorption coefﬁcients related to local ρ and K(θ) wereused in a 3D ﬂow-transport model to predict observed tracer
breakthrough3D geometry and topologya of macropore CAT scans and image analysis of four large (800 mm length,
networks (conceptual model of four ﬂow 77 mm diam.) intact dry soil columns, 3D reconstruction ofmacropore networks, solute breakthrough analysis in spacedomains)
and time
3D hydraulic property distribution and3D hydraulic tomography by inversion of travel times ofdiffusivity
transient pressure propagation and diffusivity measured inhydraulic and pneumatic tests in unsaturated fractured
sandstone core.
3D-CT scan: model parameters, e.g., θm, α, 3D X-ray CT scans combined with solute breakthroughanalysis for Cl and two ﬂuorescent tracers in three intact soil
λm, f (MIM), soil properties that control
columns: 3D pore structure, spatial dye distribution after
transport of adsorbing solutes
leaching, modelling (CDE, MIM, STM)
θ(h) dependency on micromorphology
Micromorphology of soil pore structure was studied on thinsections. Hydraulic properties were studied on undisturbed
soil samples using constant head and multi-step outﬂow
methods
θim, α, Dm (MIM)
TDR method analyzing resident concentration breakthroughof salt solution applied at step input.
Intermittent solute leaching and ﬂow interruption for soilθim, α (MIM)
columns
θ(h) and K(h) parameters (bimodal,
Parameter ﬁtting of retention and hydraulic conductivity
multimodal, or DPM)
measurements on soil coresUCθ(h) and K(h) of fractured porous media

UC

Ks, Kb hb, n⁎, Dv, d, mixing depth (MACRO
DPM)

UC

Ks, θim, α (MIM)UC

d (MACRO DPM)

Scale IIS
Ks, θim, α (MIM)S

K, θ, θim, α (MIM)

S

Ks, Kb hb, n⁎, θsma, d (DPM; MACRO)

Authors
Köhne et al. (2002a,b)
Gerke and Köhne (2002)
Baumann et al. (2000)

Gwo, (2001)

Císlerová and Votrubová et al. (2002)

Vogel and Roth (2001)

Kasteel et al. (2000)(Perret et al., 1999, 2000)

Brauchler et al. (2003)

Vanderborght et al. (2002)

Kodešová et al. (2006)(Lee et al., 2000, 2001, 2002)
Ilsemann et al. (2002)(Othmer et al., 1991; Ross and Smettem, 1993;
Durner, 1994; Coppola, 2000; Köhne et al.,
2002b, Priesack and Durner, 2006)
Upscaling of modelled liquid conﬁgurations of cross-sections (Tuller and Or, 2002; Or and Tuller, 2003)
of partially saturated fractured porous soil or rockMerdun and Quisenberry (2004)
Statistical variation in MACRO parameters calibrated insimulations of water ﬂow and solute transport between 3
soils with different textures
Statistical variation in MIM parameters with tracer BTCs, dye (Vervoort et al., 1999; Shaw et al., 2000)
stained soil fractions, and soil properties
Statistical variation in the mass transfer coefﬁcient for BTCs Jarvis et al. (2007)
in 33 intact soil columns representing a range of soil texturesand organic matter contents

Alletto et al. (2006)
Inﬁltration of tracers with a tension inﬁltrometer;
inﬁltration rates and tracer concentrations in soil samples
below inﬁltration area
Casey et al. (1998)
Tracer inﬁltration with disk inﬁltrometer for a series ofapplied tensions; inﬁltration rates and tracer concentrations
in soil samples.
Measurements of wet-end K–θ–h relations and results from Logsdon (2002)
image analysis for Mollisols (sandy loam to silty clay).
(continued on next page)

28

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Table 4 (continued)(continued)
Scale Parameters (model)

Determination approach

Authors

Scale IIS
Pore spectrum parameters of a model of
ﬂow in equivalent capillary tubes

Kung et al. (2006)

S

Kasteel et al. (2005)

S

P
F
(UC)
F-UC
F-UC

Tracer mass ﬂux breakthrough patterns under transient
unsaturated ﬂow conditions were used to quantify
equivalent pore spectra of macropore-type preferential ﬂow
pathways.
3D distribution and correlation structure
Brilliant blue dye (BB) and bromide (Br) inﬁltration, imageof BB concentration
analysis (BB), soil sampling (Br).
Statistical categorization of soil layers with cluster analysis is used to derive different transport
ﬂow regimes ranging from homogeneous mechanisms from multiple dye staining patterns in verticalsoil proﬁles from 25 plots
to PF
Ks, λc, θim/θ, α (MIM)
Dripper-TDR ﬁeld method measuring point source chemical
movement
Α (MIM), k (two-site model)
empirical relation between α and residence time, andbetween k and pore water velocity
Correlation between α and residence time Comparison of results from 316 solute transport(and v)
experiments in 35 articles, plus own experiments.
Review of multiple leaching experiments in a range of soils,
Relation between λL, v, K(θ) and soil
properties controlling lateral mixing
linking λL to soil morphological features, and analysis ofdepth dependence of λLλL.
regimeKulli et al. (2003)

Al-Jabri et al. (2002, 2006)
Maraqa (2001)
Haggerty et al., 2004
(Vanderborght et al., 1997; Vanderborght et al.,
2002; Vanderborght and Vereecken, 2007)

Scale I: soil aggregate (Ped), artiﬁcial or repacked soil column (AC), undisturbed column (UC), lysimeter (L); scale II: soil proﬁle (S), plot (P), ﬁeld (F); acronyms:
computer tomography (CT), saturated hydraulic conductivity (Ks), macroscopic capillary length (lambda(c)), immobile water fraction (θim/θ), mass exchange
coefﬁcient (α) and dispersion coefﬁcient (Dm) or dispersivity (λ), saturated and boundary hydraulic conductivities (Ks and Kb), boundary head between macropore
and micropore domains (hb), exponent (n⁎) of relation between K and water content (θ), macropore fraction (θsma), and half spacing (d) between equivalent parallel
fractures, dispersivity (Dv), ρ bulk density, K(θ) hydraulic conductivity function.
adensity, position, length, volume, wall area, hydraulic radius, tortuosity, inclination, number of branches.

Moreover, it was suggested to direct more effort to derivepedotransfer functions accounting for soil structure effects(Merdun and Quisenberry, 2004; Pachepsky et al., 2006). The
relation between lateral dispersivity, unsaturated hydraulicconductivity, and soil morphological properties, was studied
in tracer experiments carried out in columns and ﬁeld plots
with key Belgian soil types (Vanderborght et al., 2001).
Procedures based on combined analyses of soil hydraulicmeasurements, tracer BTCs, and dye stained soil fractions havebeen proposed to relate MIM parameters to soil clay content(Shaw et al., 2000), or to soil structure formation (Vervoort
et al., 1999). The latter authors conceptually linked soil
structure development to the predominant state of saturationof a soil horizon as depending on its location in pedon andlandscape position. Jarvis et al. (2007) analyzed the variationof the mass transfer coefﬁcient, d, in the DPM MACRO (Larsbo
et al., 2005) as a function of the geometric mean particle size
and the organic matter content. First, SUFI (Abbaspour et al.,1997) was utilized to identify d for chloride BTCs observed for
33 intact topsoil columns sampled from 3 different ﬁeld sites.
Multiple linear regression suggested that approximately 50–
60% of the variation in d can be explained by two easily
measured fundamental soil properties (Jarvis et al., 2007).
6.1.5. Summary on independent parameter determination, scale
I (soil column, lysimeter)
Techniques using CT, MRI and tracers were applied either
independently or simultaneously, and their results were
implemented in 3D models to predict PF and transport asfunctions of the soil structure. Additionally, laboratory procedures were tested to determine two-region model parameters,
involving region-related water contents, hydraulic functions,dispersivities, and parameters related to water and solutetransfer. The attempt was also made to statistically relate suchmodel parameters to basic soil properties. More research on thelinks between soil structure and solute transport is needed tosupport the use of DPM for management purposes. Despite theprogress that has been made, the value of these measurement
techniques for ﬁeld conditions has yet to be proven.
6.2. Scale II (plot, ﬁeld)
6.2.1. Field techniques to estimate two-region model parameters
Jaynes et al. (1995) and Clothier et al. (1995) suggested ﬁeld
methods for estimating MIM parameters. Using disk inﬁltrometers, several tracers were inﬁltrated into the soil and theirdistributions below the disk were then measured in soil
samples taken after inﬁltration. A single applied tension yielded
one set of values for the hydraulic conductivity, K, the total
water content, θ, the water content of the immobile zone, θim,and the mass transfer coefﬁcient, α (Alletto et al., 2006), while
several tracer inﬁltrations conducted at different tensions
produced K, θ, θim, and α, as near-saturated functions of thematric potential (Casey et al., 1998). Al-Jabri et al. (2002)
suggested a variant of this method for simultaneous measurement at multiple locations based on point-source dripper lines.
They used the Wooding's (1968) analytical solution for
estimating K and the capillary length, λc, while MIM transportparameters were estimated according to Jaynes et al. (1995). In
an application with several dripping rates, the saturated
hydraulic conductivity, Ks, the ratio of the immobile and thesaturated water content, θim/θs, and λc were similar, while α
was different compared to corresponding values derived using
the disk inﬁltrometer (Al-Jabri et al., 2002). The dripper line
method was further extended by the installation of a TDR probebeneath each inﬁltration spot to also measure solute breakthrough to yield the dispersion coefﬁcient in the mobile region(Al-Jabri et al., 2006). The accuracy of these methods depends
on the validity of the underlying assumptions, including piston
movement, a constant mobile water concentration equal to theJ.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

inﬁltration concentration, and negligible solute exchange with
the immobile region (e.g., Alletto et al., 2006).
Different lab and ﬁeld methods, including disk inﬁltrometer techniques, image analysis, and lab water desorption
tests, were suggested for independent estimation of hydraulicparameters for the MACRO model (Logsdon, 2002).
Kung et al. (2006) analyzed the ﬁeld-scale ﬂux breakthrough patterns of four conservative tracers sequentiallyapplied during a 10-h inﬁltration. They used a six-parameter
capillary tube model to quantify equivalent macropore
spectra of PF pathways. Results indicated that pathways
with increasingly larger equivalent pore radii became
hydraulically active during soil wetting in the course of
inﬁltration (Kung et al., 2006).
6.2.2. Dye tracer experiments
Dye tracer experiments were often used to infer solutetransport processes, test model hypotheses, and estimate
model parameters. Absolute dye concentrations can now be
calculated with reasonable accuracy (Forrer et al., 2000). The
spatial resolution (1 mm) of dye concentration maps ofhorizontal soil cross-sections can be interpolated to infer the
3D distribution of the absolute dye concentrations, sufﬁcientfor characterizing soil structure effects on transport (Kasteel
et al., 2005). An important step towards more realistictransport simulations will be to better understand anddescribe ﬂow at interfaces or transitions between different
soil layers. Kulli et al. (2003) analyzed dye staining patterns invertical soil proﬁles of 25 plots at 8 sites. Zones ofconcentrating ﬂow (attractor zone), spreading ﬂow (dispersion zone), and no alteration of ﬂow regime (transmissionzone) were identiﬁed. A cluster analysis yielded classiﬁcationof seven different transport regimes ranging from homogeneous to preferential. Soils with a similar sequence of soil
horizons did not always display a similar sequence of certaintransport features. Other factors, such as gradual or sharp
textural layer boundaries, may invoke different ﬂow patterns
in similar underlying horizons (Kulli et al., 2003).
6.2.3. Combinations of techniques: CT, dye tracing, andbreakthrough curves
Hydrogeophysical, dye tracing, and breakthrough analyses
alone have their conceptual limitations. Using CT allows
accurate macropore network reconstruction, but may notshow which part of the network will be hydraulically active. A
dye tracing proﬁle reveals the macropore network engaged indye transport, but adsorption may restrict the dye from partof the water ﬂow network. A leaching BTC provides atemporal pattern of the system response, but does not allow
conclusions about the internal spatial transport dynamics.
However, complementary information on transport processes
may be gained when using multiple techniques simultaneously (e.g., Vanderborght et al., 2002). Hydrogeophysicalﬁeld monitoring techniques (e.g. Gaur et al., 2006) should befurther developed for integrating both spatial and temporal
ﬂow-transport observations. The processes linking ‘form’ (soilstructure, layering, interfaces) and ‘functions’ (solute transport patterns) should be further explored utilizing tracer
experiments (Kulli et al., 2003; Kasteel et al., 2005) for
structure based prediction of solute transport (Vervoort et al.,
1999; Vogel et al., 2006).

29

7. Conclusions
7.1. Water ﬂow
When only simulating water ﬂow (without solute transport)
through structured soils, equilibrium single-domain models
often yield results similar to two-domain (MIM, DPM)
approaches, unless dynamic shrinkage cracks are present. An
inverse parameter identiﬁcation of two-domain PF parameters
based on water outﬂow information suffers from non-uniqueness problems. At the ﬁeld scale, stochastic methods were
utilized to analyze the effects of heterogeneity, parameter
uncertainty, and PF triggering processes, such as precipitation
rate and depth thresholds. The effects of domain-separated
boundary conditions in two ﬂow-domain approaches need tobe considered. Complex catchment area hydrology wasdescribed with distributed hydrological unit models or 3D
coupled numerical bypass ﬂow models, some of which
considered lateral PF as important hillslope ﬂow process.
Threshold values of rainfall intensities apparently also trigger
PF at the catchment area scale. However, scarce data and model
dicretization beyond the process scale make ﬂow model
validation at the catchment area scale complicated.
7.2. Solute transport
Model approaches applied at the ﬁeld, and smaller, scales
often involved SPM, MIM, and DPM or three-domain models
that can consider natural transient ﬂow conditions. Alternative
approaches included DFM, 1D stream tube models (STM),
method of moments, a stochastic multifractal approach, adrainage network reconstruction, a structure-based 3D ﬂow
transport model and hillslope models with lateral pipe ﬂow.
Several inverse, uncertainty, and parameter sensitivity quantiﬁcation concepts were tested. Inverse schemes need to makebetter use of localized important information or patterns.
Parameter uncertainty is a useful tool, e.g., for including effectsof spatial variability in ﬁeld-scale predictions. A systematic
correlation between parameter sensitivities and soil properties(for experiments conducted at a particular scale) might further
help to identify important parameters for certain conditions.
The effects of experimental conditions on PNE solute transport
were investigated at the column scale. 2D, inverse, andtransient two-region models were applied at the ﬁeld scale.
At the catchment area scale, evaluating PF effects on solutetransport has only just begun and will require innovative
approaches, such as combining tracer tests and modelling at
nested scales (plot-ﬁeld-catchment area), geophysical andremote sensing techniques, stochastic methods and scenario
analysis (‘virtual experiments’).
7.3. Parameter determination
Progress includes both, independent determination ofparameters of two-region models, and the development ofother predictive model approaches. At the column scale, thecombined use of CT, MRI and tracer techniques showedpromise in predicting PF and transport as a function of the soil
structure. Pedotransfer functions that were statisticallyrelated to fundamental soil properties, as proxies of soil
structure, may permit a rough prediction of two-region model

30

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

parameters. Despite the progress made, an independent
determination of hydraulic and transport parameters is still
very limited because of the simplifying assumptions underlying various methods, and often requires separation ofanalyzed soil samples from their natural setting. Methods
that showed potential for parameter estimation and predictivemodel capabilities at the ﬁeld scale included (i) experimental
methods, e.g., a combination of the tension inﬁltration methodwith tracer applications, an analysis of ﬂux breakthroughpatterns of sequentially applied tracers, or combined dye
tracer and hydraulic measurements for setting up a 3D soil
hydraulic structure model, and (ii) statistical evaluations, e.g.,of relations between solute transfer coefﬁcients and theresidence time, or classiﬁcation of transport regimes using
cluster analysis of dye concentration spatial patterns indifferent sequences of soil horizons. Soil structure and its
spatial variability need to be jointly considered in modelling at
the ﬁeld and catchment area scales relevant for agriculturaland environmental management.
Acknowledgements
The ﬁrst author was supported by German Research Foundation (DFG) under project KO 2899/1-1. The authors thank Nick
Jarvis for his detailed and constructive review of a draft version ofthis paper. The helpful comments by the editor-in-chief, Emil
Frind, the managing guest editors for this special issue ‘Flow
domains’, Miroslav Kutilek and Antonio Coppola, and an
anonymous reviewer are gratefully acknowledged. Steffen
Schlüter provided technical assistance in preparing Fig. 2.
References
Abbaspour, K.C., van Genuchten, M.Th., Schulin, R., Schläppi, E., 1997. A
sequential uncertainty domain inverse procedure for estimating subsurface ﬂow and transport parameters. Water Resour. Res. 33, 1879–1892.
Abbaspour, K.C., Kohler, A., Šimůnek, J., Fritsch, M., Schulin, R., 2001.
Application of a two-dimensional model to simulate ﬂow and transport
in a macroporous agricultural soil with tile drains. Europ. J. Soil Sci. 52 (3),433–447.
Abbasi, F., Šimůnek, J., Feyen, J., van Genuchten, M.Th., Shouse, P.J., 2003.
Simultaneous inverse estimation of soil hydraulic and solute transportparameters from transient ﬁeld experiments: homogeneous soil. Trans.
ASAE 46 (4), 1085–1095.
Abbasi, F., Feyen, J., van Genuchten, M.Th., 2004. Two-dimensional simulationof water ﬂow and solute transport below furrows: model calibration andvalidation. J. Hydrol. 290, 63–79.
Abrahamsen, P., Hansen, S., 2000. Daisy: an open soil–crop–atmospheresystem model. Environ. Model. Software 15, 313–330.
Aden, K., Diekkrüger, B., 2000. Modelling pesticide dynamics of four different
sites using the model system SIMULAT. Agric. Water Manag. 44 (1), 337–355.
Addiscott, T.M., Heye, P.J., Whitmore, A.P., 1986. Application of simple
leaching models in heterogeneous soils. Geoderma 38, 185–194.
Addiscott, T.M., Whitmore, A.P., 1991. Simulation of solute leaching in soils ofdiffering permeabilities. Soil Use Manag. 7, 94–102.
Ahuja, L.R., Johnson, K.E., Heathman, G.C., 1995. Macropore transport of asurface-applied bromide tracer: Model evaluation and reﬁnement. Soil
Sci. Soc. Am. J. 59, 1234–1241.
Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, J.J., Ma, L. (Eds.), 2000. The Root Zone
Water Quality Model, Water Resources Publications LLC, Highlands Ranch, CO.
Akhand, N.A., Lapen, D.R., Topp, E., Edwards, M.J., Sabourin, L., Ball Coelho, B.R.,
Duenk, P.W., Payne, M., 2006. Prediction of liquid municipal biosolid andprecipitation induced tile ﬂow in a Southern Ontario agricultural ﬁeld
using MACRO. Agricult. Water Manag. 83 (1–2), 37–50.
Al-Huthali, A., Datta-Gupta, A., 2004. Streamline simulation of countercurrent imbibition in naturally fractured reservoirs. J. Petrol. Sci.
Engineer. 43 (3–4), 271–300.
Al-Jabri, S.A., Horton, R., Jaynes, D.B., 2002. A point-source method for rapid
simultaneous estimation of soil hydraulic and chemical transport
properties. Soil Sci. Soc. Am. J. 66, 12–18.

Al-Jabri, S.A., Lee, J., Gaur, A., Horton, R., Jaynes, D.B., 2006. A dripper-TDR
method for in situ determination of hydraulic conductivity and chemicaltransport properties of surface soils. Adv. Water Resour. 29 (2), 239–249.
Alaoui, A., Germann, P., Jarvis, N., Acutis, M., 2003. Dual-porosity andkinematic wave approaches to assess the degree of preferential ﬂow inan unsaturated soil. Hydrol. Sci. J. — Journal des Sciences Hydrologiques
48 (3), 455–472.
Allaire, S.E., Gupta, S.C., Nieber, J., Moncrief, J.F., 2002a. Role of macropore
continuity and tortuosity on solute transport in soils: 1. Effects of initialand boundary conditions. J. Contam. Hydrol. 58, 299–321.
Allaire, S.E., Gupta, S.C., Nieber, J., Moncrief, J.F., 2002b. Role of macropore
continuity and tortuosity on solute transport in soils: 2. Interactions withmodel assumptions for macropore description. J. Contam. Hydrol. 58 (3–4),
283–298.
Alletto, L., Coquet, Y., Vachier, P., Labat, C., 2006. Hydraulic conductivity,
immobile water content, and exchange coefﬁcient in three soil proﬁles.
Soil Sci. Soc. Am. J. 70, 1272–1280.
Anderson, A.N., Crawford, J.W., McBratney, A.B., 2000. On diffusion in fractal
soil structures. Soil Sci. Soc. Am. J. 64, 19–24.
Armstrong, A.C., Matthews, A.M., Portwood, A.M., Leeds-Harrison, P.B., Jarvis,
N.J., 2000. CRACK-NP: a pesticide leaching model for cracking clay soils.
Agricult. Water Manag. 44, 183–199.
Barenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theoryof seepage of homogeneous liquids in ﬁssured rocks. J. Appl. Math. Mech.
24, 1286–1303.
Baumann, Th., Petsch, R., Niessner, R., 2000. Direct 3-D measurement of theﬂow velocity in porous media using magnetic resonance tomography.
Environ. Sci. Technol. 34, 4242–4248.
Baveye, P., Boast, C.W., Ogawa, S., Parlange, J.-Y., Steenhuis, T., 1998. Inﬂuenceof image resolution and thresholding on the apparent mass fractal
characteristics of preferential ﬂow patterns in ﬁeld soils. Water Resour.
Res. 34, 2783–2796.
Beckers, J., Alila, Y., 2004. A model of rapid preferential hillslope runoff
contributions to peak ﬂow generation in a temperate rain forest
watershed. Water Resour. Res. 40 (3) Art. No. W03501.
Beven, K., Germann, P., 1982. Macropores and water ﬂow in soils. Water
Resour. Res. 18, 1311–1325.
Beven, K., Binley, A., 1992. The future of distributed models: model calibration
and uncertainty prediction. Hydrol. Process 6 (3), 279–298.
Beven, K., Freer, J., 2001. Equiﬁnality, data assimilation, and uncertainty
estimation in mechanistic modelling of complex environmental systems
using the GLUE methodology. J. Hydrol. 249, 11–29.
Birk, S., Liedl, R., Sauter, M., 2006. Karst spring responses examined byprocess-based modelling. Ground Water 44 (6), 832–836.
Booltink, H.W.G., 1994. Field-scale distributed modelling of bypass ﬂow in aheavily textured clay soil. J. Hydrol. 163 (1–2), 65–84.
Bormann, H., Fass, T., Giertz, S., Junge, B., Diekkrüger, B., Reichert, B.,
Skowronek, A., 2005. From local hydrological process analysis to regionalhydrological model application in Benin: concept, results and perspectives. Phys. Chem. Earth 30 (6–7), 347–356.
Boufadel, M.C., Lu, S., Molz, F.J., Lavallée, D., 2000. Multifractal scaling of the
intrinsic permeability. Water Resour. Res. 36, 3211–3222.
Bouma, J., 1990. Using morphometric expressions for macropores to improvesoil physical analyses of ﬁeld soils. Geoderma 46, 3–11.
Brauchler, R., Liedl, R., Dietrich, P., 2003. A travel time based hydraulictomographic approach. Water Resour. Res. 39 (12), 1370. doi:10.1029/
2003WR002262.
Bronstert, A., 1994. Modellierung der Abﬂussbildung und der Bodenwasserdynamik von Hängen. Universität Karlsruhe, Inst. f. Hydrologie u.
Wasserwirtschaft, Dissertation Thesis, No. 46.
Buttle, J.M., McDonald, D.J., 2002. Coupled vertical and lateral preferential
ﬂow on a forested slope. Water Resour. Res. 38 (5), 1060. doi:10.1029/
2001WR000773.
Casey, F.X.M., Logsdon, S.D., Horton, R., Jaynes, D.B., 1998. Measurement of ﬁeld soil
hydraulic and solute transport parameters. Soil Sci. Soc. Am. J. 62, 1172–1178.
Castiglione, P., Mohanty, B.P., Shouse, P.J., Šimůnek, J., van Genuchten, M.Th.,
Santini, A., 2003. Lateral water diffusion in an artiﬁcial macroporous
system: modelling and experimental evidence. Vadose Zone J. 2, 212–221.
Cey, E., Rudolph, D., Therrien, R., 2006. Simulation of groundwater recharge
dynamics in partially saturated fractured soils incorporating spatiallyvariable fracture apertures. Water Resour. Res. 42, W09413. doi:10.1029/
2005WR004589.
Christiansen, J.S., Thorsen, M., Clausen, T., Hansen, S., Refsgaard, J.C., 2004.
Modelling of macropore ﬂow and transport processes at catchment scale.
J. Hydrol. 299 (1–2), 136–158.
Císlerová, M., Votrubová, J., 2002. CT derived porosity distribution and ﬂow
domains. J. Hydrol. 267, 186–200.
Clothier, B.E., Heng, L., Magesan, G.N., Vogeler, I., 1995. The measured mobilewater content of an unsaturated soil as a function of hydraulic regime.
Aust. J.Soil Res. 33, 397–414.

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35
Clothier, B.E., Green, S.R., Deurer, M., 2008. Preferential ﬂow and transport insoil: progress and prognosis. Europ. J. Soil Sci. 59, 2–13.
Coats, K.H., Smith, B.D., 1964. Dead-end pore volume and dispersion in
porous media. Soc. Petrol. Eng. J. 4, 73–84.
Comegna, V., Coppola, A., Sommella, A., 2001. Effectiveness of equilibrium
and physical non-equilibrium approaches for interpreting solute transport through undisturbed soil columns. J. Contam. Hydrol. 50 (1–2),
121–138.
Coppola, A., 2000. Unimodal and bimodal description of hydraulic properties
for aggregated soils. Soil Sci. Soc. Am. J. 64, 1252–1262.
Cote, C.M., Bristow, K.L., Ross, P.J., 1999. Quantifying the inﬂuence of intraaggregate concentration gradients on solute transport. Soil Sci. Soc. Am. J.
63 (4), 759–767.
Croton, J.T., Barry, D.A., 2001. WEC-C: a distributed, deterministic catchment
model — theory, formulation and testing. Environ. Model. Software 16,
583–599.
Dasgupta, S., Mohanty, B., Köhne, M., 2006a. Impacts of juniper vegetationand karst geology on subsurface ﬂow processes in the Edwards Plateau,
Texas. Vadose Zone J. 5, 1076–1085.
Dasgupta, S., Mohanty, B., Köhne, M., 2006b. Soil hydraulic conductivities andtheir spatial and temporal variations in a vertisol, Soil Sci. Soc. Am. J. 70,
1872–1881.
De Rooij, G.H., 2000. Modelling ﬁngered ﬂow of water in soils owing towetting front instability: a review. J. Hydrol. 231, 277–294.
Deurer, M., Green, S.R., Clothier, B.E., Böttcher, J., Duijnisveld, W.H.M., 2003.
Drainage networks in soils. A concept to describe bypass-ﬂow pathways.
J. Hydrol. 272 (1–4), 148–162.
Diekkrüger, B., 1996. SIMULAT - Ein Modellsystem zur Berechnung derWasser- und Stoffdynamik landwirtschaftlich genutzter Standorte. In:
Richter, O., Söndgerath, D., Diekkrüger, B. (Eds.), Sonderforschungsbereich 179 `Wasser- und Stoffdynamik in Agrarökosystemen'. Selbstverlag,
Institut für Geographie und Geoökologie der Technischen Universität
Braunschweig, Braunschweig, pp. 30–47.
Di Pietro, L., Ruya, S., Capowiez, Y., 2003. Predicting preferential water ﬂow insoils by traveling-dispersive waves. J. Hydrol. 278 (1–4), 64–75.
Dragila, M.-I., Weisbrod, N., 2004. Flow in menisci corners of capillary
rivulets. Vadose Zone J. 3, 1439–1442.
Duguid, J.O., Lee, P.C.Y., 1977. Flow in fractured porous media. Water Resour.
Res. 13, 558–566.
Dunn, S.M., Freer, J., Weiler, M., Kirkby, M.J., Seibert, J., Quinn, P.F., Lischeid, G.,
Tetzlaff, D., Soulsby, C., 2008. Conceptualizing catchment modelling: simply
learning? Hydrol. Process. 22, 2389–2393.
Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water Resour. Res. 30 (2), 211–223.
Eckersten, H., Jansson, P.-E., Johnsson, H., 1998. SOILN Model, Version 9.2,
User's Manual. Division of Hydrotechnics, Communication 98:6. Swedish
Agricultural University, Uppsala.
Ersahin, S., Papendick, R.I., Smith, J.L., Keller, C.K., Manoranjan, V.S., 2002.
Macropore transport of bromide as inﬂuenced by soil structure
differences. Geoderma 108 (3–4), 207–223.
Famiglietti, J.S., Wood, E.F., 1994. Multi-scale modelling of spatially variable
water and energy balance processes. Water Resour. Res. 30 (11),3061–3078.
Feyen, J., Jacques, D., Timmerman, A., Vanderborght, J., 1998. Modelling waterﬂow and solute transport in heterogeneous soils: a review of recent
approaches. J. Agric. Eng. Res. 70, 231–256.
Flühler, H., Durner, W., Flury, M., 1996. Lateral solute mixing processes — a
key for understanding and modelling ﬁeld-scale transport of water andsolutes. Geoderma 70 (2–4), 165–183.
Forrer, I., Papritz, A., Kasteel, R., Flühler, H., Luca, D., 2000. Quantifying dye
tracers in soil proﬁles by image processing. Europ. J. Soil Sci. 51, 313–322.
Gaudet, J.P., Jégat, H., Vachaud, G., Wierenga, P.J., 1977. Solute transfer, withexchange between mobile and stagnant water, through unsaturated
sand. Soil Sci. Soc. Am. J. 41, 665–671.
Gaur, A., Horton, R., Jaynes, D.B., Baker, J.L., 2006. Measured and predicted
solute transport in a tile drained ﬁeld. Soil Sci. Soc. Am. J. 70 (3), 872–881.
Gérard, F., Tinsley, M., Mayer, K.U., 2004. Preferential ﬂow revealed byhydrologic modelling based on predicted hydraulic properties. Soil Sci.
Soc. Am. J. 68 (5), 1526–1538.
Gerke, H.H., 2006. Preferential ﬂow descriptions for structured soils. J. Plant
Nutr. Soil Sci. — Z. Pﬂanzenernähr. Bodenk. 169 (3), 382–400.
Gerke, H.H., van Genuchten, M.Th., 1993. A dual-porosity model for
simulating the preferential movement of water and solutes in structured
porous media. Water Resour. Res. 29, 305–319.
Gerke, H.H., van Genuchten, M.T., 1996. Macroscopic representation of
structural geometry for simulating water and solute mass transfer indual-porosity media. Adv. Water Resour. 19, 343–357.
Gerke, H.H., Köhne, J.M., 2002. Estimating hydraulic properties of soil
aggregate skins from sorptivity and water retention. Soil Sci. Soc. Am. J.
66 (1), 26–36.

31

Gerke, H.H., Köhne, J.M., 2004. Dual-permeability modelling of preferential bromide leaching from a tile drained glacial till agriculturalﬁeld, J. Hydrol. 289, 239–257.
Gerke, H.H., Dusek, J., Vogel, T., Köhne, J.M., 2007. 2D dual-permeability
analyses of a bromide tracer experiment on a tile-drained ﬁeld. Vadose
Zone J. 6, 651–667.
Germann, P.F., 1985. Kinematic wave approximation to inﬁltration anddrainage into and from soil macropores. Transactions ASAE 28, 745–749.
Germann, P.F., Beven, K., 1981. Water ﬂow in soil macropores. III. A statistical
approach. J. Soil Sci. 32, 31–39.
Germann, P.F., Di Pietro, L., 1996. When is porous-media ﬂow preferential? Ahydromechanical perspective. Geoderma 74, 1–21.
Germann, P.F., Di Pietro, L., 1999. Scales and dimensions of momentum dissipation
during preferential ﬂow in soils. Water Resour. Res. 35 (5), 1443–1454.
Ghezzehei, T.A., 2004. Constraints for ﬂow regimes on smooth fracturesurfaces. Water Resour. Res. 40. doi:10.1029/2004WR003164.
Ghodrati, M., Chendorain, M., Jason Chang, Y., 1999. Characterization ofmacropore ﬂow mechanisms in soil by means of a split macropore
column. Soil Sci. Soc. Am. J. 63, 1093–1101.
Graf and Therrien 2006.
Grathwohl, P., Kraft, S., Lamprou, A., 2004. Beeinﬂussung von in-situ Reaktoren
durch grundwasserspeziﬁsche Parameter am Beispiel der in-situ Aktivkohleﬁltration. In: Weiß, H., Teutsch, G., Daus, B. (Eds.), UFZ-Bericht 13/
2004. Helmholtz Centre For Environmental Research, Leipzig, pp. 19–40.
Greco, R., 2002. Preferential ﬂow in macroporous swelling soil with internal
catchment: model development and applications. J. Hydrol. 269 (3–4),
150–168.
Grifﬁoen, J.W., 1998. Suitability of the ﬁrst-order mass transfer concept for
describing cyclic diffusive mass transfer in stagnant zones. J. Contam.
Hydrol. 34, 155–165.
Grifﬁoen, J.W., Barry, D.A., Parlange, J.Y., 1998. Interpretation of two-region
model parameters. Water Resour. Res. 34 (3), 373–384.
Gwo, J.-P., 2001. In search of preferential ﬂow paths in structured porous
media using a simple genetic algorithm. Water Resour. Res. 37, 1589–1601.
Gwo, J.P., Jardine, P.M., Wilson, G.V., Yeh, G.T., 1995. multiple-pore-region
concept to modeling mass transfer in subsurface media, J. Hydrol. 164,
217–237.
Gwo, J.P., Jardine, P.M., Wilson, G.V., Yeh, G.T., 1996. Using a multiregion
model to study the effects of advective and diffusive mass transfer on
local physical nonequilibrium and solute mobility in a structured soil.
Water Resour. Res. 32, 561–570.
Gwo, J.P., O'Brien, R., Jardine, P.M., 1998. Mass transfer in structured porous
media: embedding mesoscale structure and microscale hydrodynamics
in a two-region model. J. Hydrol. 208 (3–4), 204–222.
Haggerty, R., Gorelick, S.M., 1995. Multiple-rate mass transfer for modeling
diffusion and surface reactions in media with pore-scale heterogeneity.
Water Resour. Res. 31 (10), 2383–2400.
Haggerty, R., Gorelick, S.M., 1998. Modeling mass transfer processes in soilcolumns with pore-scale heterogeneity. Soil Sci. Soc. Am. J. 62 (1), 62–74.
Haggerty, R., Harvey, C.F., von Schwerin, C.F., Meigs, L.C., 2004. What controls theapparent timescale of solute mass transfer in aquifers and soils? A comparisonof experimental results. Water Resour. Res. 40 (1) Art. No. W01510.
Hansen, S., Jensen, H.E., Nielsen, N.E., Svendsen, H., 1991. Simulation of
nitrogen dynamics and biomass production in winter wheat using theDanish simulation model Daisy. Fertil. Res. 27, 245–259.
Hatano, R., Booltink, H.W.G., 1992. Using fractal dimensions of stained ﬂow
patterns in a clay soil to predict bypass ﬂow. J. Hydrol. 135, 121–131.
Haws, N.W., Rao, P.S.C., 2004. The effect of vertically decreasing macropore
fractions on simulations of non-equilibrium solute transport. Vadose
Zone J. 3 (4), 1300–1308.
Haws, N.W., Rao, P.S.C., Šimůnek, J., Poyer, I.C., 2005. Single-porosity and dualporosity modelling of water ﬂow and solute transport in subsurface-drained
ﬁelds using effective ﬁeld-scale parameters. J. Hydrol. 313 (3–4), 257–273.
Haws, N.W., Das, B.S., Rao, P.S.C., 2004. Dual-domain solute transfer andtransport processes: evaluation in batch and transport experiments.
J. Contam. Hydrol. 75, 257–280.
Hutson, J.L., Wagenet, R.J., 1992. LEACHM: leaching estimation and chemistrymodel version 3. Cornell Univ. Dept. of Soil, Crop and Atmospheric
Sciences, Research Series No. 92-3.
Hutson, J.L., Wagenet, R.J., 1995. A multiregion model describing water ﬂow andsolute transport in heterogeneous soils. Soil Sci. Soc. Am. J. 59, 743–751.
Ilsemann, J., van der Ploeg, R.R., Horton, R., Bachmann, J., 2002. Laboratory
method for determining immobile soil water content and mass exchange
coefﬁcient. J. Plant Nutr. Soil Sci. 165, 332–338.
Jabro, J.D., Toth, J.D., Fox, R.H., 1998. Evaluation and comparison of ﬁve
simulation models for estimating water drainage ﬂuxes under corn.
J. Environ. Qual. 27, 1376–1381.
Jansson, P.-E., 1991. Simulation Model for Soil Water and Heat Conditions.
Description of SOIL Model. Swedish University of Agricultural Sciences,
Department of Soil Sciences, Uppsala, Sweden. Report No. 165, 73 pp.

32

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Jansson, P.-E., Karlberg, L., 2004. Coupled Heat and Mass Transfer Model for
Soil–Plant–Atmosphere Systems. Royal Institute of Technology, Department of Civil and Environmental Engineering, Stockholm. 435 pp.
Jarvis, N.J., 1994. The MACRO model (Version 3.1). Technical Description andSample Simulations. Reports and Dissert, 19. Dept. Soil Sci., Swedish Univ.
Agric. Sci., Uppsala, Sweden.
Jarvis, N.J., 1998. Modelling the impact of preferential ﬂow on non-point
source pollution. In: Selim, H.H., Ma, L. (Eds.), Physical Non-Equilibriumin Soils: Modelling and Application. Ann Arbor Press, pp. 195–221.
Jarvis, N.J., 2007. Review of non-equilibrium water ﬂow and solute transport
in soil macropores: principles, controlling factors and consequences forwater quality. Eur. J. Soil Sci. 58, 523–546.
Jarvis, N., Larsbo, M., Roulier, S., Lindahl, A., Persson, L., 2007. The role of soil
properties in regulating non-equilibrium macropore ﬂow and solutetransport in agricultural topsoils. Eur. J. Soil Sci. 58, 282–292.
Jaynes, D.B., Logsdon, S.D., Horton, R., 1995. Field method for measuring mobile/immobile water content and solute transfer rate. Soil Sci. Soc. Am. J. 59,
352–356.
Javaux, M., Kasteel, R., Vanderborght, J., Vanclooster, M., 2006. Interpretationof dye transport in a macroscopically heterogeneous, unsaturated subsoil
with a one-dimensional model. Vadose Zone Journal 5 (2), 529–538.
Jørgensen, P.R., Helstrup, T., Urup, J., Seifert, D., 2004. Modelling of nonreactive solute transport in fractured clayey till during variable ﬂow rate
and time. J. Contam. Hydrol. 68 (3–4), 193–216.
Jones, J.P., Sudicky, E.A., McLaren, R.G., 2008. Application of a fully-integrated
surface-subsurface ﬂow model at the watershed-scale: a case study.
Water Resour. Res. 44, W03407. doi:10.1029/2006WR005603.
Ju, S.H., Kung, K.J.S., 1997. Impact of funnel ﬂow on contaminant transport insandy soils: Numerical simulation. Soil Sci. Soc. Am. J. 61 (2), 409–415.
Jury, W.A., Roth, K., 1990. Transfer Functions and Solute Movement ThroughSoil: Theory and Application. Birkhäuser, Basel, Switzerland.
Kamra, S.K., Lennartz, B., van Genuchten, M.Th., Widmoser, P., 2001.
Evaluating non-equilibrium solute transport in small soil columns.
J. Contam. Hydrol. 48 (3–4), 189–212.
Kasteel, R., Vogel, H.-J., Roth, K., 2000. From local hydraulic properties toeffective transport in soil. Europ. J. Soil Sci. 51, 81–91.
Kasteel, R., Vogel, H.-J., Roth, K., 2002. Effect of non-linear adsorption on thetransport behaviour of Brilliant Blue in a ﬁeld soil. Eur. J. Soil Sci. 53, 231–240.
Kasteel, R., Burkhardt, M., Giesa, S., Vereecken, H., 2005. Characterization of ﬁeld
tracer transport using high-resolution images. Vadose Zone Journal 4,101–111.
Kätterer, T., Schmied, B., Abbaspour, K.C., Schulin, R., 2001. Single- and dualporosity modelling of multiple tracer transport through soil columns:
effects of initial moisture and mode of application. Europ. J. Soil Sci. 52 (1),25–36.
Kjaergaard, C., Poulsen, T.G., Moldrup, P., de Jonge, L.W., 2004. Colloid
mobilization and transport in undisturbed soil columns. I. Pore structure
characterization and tritium transport. Vadose Zone J. 3 (2), 413–423.
Kodešová, R., Kodeš, V., Žigová, A., Šimůnek, J., 2006. Impact of Plant Rootsand Soil Organisms on Soil Micromorphology and Hydraulic Properties,
61. Biologia, Bratislava, pp. 339–343. /Suppl. 19.
Kohler, A., Abbaspour, K.C., Fritsch, M., van Genuchten, M.Th., Schulin, R.,
2001. Simulating unsaturated ﬂow and transport in a macroporous soil totile drains subject to an entrance head: model development andpreliminary evaluation. J. Hydrol. 254, 67–81.
Kohler, A., Abbaspour, K.C., Fritsch, M., Schulin, R., 2003. Using simple bucket
nodels to analyze solute export to subsurface drains by preferential ﬂow.
Vadose Zone J. 2, 68–75.
Köhne, J.M., Gerke, H.H., Köhne, S., 2002a. Effective diffusion coefﬁcients ofsoil aggregates with surface skins. Soil Sci. Soc. Am. J. 66 (5), 1430–1438.
Köhne, J.M., Köhne, S., Gerke, H.H., 2002b. Estimating the hydraulic functionsof dual-permeability models from bulk soil data. Water Resour. Res.
doi:10.1029/2001 WR000492.
Köhne, J.M., Mohanty, B.P., Šimůnek, J., Gerke, H.H., 2004a. Evaluation of a secondorder water transfer term for variably saturated dual-permeability ﬂow
models. Water Resour. Res. vol. 40, W07409. doi:10.1029/2004WR003285.
Köhne, J.M., Köhne, S., Mohanty, B.P., Šimůnek, J., 2004b. Inverse mobile–immobile modelling of transport during transient ﬂow: effects of
between-domain transfer and initial water content. Vadose Zone J. 3,
1309–1321.
Köhne, J.M., Mohanty, B.P., 2005. Water ﬂow processes in a soil column with acylindrical macropore: experiment and hierarchical modelling. Water
Resour. Res. vol. 41. doi:10.1029/2004WR003303 W03010.
Köhne, J.M., Gerke, H.H., 2005. Spatial and temporal dynamics of preferential
bromide movement towards a tile drain. Vadose Zone J. 4, 79–88.
Köhne, J.M., Mohanty, B.P., Šimůnek, J., 2006a. Inverse dual-permeability
modelling of preferential water ﬂow in a soil column and implications forﬁeld-scale solute transport. Vadose Zone J. 5, 59–76.
Köhne, S., Lennartz, B., Köhne, J.M., Šimůnek, J., 2006b. Bromide transport at a tiledrained ﬁeld site: experiment, and one- and two-dimensional equilibriumand non-equilibrium numerical modelling. J. Hydrol. 321 (1–4), 390–408.

Kramers, G., van Dam, J.C., Ritsema, C.J., Stagnitti, F., Oostindie, K., Dekker, L.W.,
2005. A new modelling approach to simulate preferential ﬂow andtransport in water repellent porous media: Parameter sensitivity, andeffects on crop growth and solute leaching. Austr. J. Soil Res. 43 (3),371–382.
Kulli, B., Stamm, Ch., Papritz, A., Flühler, H., 2003. Discrimination of ﬂow
regions on the basis of stained inﬁltration patterns in soil proﬁles. Vadose
Zone J. 2, 338–348.
Kung, K.-J.S., Kladivko, E.J., Helling, C.S., Gish, T.J., Steenhuis, T.S., Jaynes, D.B.,
2006. Quantifying the pore size spectrum of macropore-type preferential
pathways under transient ﬂow. Vadose Zone J. 5, 978–989.
Langner, H.W., Gaber, H.M., Wraith, J.M., Huwe, B., Inskeep, W.P., 1999.
Preferential ﬂow through intact soil cores: effects of matric head. Soil Sci.
Soc. Am. J. 63 (6), 1591–1598.
Lehmann, P., Hinz, C., McGrath, G., Tromp-van Meerveld, H.J., McDonnel, J.J.,
2007. Rainfall threshold for hillslope outﬂow: an emergent property of
ﬂow pathway connectivity. Hydrol. Earth System Sci. 11, 1047–1063.
Larsbo, M., Jarvis, N.J., 2003. MACRO5.0. A model of water ﬂow and solutetransport in macroporous soil. Technical Description. Emergo 2003, 6.
Swedish Univ. of Agricult. Sci., Dept. of Soil Sci., Uppsala, Sweden. 48 pp.
Larsbo, M., Roulier, S., Stenemo, F., Kasteel, R., Jarvis, N., 2005. An improved
dual-permeability model of water ﬂow and solute transport in thevadose zone. Vadose Zone J. 4, 398–406.
Larsbo, M., Jarvis, N.J., 2005. Simulating solute transport in a structured ﬁeld
soil: uncertainty in parameter identiﬁcation and predictions. J. Environ.
Quality 34, 621–634.
Larsbo, M., Jarvis, N., 2006. Information content of measurements from tracer
microlysimeter experiments designed for parameter identiﬁcation indual-permeability models. J. Hydrol. 325 (1–4), 273–287.
Larsson, M.H., Jarvis, N.J., 1999a. Evaluation of a dual-porosity model topredict ﬁeld-scale solute transport in a macroporous soil. J. Hydrol. 215(1–4), 153–171.
Larsson, M.H., Jarvis, N.J., 1999b. A dual-porosity model to quantify macropore
ﬂow effects on nitrate leaching. J. Environ. Quality 28, 1298–1307.
Lee, J., Horton, R., Jaynes, D.B., 2000. A time domain reﬂectometry method tomeasure immobile water content and mass exchange coefﬁcient. Soil Sci.
Soc. Am. J. 64, 1911–1917.
Lee, J., Horton, R., Noborio, K., Jaynes, D.B., 2001. Characterization ofpreferential ﬂow in undisturbed, structured soil columns using a vertical
TDR probe. J. Contam. Hydrol. 51, 131–144.
Lee, J., Horton, R., Jaynes, D.B., 2002. The feasibility of shallow time domain
reﬂectometry probes to describe solute transport through undisturbed
soil columns. Soil Sci. Soc. Am. J. 66, 53–57.
Lewandowska, J., Szymkiewicz, A., Burzyński, K., Vauclin, M., 2004. Modellingof unsaturated water ﬂow in double porosity soils by the homogenization approach. Advances in Water Resources 27, 283–296.
Lewandowska, J., Szymkiewicz, A., Gorczewska, W., Vauclin, M., 2005.
Inﬁltration in a double-porosity medium: experiments and comparison
with a theoretical model. Water Resour. Res. 41 (2) Art. No. W02022.
Li, M.H., Cheng, H.P., Yeh, G.T., 2005. An adaptive multigrid approach for thesimulation of contaminant transport in the 3D subsurface. Comput.
Geosci. 31 (8), 1028–1041.
Li, Y., Ghodrati, M., 1995. Transport of nitrate in soils as affected by earthworm
activities. J. Environ. Qual. 24, 432–438.
Li, Y., Ghodrati, M., 1997. Preferential transport of solute through constructed
macropores. Soil Sci. Soc. Am. J. 61, 1308–1317.
Lin, H.S., Kogelmann, W., Walker, C., Bruns, M.A., 2006. Soil moisture patterns
in a forested catchment: a hydropedological perspective. Geoderma 131,
345–368.
Lin, H., Zhou, X., 2008. Evidence of subsurface preferential ﬂow using soilhydrologic monitoring in the Shale Hills catchment. Eur. J. Soil Sci. 59, 34–49.
Liu, H.H., Zhang, R., Bodvarsson, G.S., 2005. An active region model for
capturing fractal ﬂow patterns in unsaturated soils: model development.
J. Contam. Hydrol. 80 (1–2), 18–30.
Logsdon, S.D., 1995. Flow mechanisms through continuous and buried
macropores. Soil Sci. 160, 237–242.
Logsdon, S.D., Keller, K.E., Moorman, T.B., 2002. Measured and predicted
solute leaching from multiple undisturbed soil columns. Soil Sci. Soc. Am.
J. 66 (3), 686–695.
Logsdon, S.D., 2002. Determination of preferential ﬂow model parameters.
Soil Sci. Soc. Am. J. 66, 1095–1103.
Ludwig, R., Gerke, H.H., Wendroth, O., 1999. Describing water ﬂow inmacroporous ﬁeld soils using the modiﬁed macro model. J. Hydrol. 215 (1),
135–152 1999.
Ma, L., Ahuja, L.R., Malone, R.W., 2007. Systems modelling for soil and water
research and management: current status and needs for the 21st century.
Trans. ASABE 50 (5), 1705–1713.
Malone, R.W., Shipitalo, M.J., Ma, L., Ahuja, L.R., Rojas, K.W., 2001. Macropore
component assessment of the Root Zone Water Quality Model (RZWQM)
using no-till soil blocks. Trans. ASAE 44, 843–852.

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35
Malone, R.W., Ahuja, L.R., Ma, L.W., Wauchope, R.D., Ma, Q.L., Rojas, K.W.,
2004. Application of the Root Zone Water Quality Model (RZWQM) to
pesticide fate and transport: an overview. Pest Manag. Sci. 60 (3),205–221.
Maraqa, M.A., 2001. Prediction of mass-transfer coefﬁcient for solutetransport in porous media. J. Cont. Hydrol. 50, 1–19.
Mayer, K.U., Frind, E.O., Blowes, D.W., 2002. A numerical model for theinvestigation of reactive transport in variably saturated media using ageneralized formulation for kinetically controlled reactions. Water
Resour. Res. 38, 1174–1194.
McGuire, K.J., Weiler, M., McDonnell, J.J., 2007. Integrating tracer experimentswith modeling to assess runoff processes and water transit times. Adv.
Water Resour. 30, 824–837.
Merdun, H., Quisenberry, V.L., 2004. Simulation of water and solute transportwith MACRO model in Cecil loamy sand soil. Australian J. Soil Res. 42 (8),939–951.
Messing, I., Jarvis, N.J., 1990. Seasonal variation in ﬁeld-saturated hydraulicconductivity in two swelling clay soils in Sweden. J. Soil Sci. 41, 229–237.
Messing, I., Jarvis, N.J.,1993. Temporal variation in the hydraulic conductivity of atilled clay soil as measured by tension inﬁltrometers. J. Soil Sci. 44, 11–24.
Meyer-Windel, S., Lennartz, B., Widmoser, P., 1999. Bromide and herbicidetransport under steady-state and transient ﬂow conditions. J. Eur. Soil Sci.
50, 23–33.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., van Genuchten, M.T., 1997.
New piecewise-continuous hydraulic functions for modelling preferential ﬂow in an intermittent-ﬂood-irrigated ﬁeld. Water Resour. Res. 33(9), 2049–2063.
Molina, J.A.E., Richards, K., 1984. Simulation models of the nitrogen andcarbon cycle in the soil–water–plant system, NCSWAP: guide for the
preparation of input data ﬁles and execution of NCSWAP. Soil Ser., 116.
Dept. of Soil Sci., Univ. of Minnesota, St. Paul.
Mori, Y., Maruyama, T., Mitsuno, T., 1999. Soft X-ray radiography of drainage
patterns of structured soils. Soil Sci. Soc. Am. J. 63, 733–740.
Mulungu, D.M.M., Ichikawa, Y., Shiiba, M., 2005. A physically based
distributed subsurface–surface ﬂow dynamics model for forested
mountainous catchments. Hydrol. Processes 19 (20), 3999–4022.
Neuman, S.P., 1972. Finite Element Computer Programs for Flow in Saturated–
Unsaturated Porous Media. 2nd Annu. Rep. Project A10-SWC-77.
Hydraulic Eng. Lab., Technion, Haifa, Israel.
Nieber, J.L., Misra, D., 1995. Modeling ﬂow and transport in heterogeneous,dual-porosity drained soils. Irrig. Drain. Syst. 9, 217–237.
Nieber, J.L., Friedel, M.J., Munir, H.M., 1994. VARSAT2D: Finite-Element
Analysis of Variably Saturated Two Dimensional Flow. Information
Circular, vol. 9373. Bureau of Mines, USDI, Washington, DC.
Noguchi, S., Tsuboyama, Y., Sidle, R.C., Hosoda, I., 1999. Morphologicalcharacteristics of macropores and the distribution of preferential ﬂow
pathways in a forested slope segment. Soil Sci. Soc. Am. J. 63, 1413–1423.
Novák, V., Šimůnek, J., van Genuchten, M.Th., 2000. Inﬁltration of water into
soils with cracks, ASCE J. Irrig. And Drain. Engineering 126 (1), 41–47.
Olsson, J., Persson, M., Albergel, J., Berndtsson, R., Zante, P., Öhrström, P., Nasri, S.,
2002. Multiscaling analysis and random cascade modelling of dye inﬁltration.
Water Resour. Res. 38 (11), 1263. doi:10.1029/2001WR000880.
Olsson, J., Persson, M., Jinno, K., 2007. Analysis and modeling of solutetransport dynamics by breakdown coefﬁcients and random cascades.
Water Resour. Res. 43. doi:10.1029/2005WR004631 W03417.
Or, D., Tuller, M., 2000. Flow in unsaturated fractured porous media: hydraulicconductivity of rough surfaces. Water Resour. Res. 36, 1165–1177.
Or, D., Tuller, M., 2003. Hydraulic conductivity of partially saturated fractured
porous media: ﬂow in a cross-section. Adv. Water Resour. 26, 883–898.
Othmer, H.B., Diekkrüger, B., Kutilek, M., 1991. Bimodal porosity andunsaturated hydraulic conductivity. Soil Science 152, 139–150.
Refsgaard, J.C., Storm, B., 1995. MIKE SHE. In: Singh, V.P. (Ed.), Computer
Models of Watershed Hydrology. Water Resources Publications, Colorado, USA, pp. 809–846.
Ross, P.J., Smettem, R.J., 1993. Describing soil hydraulic properties with sumsof simple functions. Soil Sci. Soc. Am. J. 57, 26–29.
Pachepsky, Y.A., Rawls, W.J., Lin, H.S., 2006. Hydropedology and pedotransfer
functions. Geoderma 131 (3–4), 308–316.
Panday, S., Huyakorn, P.S., 2004. A fully coupled physically-based spatiallydistributed model for evaluating surface/subsurface ﬂow. Adv. Water
Resour. 27, 361–382.
Park, Y.-J., Sudicky, E.A., Panday, A., Sykes, J.F., Guvanasen, V., 2008.
Application of implicit sub-time stepping to simulate ﬂow and transport
in fractured porous media. Adv. Water Resour. 31 (7), 995–1003.
Parker, J.C., Valocchi, A.J., 1986. Constraints on the validity of equilibrium andﬁrst-order kinetic transport models in structured soils. Water Resour.
Res. 223, 399–407.
Perret, J., Prasher, S.O., Kantzas, A., Langford, C., 1999. Three-dimensional
quantiﬁcation of macropore networks in undisturbed soil columns. Soil
Sci. Soc. Am. J. 63, 1530–1543.

33

Perret, J., Prasher, S.O., Kantzas, A., Langford, C., 2000. A two-domain
approach using CAT scanning to model solute transport in soil. J. Environ.
Qual. 29 (3), 995–1010.
Pot, V., Šimůnek, J., Benoit, P., Coquet, Y., Yra, A., Martínez-Cordón, M.-J., 2005.
Impact of rainfall intensity on the transport of two herbicides inundisturbed grassed ﬁlter strip soil cores. J. Contam. Hydrol. 81, 63–88.
Priesack, E., Durner, W., 2006. Closed-form expression for the multi-modal
unsaturated conductivity function. Vadose Zone J. 5, 121–124.
Pruess, K., 2004. The TOUGH Codes—a family of simulation tools for
multiphase ﬂow and transport processes in permeable media. Vadose
Zone J. 3, 738–746.
Rao, P.S.C., Jessup, R.E., Rolston, E., Davidson, J.M., Kilcrease, D.P., 1980a.
Experimental and mathematical description of nonadsorbed solutetransfer by diffusion in spherical aggregates. Soil Sci. Soc. Am. J. 44,
684–688.
Rao, P.S.C., Rolston, D.E., Jessup, R.E., Davidson, D.M., 1980b. Solute transport
in aggregated porous media: theoretical and experimental evaluation.
Soil Sci. Soc. Am. J. 44, 1139–1146.
Rappoldt, C., Verhagen, J.H.G., 1999. Equivalent cylinder systems representing
the soil matrix in diffusion-reaction models for an aggregated soil.
Transport Porous Media 37 (1), 1–24.
Rieu, M., Sposito, G., 1991. Fractal fragmentation, soil porosity, and soil water
properties. II. Applications. Soil Sci. Soc. Am. J. 55, 1239–1244.
Rosqvist, N.H., Dollar, L.H., Fourie, A.B., 2005. Preferential ﬂow in municipal
solid waste and implications for long-term leachate quality: valuation oflaboratory-scale experiments. Waste Manag. Res. 23, 1–14.
Ross, P.J., Smettem, K.R.J., 2000. A simple treatment of physical nonequilibrium water ﬂow in soils. Soil Sci. Soc. Am. J. 64 (6), 1926–1930.
Roulier, S., Jarvis, N., 2003. Modelling macropore ﬂow effects on pesticideleaching: inverse parameter estimation using microlysimeters. J. Environ.
Qual. 32 (6), 2341–2353.
Schwartz, R.C., Juo, A.S.R., McInnes, K.J., 2000. Estimating parameters for adual-porosity model to describe non-equilibrium, reactive transport in aﬁne-textured soil. J. Hydrol. 229 (3–4), 149–167.
Seo, Y., Lee, J., 2005. Characterizing preferential ﬂow of nitrate and phosphatein soil using time domain reﬂectometry. Soil Sci. 170 (1), 47–54.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 2001. ROSETTA: a computer
program for estimating soil hydraulic parameters with hierarchicalpedotransfer functions. J. Hydrol. 251, 163–176.
Scorza Júnior, R.P., Jarvis, N.J., Boesten, J.J.T.I., van der Zee, S.E.A.T.M., Roulier,S., 2007. Testing MACRO (version 5.1) for pesticide leaching in a Dutch
clay soil. Pest Management Science 63, 1011–1025.
Shaw, J.N., West, L.T., Radcliffe, D.E., Bosch, D.D., 2000. Preferential ﬂow andpedotransfer functions for transport properties in sandy Kandiudults.
Soil Sci. Soc. Am. J. 64, 670–678.
Sidle, R.C., Noguchi, S., Tsuboyama, Y., Laursen, K., 2001. A conceptual modelof preferential ﬂow systems in forested hillslopes: evidence of selforganization. Hydrol. Processes 15 (10), 1675–1692.
Šimůnek, J., Vogel, T., van Genuchten, M.Th., 1994. The SWMS_2D Code forsimulating water ﬂow and solute Transport in two dimensional variably
saturated media (Version 1.21). Research Report No. 132, U.S. Salinity
Laboratory, USDA, ARS, Riverside, California.
Šimůnek, J., Huang, K., van Genuchten, M.Th., 1995. The SWMS_3D Code forSimulating Water Flow and Solute Transport in Three-Dimensional
Variably Saturated Media. Version 1.0, Research Report No. 139, U.S.
Salinity Laboratory, USDA, ARS, Riverside, California, 155 pp.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1998. The HYDRUS-1D software
package for simulating the movement of water, heat, and multiple solutes
in variably-saturated media. Version 2.0. IGWMC-TPS-70. International
Ground Water Modelling Center, Colorado School of Mines, Golden, CO.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1999. The HYDRUS-2D software
package for simulating two-dimensional movement of water, heat, andmultiple solutes in variably-saturated media. Version 2.0. IGWMC-TPS53. International Ground Water Modelling Center, Colorado School ofMines, Golden, CO.
Šimůnek, J., Wendroth, O., Wypler, N., van Genuchten, M.T., 2001. Nonequilibrium water ﬂow characterized by means of upward inﬁltration
experiments. Europ. J. Soil Sci. 52 (1), 13–24.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.Th., Gärdenäs, A., 2003. Review andcomparison of models for describing non-equilibrium and preferential
ﬂow and transport in the vadose zone. J. Hydrol. 272, 14–35.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2005. The HYDRUS-1D software
package for simulating the one-dimensional movement of water, heat,and multiple solutes in variably-saturated media. Version 3.0, HYDRUS
Software Series 1. Department of Environmental Sciences, University ofCalifornia Riverside, Riverside, CA. 270 pp.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2006. The HYDRUS Software
Package for Simulating Two- and Three-Dimensional Movement of
Water, Heat, and Multiple Solutes in Variably-Saturated Media, Technical
Manual, Version 1.0. PC Progress, Prague, Czech Republic, p. 241.

34

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35

Šimůnek, J., van Genuchten, M.Th., 2008. Modelling nonequilibrium ﬂow andtransport with HYDRUS. Vadose Zone J. 7 (2), 782–797 Special Issue“Vadose Zone Modelling”.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.Th., 2008a. The
HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media.
Version 4.0, HYDRUS Software Series 4. Department of Environmental
Sciences, University of California Riverside, Riverside, CA, USA. 315 pp.
Šimůnek, J., van Gencuhten, M.Th., Šejna, M., 2008b. Development andapplications of the HYDRUS and STANMOD Software Packages andRelated Codes. Vadose Zone J. 7, 587–600.
Šimůnek, J., Köhne, J.M., Kodešová, R., Šejna, M., 2008c. Simulating nonequilibrium
movement of water, solutes and particles using HYDRUS — a review of recent
applications. Soil & Water Res. 3, S42–S51 2008 (Special Issue 1).
Sinkevich, M.G., Walter, M.T., Lembo, A.J., Richards, B.K., Peranginangin, N.,
Aburime, S.A., Steenhuis, T.S., 2005. A GIS-based ground water contamination risk assessment tool for pesticides. Ground Water Monitor.
Remediation 25 (4), 82–91.
Slough, K.J., Sudicky, E.A., Forsyth, P.A., 1999. Grid reﬁnement for modelling
multiphase ﬂow in discretely fractured porous media. Adv. Water Resour.
23, 261–269.
Soulsby, C., Neal, C., Laudon, H., Burns, D., Merot, P., Bonell, M., Dunn, S.M.,
Tetzlaff, D., 2008. Catchment Data for Process Conceptualization: Simply
Not Enough? Hydrol. Process. 22, 2057–2061.
Stagnitti, F., Allinson, G., Morita, M., Nishikawa, M., Ii, H., Hirata, T., 2000.
Temporal moments analysis of preferential solute transport in soils.
Environ. Model. Assessment 5 (4), 229–236.
Stagnitti, F., Li, L., Barry, A., Allinson, G., Parlange, J.-Y., Steenhuis, W.,
Lakshmanan, E., 2001. Modelling solute transport in structured soils:
performance evaluation of the ADR and TRM models. Math. Comput.
Model. 34, 433–440.
Sudicky, E.A., McLaren, R., 1992. Numerical analysis of solute migration
through fractured clayey deposits into underlying aquifers. Water
Resour. Res. 28 (2), 515–526.
Szymkiewicz, A., Lewandowska, J., Angulo-Jaramillo, R., Butlańska, J., 2008. Twoscale modelling of unsaturated water ﬂow in a double-porosity medium
under axi-symmetric conditions. Can. Geotech. J. 45 (2), 238–251.
doi:10.1139/T07-096.
Tetzlaff, D., McDonnell, J.J., Uhlenbrook, S., McGuire, K.J., Bogaart, P.W., Naef,
F., Baird, A.J., Dunn, S.M., Soulsby, C., 2008. Conceptualizing catchment
processes: simply too complex? Hydrol. Proc. 22, 1727–1730.
Therrien, R., Sudicky, E.A., 1996. Three-dimensional analysis of variablysaturated ﬂow and solute transport in discretely-fractured porous media.
J. Contam. Hydrol. 23, l–44.
Therrien, R., McLaren, R.G., Sudicky, E.A., Panday, S.M., 2005. HydroGeoSphere: A
Three-Dimensional Numerical Model Describing Fully-Integrated Subsurface
and Surface Flow and Solute Transport. Groundwater Simulations Group,
Univ. of Waterloo, Waterloo, Ont., Canada. 322 pp.
Tiemeyer, B., Moussa, R., Lennartz, B., Voltz, M., 2007. MHYDAS-DRAIN: a spatiallydistributed model for small, artiﬁcially drained lowland catchments. Eco.
Model. 209 (1), 2–20.
Tokunaga, T.K., Wan, J., Sutton, S.R., 2000. Transient ﬁlm ﬂow on roughfracture surfaces. Water Resour. Res. 36, 1737–1746.
Toride, N., Leij, F.J., van Genuchten, M.Th., 1999. The CXTFIT code for estimatingtransport parameters from laboratory or ﬁeld tracer experiments. Version
2.1. Res. Rep. No 137. U.S. Salinity Laboratory. ARS-USDA, Riverside, CA.
Tromp-van Meerveld, H.J., McDonnell, J.J., 2006a. Threshold relations insubsurface stormﬂow 1: a 147 storm analysis of the Panola hillslope
trench. Water Resour. Res. 42. doi:10.1029/2004WR003778.
Tromp-van Meerveld, H.J., McDonnell, J.J., 2006b. Threshold relations insubsurface stormﬂow 2: the ﬁll and spill hypothesis: an explanation forobserved threshold behaviour in subsurface stormﬂow. Water Resour.
Res. 42. doi:10.1029/2004WR003800.
Tsutsumi, D., Sidle, R.C., Kosugi, K., 2005. Development of a simple lateralpreferential ﬂow model with steady state application in hillslope soils.
Water Resour. Res. 41 (12) Art. No. W12420.
Tuller, M., Or, D., 2002. Unsaturated hydraulic conductivity of structured
porous media: a review of liquid conﬁguration-based models. Vadose
Zone J. 1, 14–37.
Uchida, T., Kosugi, K., Mizuyama, T., 2001. Effects of pipeﬂow on hydrologicalprocesses and its relation to landslides: a review of pipeﬂow studies inforested headwater catchments. Hydrol. Processes 15, 2151–2174.
Van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabat, P., Van Walsum,
P.E.V., Groenendijk, P., Van Diepen, C.A., 1997. Theory of SWAP Version
2.0. Simulation of water ﬂow. Solute Transport and Plant Growth in the
Soil–Water–Atmosphere–Plant Environment, Technical Document 45,
DLO Winand Staring Centre, Report 71. Department of Water Resources,
Wageningen Agricultural University, Wageningen. 167 pp.
van Dam, J.C., 2000. Simulation of ﬁeld-scale water ﬂow and bromidetransport in a cracked clay soil. Hydrol. Proc. 14 (6), 1101–1117.

Van der Hoven, S.J., Solomon, D.K., Moline, G.R., 2002. Numerical simulationof unsaturated ﬂow along preferential pathways: implications for the useof mass balance calculations for isotope storm hydrograph separation.
J. Hydrol. 268 (1), 214–233.
Van Genuchten, M.Th., 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Genuchten, M.Th., Wierenga, P.J., 1976. Mass transfer studies in sorbing
porous media. I. analytical solutions. Soil Sci. Soc. Am. J. 40, 473–480.
van Genuchten, M.Th., Dalton, F.N., 1986. Models for simulating salt
movement in aggregated ﬁeld soils. Geoderma 38, 165–183.
Vanderborght, J., Vereecken, H., 2007. Review of dispersivities for transportmodeling in soils. Vadose Zone J. 6, 29–52.
Vanderborght, J., Mallants, D., Vanclooster, M., Feyen, J., 1997. Parameter
uncertainty in the mobile–immobile solute transport model. J. Hydrol.
190 (1–2), 75–101.
Vanderborght, J., Timmerman, A., Feyen, J., 2000. Solute transport for steadystate and transient ﬂow in soils with and without macropores. Soil Sci.
Soc. Am. J. 64, 1305–1317.
Vanderborght, J., Vanclooster, M., Timmerman, A., Seuntjens, P., Mallants, D.,
Kim, D.-J., et al., 2001. Overview of inert tracer experiments in key
Belgian soil types: relation between transport and soil morphologicaland hydraulic properties. Water Resour. Res. 37, 2873–2888.
Vanderborght, J., Gähwiller, P., Flühler, H., 2002. Identiﬁcation of transport
processes in soil cores using ﬂuorescent tracers. Soil Sci. Soc. Am. J. 66,
774–787.
Vanderborght, J., Kasteel, R., Vereecken, H., 2006. Stochastic continuum
transport equations for ﬁeld-scale solute transport: overview of theoreticaland experimental results. Vadose Zone J. 5, 184–203.
VanderKwaak, J.E., 1999. Numerical Simulation of Flow and Chemical
Transport in Integrated Surface–Subsurface Hydrologic Systems, Ph.D.
Thesis, 217 pp., Univ. of Waterloo, Waterloo, Ont., Canada.
VanderKwaak, J.E., Loague, K., 2001. Hydrologic-response simulations for theR-5 catchment with a comprehensive physics-based model. Water
Resour. Res. 37, 999–1013.
Vereecken, H., Hubbard, S., Binley, A., Ferre, T., 2004. Hydrogeophysics: anintroduction from the guest editors. Vadose Zone J. 3, 1060–1062.
Vervoort, R.W., Radcliffe, D.E., West, L.T., 1999. Soil structural development
and preferential solute ﬂow. Water Resour. Res. 35, 913–928.
Vogel, H.-J., Roth, K., 1998. A new approach for determining effective soil
hydraulic functions. Europ. J. Soil Sci. 49, 547–556.
Vogel, H.-J., Roth, K., 2001. Quantitative morphology and network representation of soil pore structure. Adv. Water Resour. 24, 233–242.
Vogel, H.-J., 2002. Topological characterization of porous media. In: Mecke, K.,
Stoyan, D. (Eds.), Morphology and Condensed Matter — Physics andGeometry of Spatially Complex Systems. Lecture Notes in Physics, 600,
pp. 75–92.
Vogel, H.J., Roth, K., 2003. Moving through scales of ﬂow and transport in soil.
J. Hydrol. 272, 95–106.
Vogel, H.J., Cousin, I., Ippisch, O., Bastian, P., 2006. The dominant role ofstructure for solute transport in soil: experimental evidence andmodelling of structure and transport in a ﬁeld experiment. Hydrol.
Earth System Sci. 10 (4), 495–506.
Vogel, H.J., Ippisch, O., 2008. Estimation of a critical spatial discretization limit
for solving the Richards' equation at large scales. Vadose Zone J. 7 (1),
112–114.
Vogel, T., Huang, K., Zhang, R., van Genuchten, M.Th., 1996. The HYDRUS codefor simulating one-dimensional water ﬂow, solute transport, and heatmovement in variably-saturated media (version 5.0). Research Report
No. 140. U. S. Salinity Laboratory, USDA-ARS, Riverside, CA.
Vogel, T., Gerke, H.H., Zhang, R., van Genuchten, M.Th., 2000. Modelling ﬂow
and transport in a two-dimensional dual-permeability system withspatially variable hydraulic properties. J. Hydrol. 238, 78–89.
Vrugt, J.A., Bouten, W., Weerts, A.H., 2001. Information content of data for
identifying soil hydraulic parameters from outﬂow experiments. Soil Sci.
Soc. Am. J. 65, 19–27.
Vrugt, J.A., Schoups, G., Hopmans, J.W., Young, C., Wallender, W.W., Harter, T.,
Bouten, W., 2004. Inverse modelling of large-scale spatially distributed
vadose zone properties using global optimization. Water Resour. Res. 40 (6)Art. No. W06503.
Wallach, R., Parlange, J.-Y., 1998. Modeling transport in a single crack by thedual-porosity concept with a boundary layer at the interface. J. Contam.
Hydrol. 34, 121–138.
Wallach, R., Parlange, J.Y., 2000. Applying the boundary layer concept tomodel transport of dissolved chemicals in preferential ﬂow paths. Water
Resour. Res. 36 (10), 2845–2851.
Wang, J.S.Y., Narasimhan, T.N., 1985. Hydrologic mechanisms governing ﬂuid
ﬂow in a partially saturated fractured porous medium. Water Resour. Res.
21 (12), 1861–1874.
Wang, Z., Jury, W.A., Tuli, A., Kim, D.J., 2004. Unstable ﬂow during redistribution:
controlling factors and practical implications. Vadose Zone J. 3 (2), 549–559.

J.M. Köhne et al. / Journal of Contaminant Hydrology 104 (2009) 4–35
Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs.
Soc. Petrol. Eng. J. 3, 245–255.
Weiler, M., 2005. An inﬁltration model based on ﬂow variability in macropores:
development, sensitivity analysis and applications. J. Hydrol. 310 (1–4),
294–315.
Weiler, M., Naef, F., 2003. Simulating surface and subsurface initiation ofmacropore ﬂow. J. Hydrol. 273, 139–154.
Weiler, M., McDonnell, J., 2004. Virtual experiments: a new approach for
improving process conceptualization in hillslope hydrology. J. Hydrol.
285, 3–18.
Weiler, M., McDonnell, J.J., 2006. Testing nutrient ﬂushing hypotheses at thehillslope scale: a virtual experiment approach. J. Hydrol. 319, 339–356.
Weiler, M., McDonnell, J.R.J., 2007. Conceptualizing lateral preferential ﬂow
and ﬂow networks and simulating the effects on gauged and ungauged
hillslopes. Water Resour. Res. 43, W03403. doi:10.1029/2006WR004867.
Wigmosta, M.S., Vail, L.W., Lettenmaier, D.P., 1994. A distributed hydrologyvegetation model for complex terrain. Water Resour. Res. 30, 1665–1679.
Wilson, G.V., Jardine, P.M., Gwo, J.P., 1992. Modeling the hydraulic propertiesof a multi-region soil. Soil Sci. Soc. Am. J. 56, 1731–1737.
Wooding, R., 1968. Steady inﬁltration from a shallow circular pond. Water
Resour. Res. 4 (6), 1259–1273.
Wöhling, Th., Vrugt, J., Barkle, G.F., 2008. Comparison of three multiobjective
optimization algorithms for inverse modeling of vadose zone hydraulicproperties. Soil Sci. Soc. Am. J. 72, 305–319.

35

Wu, Y.S., Pruess, K., 2000. Numerical simulation of non-isothermal multiphase tracer transport in heterogeneous fractured porous media. Adv.
Water Resour. 23 (7), 699–723.
Young, D.F., Ball, W.P., 1997. Effects of column conditions on the ﬁrst-order rate
modelling of nonequilibrium solute breakthrough: cylindrical macroporesversus spherical media. Water Resour. Res. 33 (5), 1149–1156.
Young, D.F., Ball, W.P., 1998. Estimating diffusion coefﬁcients in lowpermeability porous media using a macropore column. Environ. Sci. &Technol. 32 (17), 2578–2584.
Zehe, E., Maurer, T., Ihringer, J., Plate, E., 2001. Modelling water ﬂow and masstransport in a loess catchment. Phys. Chem. Earth Part B-Hydrol. Oceans
Atmosph. 26 (7–8), 487–507.
Zehe, E., Blöschl, G., 2004. Predictability of hydrologic response at the plot
and catchment scales: role of initial conditions. Water Resour. Res. 40,
W10202. doi:10.1029/2003WR002869.
Zehe, E., Becker, R., Bardossy, A., Plate, E., 2005. Uncertainty of simulated
catchment runoff response in the presence of threshold processes: roleof initial soil moisture and precipitation. J. Hydrol. 315 (1–4), 183–202.
Zhang, P.F., Šimůnek, J., Bowman, R.S., 2004. Nonideal transport of solute andcolloidal tracers through reactive zeolite/iron pellets. Water Resour. Res.
40 (4), W04207.
Zumr, D., Dohnal, M., Hrncir, M., Cislerovka, M., Vogel, T., Dolezal, F., 2006.
Simulation of soil water dynamics in structured heavy soils with respect
to root water uptake. Biologia 61, S320–S323 Suppl. 19.