Journal of Hydrology 272 (2003) 14–35
www.elsevier.com/locate/jhydrol

Review and comparison of models for describing non-equilibrium
and preferential flow and transport in the vadose zone
Jirka Šimůneka,*, Nick J. Jarvisb,1, M.Th. van Genuchtena, Annemieke GärdenäsbabGeorge E. Brown Jr. Salinity Laboratory, USDA-ARS, 450 West Big Springs Road, Riverside, CA 92507, USA
Department of Soil Sciences, Swedish University of Agricultural Sciences, P.O. Box 7014, 750 07 Uppsala, Sweden
Revised 19 April 2002

Abstract
In this paper, we review various approaches for modeling preferential and non-equilibrium flow and transport in the vadosezone. Existing approaches differ in terms of their underlying assumptions and complexity. They range from relatively simplistic
models to more complex physically based dual-porosity, dual-permeability, and multi-region type models. A relatively simple
dual-porosity flow model results when the Richards equation is combined with composite (double-hump type) equations for thehydraulic properties to account for both soil textural (matrix) and soil structural (fractures, macropores, peds) effects on flow.
The simplest non-equilibrium flow model, a single-porosity model which distinguishes between actual and equilibrium watercontents, is based on a formulation by Ross and Smettem [Soil Sci. Soc. Am. J. 64 (2000) 1926] that requires only one
additional parameter to account for non-equilibrium. A more complex dual-porosity, mobile– immobile water flow model
results when the Richards or kinematic wave equations are used for flow in the fractures, and immobile water is assumed to exist
in the matrix. We also discuss various dual-permeability models, including the formulation of Gerke and van Genuchten [Water
Resour. Res. 29 (1993a) 305] and the kinematic wave approach as used in the MACRO model of Jarvis [Technical Description
and Sample Simulations, Department of Soil Science, Swedish University of Agricultural Science, Uppsala, Sweden (1994)51]. Both of these models invoke terms accounting for the exchange of water and solutes between the matrix and the fractures.
Advantages and disadvantages of the different models are discussed, and the need for inter-code comparison is stressed,
especially against field data that are sufficiently comprehensive to allow calibration/validation of the more complex models andto distinguish between alternative modeling concepts. Several examples and comparisons of equilibrium and various nonequilibrium flow and transport models are also provided.
Published by Elsevier Science B.V.
Keywords: Preferential flow; Non-equilibrium flow; Vadose zone; Dual-porosity models; Dual-permeability models; Solute transport; Pedon
scale

1. Introduction* Corresponding author. Tel.: þ1-909-369-4865; fax: þ 1-909342-4964.
E-mail addresses: jsimunek@ussl.ars.usda.gov (J. Šimůnek),
nick.jarvis@mv.slu.se (N.J. Jarvis).
1Tel.: þ46-18-67-2465; fax: þ 46-18-67-2795.
0022-1694/03/$ - see front matter Published by Elsevier Science B.V.
PII: S 0 0 2 2 - 1 6 9 4 ( 0 2 ) 0 0 2 5 2 - 4

The problem of non-equilibrium and/or preferential flow and transport has received much attention in
the soil and agricultural sciences because of its
implications in accelerating the movement of agricultural contaminants (fertilizers, pesticides,

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

pathogenic microorganisms, toxic trace elements)
through the unsaturated zone to underlying groundwater. The potentially rapid migration of radionuclides from low- and high-level nuclear waste
disposal facilities, and the preferential movement of
non-aqueous liquids or other pollutants from underground storage tanks, waste disposal sites and minetailings, has also become a concern for hydrologists,
geophysicists, and environmental scientists.
Preferential and non-equilibrium flow and transport are probably the most frustrating processes interms of hampering accurate predictions of contaminant transport in soils and fractured rocks. Preferential flow, as opposed to uniform flow, results inirregular wetting of the soil profile as a directconsequence of water moving faster in certain partsof the soil profile than in others. Hendrickx and Flury(2001) defined preferential flow as ‘all phenomena
where water and solutes move along certain pathways, while bypassing a fraction of the porous
matrix’. Thus, an important characteristic of preferential flow is that during wetting, part of the moisture
front can propagate quickly to significant depths
while bypassing a large part of the matrix pore-space.
Water and solutes may move to far greater depths,
and much faster, than would be predicted with theRichards equation using area-averaged moisture
contents and pressure heads (Beven, 1991). Anotherimportant characteristic of preferential (non-uniform)
flow is its non-equilibrium nature. Even for uniform
flow conditions, most of the water and its dissolved
solutes generally move through the largest continuous pores that are filled with water at a particulartension. This is reflected in the shape of the highly
non-linear hydraulic conductivity function, which
typically shows dramatic increases with increasingwater contents, particularly as the larger pores
become active. While conditions at or close toequilibrium exist between the different types of poresin a soil during uniform flow, this is generally not thecase during preferential flow. Following the ideas
developed by Skopp (1981) and Flühler et al. (1996),
Jarvis (1998) considered non-equilibrium to be themost important feature of preferential flow bydefining it as a flow regime in which ‘for various
reasons, infiltrating water does not have sufficient
time to equilibrate with slowly moving resident waterin the bulk of the soil matrix’.

15

The presence of macropores and other structuralfeatures, development of flow instabilities (i.e.
fingering) caused by profile heterogeneities or water
repellency (Hendrickx et al., 1993), and funneling offlow due to the presence of sloping soil layers that
redirect downward water flow are probably the mostimportant causes of preferential flow. While the latter
two processes, i.e. flow instability and funneling, areusually caused by textural differences and otherfactors at scales significantly larger than the porescale, macropore flow and transport are usually
generated at pore or slightly larger scales, including
scales where soil structure first manifests itself (i.e.
the pedon scale). Following a recent report by NRC(2001), we will use the term granular medium to referto those media in which flow and transport are notsignificantly affected by the presence of macroporesand fractures, and in which infiltration produces
relatively uniform wetting fronts.
The last two decades have seen the development ofa relatively large number of models that consider nonequilibrium flow. In this paper, we review several
popularly used approaches to modeling preferential
and non-equilibrium flow and transport in the vadosezone starting with relatively simplistic approaches,
followed by models of increasing complexity. While
Hendrickx and Flury (2001) discussed preferential
flow mechanisms and processes at various scales,from pore to areal scales, in this paper we will focus
mainly on models describing macropore and nonequilibrium flow resulting from processes at the poreto pedon scales.

2. Conceptual models
Preferential flow in structured media (both macroporous soils and fractured rocks) can be describedusing a variety of dual-porosity, dual-permeability,
multi-porosity, and/or multi-permeability models(Pruess and Wang, 1987; Gerke and van Genuchten,
1993a; Gwo et al., 1995; Jarvis, 1998). Dual-porosityand dual-permeability models both assume that theporous medium consists of two interacting regions,
one associated with the inter-aggregate, macropore, or
fracture system, and one comprising micropores (orintra-aggregate pores) inside soil aggregates or therock matrix. While dual-porosity models assume that

16

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

water in the matrix is stagnant, dual-permeability
models allow for water flow in the matrix as well.
Dual-porosity models have long been applied tosolute transport studies. Especially popular early on
were dual-porosity models in which distinct mobile
and immobile flow regions are assumed to be present(van Genuchten and Wierenga, 1976). Dual-permeability models in which water can move in boththe inter- and intra-aggregate pore regions are now
also becoming more popular (Pruess and Wang, 1987;
Gerke and van Genuchten, 1993a; Jarvis, 1994).
Available dual-permeability models differ mainly in
how they implement water flow in and between thetwo pore regions, especially with respect to the degreeof simplification and empiricism. Approaches tocalculating water flow in macropores or interaggregate pores range from those invoking Poiseuille’s equation (Ahuja and Hebson, 1992), the Green
and Ampt or Philip infiltration models (Ahuja andHebson, 1992; Chen and Wagenet, 1992), thekinematic wave equation (Germann, 1985; Germann
and Beven, 1985; Jarvis, 1994), and the Richards
equation (Gerke and van Genuchten, 1993a). Multiporosity and/or multi-permeability models are basedon the same concept as dual-porosity and dualpermeability models, but include additional interacting pore regions (Gwo et al., 1995; Hutson andWagenet, 1995). They can be simplified immediately
to dual-porosity/permeability models.

3. Equilibrium flow and transport models
Process-based models for flow and transport in
granular (or single-porosity) media are generally
based on the Richards equation for variably saturated
water flow and the convection– dispersion equation
for solute transport, i.e.:



›u
›
›h
¼
KðhÞ
þ1 2S
ð1Þ
›z
›t
›z


›uc
›rs
›
›c
›qc
þ
¼
2 mðuc þ rsÞ
uD
2
›z
›t
›t
›z
›z
þ gu þ gr

ð2Þ

respectively, where z is the vertical coordinate(positive upwards), t is time, S is a sink term, h is

the pressure head, u is the water content, K is theunsaturated hydraulic conductivity function, c and s
are solute concentrations in the liquid and solid
phases, respectively, q is the volumetric flux density,
m is a first-order rate constant, g is a zero-order rate
constant, r is the soil bulk density, and D is thedispersion coefficient. For convenience, we will lump
below all zero- and first-order rate reactions into one
term, f. Eq. (1) is invariably accompanied by a set of
constitutive relations characterizing the unsaturated
soil hydraulic properties, specifically the retention
curve, uðhÞ; and the unsaturated hydraulic conductivity function, KðhÞ: Commonly used soil hydraulicmodels for granular media were reviewed by Leij et al.
(1997).
Although Eqs. (1) and (2) by themselves generally
describe uniform flow and transport, they often also
appear in various modifications, and/or have beencombined with other expressions or assumptions, toyield alternative approaches accounting for nonequilibrium preferential flow. Eqs. (1) and (2) are,for example, often used to describe flow and transport
in individual pore domains of dual- and multipermeability models.
3.1. Composite retention and hydraulic conductivity
functions
A very simple dual-porosity flow model results
when the Richards equation is combined with doublehump type composite equations for the hydraulicproperties. For example, Othmer et al. (1991) andDurner (1994) divided the porous medium into two(or more) overlapping regions and suggested to use
for each of these regions a van Genuchten – Mualem
type function (van Genuchten, 1980) of the soil
hydraulic properties. Linear superposition of thefunctions for each particular region gives the functions for the entire multi-modal pore system (Durner
et al., 1999):kXuðhÞ 2 ur
1
¼
wi
ð3Þ
ni mi
us 2 ur
ð1þlai hl Þ
i¼1
2! l X ki mikw a ½1 2 ð1 2 S1=m
X
ei Þ 
i¼1 i iwi SeiKðuÞ ¼ Ks
X k
2
i¼1
wa
i¼1 i iSe ðhÞ ¼

ð4Þ

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

where Se is the effective water content, and ur and us
denote the residual and saturated water contents,
respectively. The integer k denotes the number ofoverlapping subregions, wi are the weighting factors
for the sub-curves, and ai ; ni ; mi ð1 2 1=ni Þ; and l are
empirical parameters of the sub-curves. The hydrauliccharacteristics defined by Eqs. (3) and (4) contain4 þ 2k unknown parameters: ur ; us ; ai ; ni ; l, and Ks : Of
these, ur ; us ; and Ks have a physical meaning, whereas
ai ; ni ; and l are essentially empirical parameters
determining the shape of the retention and hydraulicconductivity functions (van Genuchten, 1980).
An example of the composite retention andhydraulic conductivity functions for two overlapping
porous media is shown in Fig. 1. Note that thepressure head axes are on a log scale, which causes thenear-saturated values to be significantly enlarged. In
this example, the fracture domain represents only
2.5% of the entire pore-space, but it accounts foralmost 90% of the hydraulic conductivity close tosaturation. Curves similar to those in Fig. 1 have beenused also for fractured rock by Peters and Klavetter(1988), Pruess and Wang (1987), and Flint et al.
(2001), among others.
One example of piecewise continuous (‘cut andjoin’ type) functions that account for a rapid increasein the hydraulic conductivity near saturation was
proposed by Mohanty et al. (1997) who used a van
Genuchten – Mualem model for the capillary-dominated flow domain and the following exponentialfunction to account for the non-capillary-dominatedflow domainpKðhÞ ¼ K p þ K p ½eðh2h Þd 2 1

hp , h # 0

ð5Þ

17

where hp is the critical or breakpoint soil waterpressure head where flow changes from capillarydominated to non-capillary-dominated flow, K p is thehydraulic conductivity corresponding to hp ; and d is afitting parameter representing effective macroporosity
or other structural features contributing to noncapillary-dominated flow. Very similar exponentialfunctions to account for the increase in hydraulicconductivity due to macropores in the near-saturated
water content range have been used by Smettem et al.
(1991) and Ross and Smettem (2000), among others.
A slightly different description was used byŠimůnek et al. (1998, 1999a) in the HYDRUS models
and its predecessors to account for rapid increases in
the conductivity near saturation as a result ofmacropore flow. Capillary-dominated flow was
described using the van Genuchten – Mualem modelfor the unsaturated hydraulic conductivity at lower
water contents until some point (hk ; Kk ), while a linear
increase in conductivity was assumed between thispoint and saturation to account for non-capillary,
macropore-dominated flow (Vogel and Cı́slerová,1988). Examples of applications that used composite
functions to account for increases in the hydraulicconductivity close to saturation are given by Smettem
et al. (1991), Mohanty et al. (1997, 1998), de Vos et al.
(1999), and Šimůnek and de Vos (1999).
Although composite hydraulic functions of thetype discussed earlier can account for a significant
increase in the hydraulic conductivity near saturation,
they do not by themselves lead to preferential flow.
When used in a single-domain model based on the
Richards equation, simulations will still produce
uniform wetting fronts, although possibly with someFig. 1. Example of composite retention (left) and hydraulic conductivity (right) functions (ur ¼ 0:00; us ¼ 0:50; a1 ¼ 0:01 cm21 ; n1 ¼ 1:50;
l ¼ 0:5; Ks ¼ 1 cm d21 ; w1 ¼ 0:975; w2 ¼ 0:025; a2 ¼ 1:00 cm21 ; n2 ¼ 5:00).

18

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

accelerated advance of the front for surface ponding
conditions because of the higher saturated conductivity. Such models, however, cannot predict preferential flow since there is no mechanism to account forlateral non-equilibrium in terms of water moving
primarily through larger pores or fractures andbypassing the matrix pore-space. Models based oncomposite hydraulic property functions implicitly
assume instantaneous equilibrium (i.e. instantaneous
exchange of water) between the two pore systems.
Also, double-hump type equations are typical of soils
having bimodal pore-size distributions, which can, butdoes not necessarily have to, be a dual-porosity
medium exhibiting non-equilibrium preferential flow.

4. Non-equilibrium flow and transport models
4.1. Single-porosity models
A relatively simple non-equilibrium flow model
was recently proposed by Ross and Smettem (2000).
Their model is based on observations in the literature(Schultze et al., 1999; Wildenschild et al., 2001) thatwater contents may sometimes differ greatly fromthose given by the equilibrium retention curve,
especially during rapid flow, and that equilibrium
may be reached only slowly as time proceeds. Ross
and Smettem (2000) took this non-equilibrium
process into account by combining the classical
Richards equation (1) with a kinetic description for
the water content approach towards equilibrium. By
doing so they removed from the Richards equation theassumption that water contents and pressure heads aretightly coupled through the equilibrium retention
curve. Since decoupling of the pressure head and thewater content in Eq. (1) now leads to two independent
variables (h and u ), they replaced the equilibrium
coupling assumption by an additional differential
equation:›u¼ f ð u; ue Þ›tð6Þwhere f ðu; ue Þ is a known function of the actual andequilibrium water contents. Ross and Smettem (2000)
assumed a simple linear driving function for f, i.e.
f ðu; ue Þ ¼ ðue 2 uÞ=t; where t is an equilibration time
constant. Substituting this function into Eq. (6),approximating the resulting expression by an implicit
finite difference equation, and incorporating this into
their numerical model produced the following
approximation (Ross and Smettem, 2000):u jþ1 ¼ u j þ ðu jþ1
2 u j Þ½1 2 expð2Dt=tÞeð7Þwhere the superscripts j þ 1 and j represent the new
and old time levels in the numerical time discretization and Dt is time step. Although one could
undoubtedly formulate more sophisticated non-equilibrium models, an important advantage of Eq. (7) isthat it requires only one additional parameter toaccount for non-equilibrium. Another advantage isthat it can be implemented easily into existing
variably saturated flow models, especially those
based on the mixed formulation of the Richards
equation. However, when we incorporated Eq. (7) inthe HYDRUS models (Šimůnek et al., 1998, 1999a)we experienced some numerical instabilities for
relatively large values of the time constant t when
the non-equilibrium moisture front moves rapidly
through the medium (i.e. when large changes in thepressure head produce only small changes in the watercontent). One disadvantage of this approach is that thelevel of preferential flow being produced is independent of the current water content (i.e. antecedent watercontent).
Ross and Smettem (2000) used Eq. (7) tosuccessfully simulate constant flux infiltration andoutflow responses for six large undisturbed cores. We
also used the model to successfully describe upward
infiltration data on small-undisturbed soil samples(Šimůnek et al., 2001) that exhibited severe nonequilibrium between measured pressure heads andvolumes of infiltrated water.
When using Ross and Smettem (2000) model, aninteresting question is how to best evaluate thehydraulic conductivity, i.e. as a function of
the pressure head, h, or the water content, u. Note,again, that u and h now are uncoupled and thus therelation between them is not unique. For the
calculations discussed later, we followed Ross andSmettem (2000) by calculating the hydraulic conductivity as a function of the pressure head. Thisimplicitly introduces the assumption that water will
move, contrary to assumptions inherent in the
Richards equation, first through the larger pores that

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

dominate the overall hydraulic conductivity. Calculating the hydraulic conductivity as a function of thewater content would significantly reduce the nonequilibrium features of this formulation.
Fig. 2 compares pressure head and water content
profiles calculated using the equilibrium (Eq. (1))
and Ross and Smettem non-equilibrium (Eq. (7))
water flow models. Van Genuchten (1980) soil
hydraulic parameters for a loamy soil (ur ¼ 0:02;
us ¼ 0:35; a1 ¼ 0:041 cm21 ; n1 ¼ 1:96; l ¼ 0:5;
Ks ¼ 0:000722 cm s21 ) were used in the simulations.
Ponding ðh ¼ 1 cmÞ was assumed at the soil surface
and the initial pressure head was set equal to
2 150 cm. A value of 3600 s was used for the time
constant, t, in the non-equilibrium model. Depending

19

upon the value of t, this model can produce wetting
fronts that will travel initially much faster, and to
greater depths (Fig. 2(b)), than those produced withthe equilibrium flow model (Fig. 2(a)). The wetting
front reached a depth of 35 cm after 1800 s for the
non-equilibrium model (Fig. 2(b)), as opposed to adepth of only 15 cm for the equilibrium model(Fig. 2(a)). This much deeper penetration is causedby the fact that flow initially is restricted to only asmall part of the flow domain (Fig. 2(b)). The
wetting front will eventually advance much slower
since most of the infiltrating water is later used forincreasing the water content behind the wetting front
towards equilibrium (Fig. 2(b)), rather than for more
downward vertical flow. For this reason, the moistureFig. 2. Pressure head and water content profiles calculated with the (a) equilibrium (Eq. (1)) and (b) non-equilibrium (Eq. (7)) (Ross andSmettem, 2000) water flow models. Thin lines in (b) represent final profiles ðt ¼ 7200 sÞ when the equilibration time constant, t, is decreased orincreased by 25%.

20

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

front advanced only 20 cm between 1800 and 7200 s,
for both the equilibrium and non-equilibrium models.
In order to demonstrate the sensitivity of the Ross
and Smettem model predictions to the equilibration
time constant, t, we included into Fig. 2(b) the final
pressure head and water content profiles ðt ¼ 7200 sÞ
obtained when we decreased or increased t by 25%(i.e. t ¼ 2700 and 4500 s, respectively).
When the non-equilibrium model of Ross andSmettem (2000) is coupled with the classical convection –dispersion equation (2), the solute front will
move much slower than the moisture front because of
the assumption that incoming water completely mixes
with water initially stored in the profile (Fig. 3(b)).
This complete mixing produces concentration profiles
that are very similar to those obtained using theequilibrium description for water flow (Fig. 3(a)).
Hence, to produce rapidly moving solute fronts, onemust also invoke assumptions for non-equilibriumtransport, such as the presence of immobile water(Fig. 3(c)). This feature will be discussed in a later
section. For the calculations shown in Fig. 3, waterflow was the same as for Fig. 2, the dispersivity was
assumed to be 2 cm, and the immobile water content
for case C was set at 0.05 cm3 cm23, with no transfer
of solute between the mobile and immobile regions.

macropores, and the subscript m to represent the soil
matrix, intra-aggregate pores, or the rock matrix.
The dual-porosity formulation for water flow canbe based on a mixed formulation of the Richards
equation (1) to describe water flow in the fractures anda mass balance equation to describe moisturedynamics in the matrix as follows:



›uf
›
›h
¼
KðhÞ
þ 1 2 S f 2 Gw ;
›z
›t
›z
ð9Þ
›um
¼ 2Sm þ Gw›twhere Sf and Sm are sink terms for both regions, and
Gw is the transfer rate for water from the inter- to theintra-aggregate pores.
An alternative approach for flow in the macropores
was suggested by Germann (1985) and Germann andBeven (1985). They used a kinematic wave equation
to describe gravitational movement of water inmacropores as follows:

›q
›q
þC
þ cr uf ¼ 0
›t
›z

ð10Þ

where r is the macropore sorbance, and where thevolumetric flux density, q, and the kinematic wave
velocity,C, are defined by:

4.2. Dual-porosity models

q ¼ buaf

ð11Þ

Dual-porosity models assume that water flow isrestricted to the fractures (or inter-aggregate pores andmacropores), and that water in the matrix (intraaggregate pores or the rock matrix) does not move atall. Thus, intra-aggregate pores represent immobile
pockets that can exchange, retain and store water, but
do not permit convective flow. This conceptualization
leads to two-region, dual-porosity type flow andtransport models (Phillip, 1968; van Genuchten and
Wierenga, 1976) that partition the liquid phase into
mobile (flowing, inter-aggregate), uf ; and immobile(stagnant, intra-aggregate), um ; regions:

C¼

›q
›uf

ð12Þ

u ¼ uf þ um

ð8Þ

with some exchange of water and/or solutes possible
between the two regions, usually calculated by meansof a first-order process. We will use here the subscript
f to represent fractures, inter-aggregate pores, orin which a is a kinematic exponent and b a macropore
conductance parameter. One advantage of thisapproach is that no water retention properties of themacropore region are required, and hence that the
number of parameters can be reduced. Such data areexperimentally very difficult to determine. One
disadvantage is that the kinematic wave equation islimited to vertical gravity-driven flow, since capillarity is ignored. Thus, application to two-dimensional
flow problems (i.e. lateral or downslope flows) facessome numerical and conceptual problems. The
kinematic wave approach also cannot explain the
non-equilibrium upward infiltration reported byŠimůnek et al. (2001). Contrary to kinematic wave
models, dual-porosity models based on the Richards
equation do account for both gravitational flow andcapillary-driven flow. Hence, Eq. (9) can describe not

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35
Fig. 3. Concentration profiles (top) and breakthrough curves at 10, 20, and 30 cm (bottom) calculated assuming (a) equilibrium water flow (Eq. (1)) and solute transport (Eq. (2)), (b)non-equilibrium water flow (Eq. (7)) (Ross and Smettem, 2000) and equilibrium solute transport (Eq. (2)), and (c) both non-equilibrium water flow (Eq. (7)) and non-equilibrium
solute transport (Eq. (13)).

21

22

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

only vertical downward flow during sudden rainfall
events, but also upward flow during evaporation
periods.
The dual-porosity formulation for solute transport
is similarly based on the convection – dispersion andmass balance equations as follows
›uf c f›f rs f
›
›c
›qcf
þ
¼
2 f f 2 Gs ;uf D f f 2
›z
›t
›t
›z
›z
›um c m
›ð1 2 f Þrsm
þ
ð13Þ
›t
›t
¼ 2fm þ Gs
for the fractures (macropores) and matrix, respectively, where f is the dimensionless fraction of sorption
sites in contact with the fractures (mobile water), and
Gs is the solute transfer rate between the two regions.
Although dual-porosity models have been popularly
used for solute transport studies (van Genuchten,1981), their application to water flow problems has
been attempted far less frequently.
As an example, Fig. 4 shows computed water andsolute distributions during infiltration obtained with adual-porosity model. The soil hydraulic parameters ofthe macropore domain were taken as follows: ur ¼
0:0; us ¼ 0:200; a1 ¼ 0:041 cm21 ; n ¼ 1:964; l ¼
0:5; Ks ¼ 0:000722 cm s21 ; while the immobile
matrix domain was assumed to have a saturated
water content, usm ; of 0.15. Initial conditions were set
equal to the pressure head of 2 150 cm. We assumed
that water mass transfer was described with Eq. (20),in which the mass transfer constant v was set at
0.00001 s21. For simplicity, we considered only
convective solute mass transfer between the twopore regions (i.e. no diffusive transfer), with thedispersivity again fixed at 2 cm. While for ponded
surface conditions water in the fracture domain
quickly reached full saturation (Fig. 4(a)), the watercontent of the matrix increased only gradually withtime. Consequently, the total water content, defined as
the sum of water contents in both the fracture andmatrix domains, also increased only gradually. The
total water content would be the quantity measured
with most field water content measurement devices,such as TDR or neutron probe. Pressure headmeasurements using tensiometers are on the other
hand, often dominated by the wetter fracture domainthat reaches equilibrium relatively quickly. Similarly

to the model of Ross and Smettem (2000), as well asother more complex non-equilibrium models, thedual-porosity model can therefore explain often
observed non-equilibrium between pressure heads
and water contents (Šimůnek et al., 1999b, 2001,
among others). Similar non-equilibrium profiles as forthe water content were also obtained for the soluteconcentration (Fig. 4(b)).
4.3. Dual-permeability models
Different types of dual-permeability approaches
may be used to describe flow and transport instructured media. Several assume similar governing
equations to describe flow in the fracture and matrix
regions, while others use different formulations for thetwo regions. A typical example of the first approach is
the work of Gerke and van Genuchten (1993a, 1996)
who applied Richards equations to each of two poreregions. The flow equations for the fracture (subscriptf) and matrix (subscript m) pore systems are,respectively,


›uf
›
›hf
G
¼
K
þ K f 2 Sf 2 w›z f ›z
›t
w

ð14Þ

and


›um
›
›hm
Gw
¼
K
þ Km 2 Sm þ›z m ›z
›t
12w

ð15Þ

where w is the ratio of the volumes of the fracture(inter-aggregate) and the total pore systems, ufs =us :This approach is relatively complicated in thatthe model requires characterization of water retention and hydraulic conductivity functions (potentially of different form) for both pore regions,as well as the hydraulic conductivity function of
the fracture – matrix interface (to be discussed
later). Note that the water contents uf and um in
Eqs. (14) and (15) have different meanings than in
Eq. (7) where they represented absolute value ofthe water content in the total pore-space, while
here they denote water contents in the poresubspace (fracture or matrix).
A similar expression to Eq. (11) for the volumetricflux density in macropores, q, was suggested by Jarvis(1994) for use in the dual-permeability model

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35
Fig. 4. Water content (a) and concentration (b) profiles in the fracture domain, matrix domain, and both domains combined, as well as water (c) (Eq. (20)) and solute (d) (Eq. (27))
mass transfer terms calculated with the dual-porosity model.

23

24

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

MACRO:

uf aq¼b¼ bðSfe Þa
usf

ð16Þ

where usf is the saturated water content of themacropores (i.e. macroporosity). Since the macropore
water content in Eq. (16) is scaled by usf ; the fracture
conductance parameter, b, has the more direct
physical meaning of the hydraulic conductivity ofmacropores at full saturation. We note that the
MACRO model assumes applicability of the Richards
equation to flow in the soil matrix.
An alternative approach for flow in the matrix andfractures was used by Ahuja and Hebson (1992) in the
RZWQM model. They described infiltration andredistribution in the soil matrix with the Green –Ampt and Richards equations, respectively, andinfiltration into cylindrical macropores using Poiseuille’s law assuming gravitational flow:
q¼

Np rgprp48n

ð17Þ

where Np is the number of pores per unit area, r is thedensity of water, g is the gravitational constant, rp isthe radius of the macropores, and n the dynamic
viscosity of water. Ahuja and Hebson (1992) used asimilar Poiseuille-type expression for planar fractures.
One disadvantage of using pore-scale physical laws isthat macropores probably rarely flow at full saturation, while their shape, form and continuity will notnormally correspond to the assumed simplified
geometry.
Analogous to Eqs. (14) and (15), the dualpermeability formulation for solute transport can bebased on convection– dispersion type equations fortransport in both the fracture and matrix regions asfollows (Gerke and van Genuchten, 1993a):›uf c f›f rs f
þ
›t
›t


›
›c
›qcf
G
¼
2 ff 2 s
uf D f f 2
›z
›z
›z
w

Eq. (18) assumes a complete convective – dispersivetype transport description for the fractures. Several
authors simplified transport in the macropore domain,for example, by considering only piston displacementof solutes (Ahuja and Hebson, 1992; Jarvis, 1994).
Fig. 5 shows an application of the dual-permeability model of Gerke and van Genuchten(1993a) to infiltration into a macroporous soil profile.
We used the following soil hydraulic parameters forthe matrix (ur ¼ 0:105; us ¼ 0:50; a1 ¼ 0:005 cm21 ;n ¼ 1:5; l ¼ 0:5; Ks ¼ 1:05 cm d21 ), the fracture(ur ¼ 0:0; us ¼ 0:50; a1 ¼ 0:1 cm21 ; n ¼ 2:0; l ¼
0:5; Ks ¼ 2000:0 cm d21 ), and the matrix – fractureinterface (w ¼ 0:05; b ¼ 3; g ¼ 0:4; a ¼ 1:0; andKas ¼ 0:01 cm d21 ; see Eq. (24)) domains. We
assumed that water mass transfer was described with
Eq. (23). Similarly as for the dual-porosity example,
we again considered only convective solute masstransfer between the two pore regions, with adispersivity of 2 cm. The initial pressure head forboth regions was equal to 2 100 cm, while water wasapplied only to the soil matrix until ponding occurred.
The soil matrix at the surface quickly became
saturated since the applied water flux (50.0 cm d21)
greatly exceeded the saturated hydraulic conductivityof the matrix ðKs ¼ 1:05 cm d21 Þ (Fig. 5(a)). Water
started to flow into the fracture domain (Fig. 5(a))
where it moved vertically downwards (Fig. 5(c)) andinfiltrated into the matrix domain (Fig. 5(d)). The two
distinct moisture fronts in Fig. 5(b) (and concentration
fronts in Fig. 5(e)) are due to lateral transfer from thefracture domain (the leading edge), and flow in thematrix domain (the second front). Again, in order todemonstrate sensitivity of the model predictions to themass transfer coefficient (see Eq. (24)), we included
into Fig. 5 final profiles ðt ¼ 0:08 dÞ calculated when
the conductivity of the matrix– fracture interface wasdecreased or increased by 25% (i.e. Kas ¼ 0:0075 and0.0125 cm d21, respectively).
4.4. Fine-textured soils

ð18Þ

›um c m
›ð1 2 f Þrsm
þ
›t
›t


›
›c
›qcm
Gs
¼
2 fm 2
ð19Þ
um D m m 2
›z
›z
›z
12w

Another class of models has been developed forfine-textured swelling clay soils containing shrinkage
cracks. Models of this type may take into accountshrinking and swelling as a function of the watercontent, and use statistical properties of cracks (e.g.
depth and width) as parameter input. Thus, they do not

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35
Fig. 5. Infiltration and mass exchange fluxes (a), water contents in the matrix (b) and fracture (c) domains, water mass exchange rates (d) and concentrations in the matrix (e) andfracture (f) domains calculated with the dual-permeability model of Gerke and van Genuchten (1993a,b). Thin lines represent final profiles ðt ¼ 0:08 dÞ when the conductivity of theinterface is decreased or increased by 25%.

25

26

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

require detailed knowledge of the crack system and/or
the spatial distribution of macropores (Slawinski et al.,
1996; van Dam et al., 1997). These models typicallyassume that water and solutes can move instantaneously to specified bypass depths once the infiltration capacity of the soil matrix is exceeded by therainfall rate and a critical depth of water has formed at
the soil surface (Verburg et al., 1996; van Dam et al.,
1997; Novák et al., 2000). Downward bypassing or
short-circuiting will not occur if the matrix absorbs allwater through the soil surface. Since in all of thesemodels water flow in the matrix is described using the
Richards equation, while water in the cracks moves
from the soil surface to some specified depth (albeit
instantaneously), these models can be viewed as a
subgroup of the dual-permeability models. Water
infiltrating from the cracks into the soil matrix isestimated using either Darcy’s law (van Dam et al.,1997) or the Green – Ampt approach (Novák et al.,2000) and then added as a source term to the Richards
equation describing flow in the soil matrix.
Fig. 6 compares various fluxes for two infiltration
scenarios, one with and one without the presence ofsoil drying cracks as calculated with the fracture
program of Novák et al. (2000). Because the applied
water flux (25 cm d21 for 2 h) at the soil surface(curve 1) is much higher than the saturated hydraulicconductivity of the matrix ðKs ¼ 5 cm d21 Þ; the soil
surface ponds quickly, a critical layer of water forms
on the soil surface, and water either further accumulates at the soil surface or starts to run off when no
cracks are present. Alternatively, when cracks arepresent, water starts to flow into the cracks (curve 5),
from which it infiltrates into the soil (curve 4). As
expected, when water is allowed to accumulate on thesoil surface (no cracks present; curve 2), it takes much
longer for this water to infiltrate than when it isallowed to flow into cracks (curve 3). If water is notallowed to accumulate on the soil surface and no
cracks are present (curve 3), then most of this water islost by surface run off (curve 5). When cracks arepresent, cumulative infiltration (curve 6) proceeds
much faster than when water is allowed to accumulate
at the soil surface (curve 2), and is much larger thanwhen water is allowed to run off (curve 3). We refer toNovák et al. (2000) for a more detailed discussion ofthis example.
4.5. Multiple-porosity/permeability models
Multiple-porosity/permeability models are conceptually similar to dual-porosity/permeability models,
except that they implement additional overlapping
pore regions. This allows for greater flexibility, albeit

Fig. 6. Actual (a) and cumulative (b) fluxes versus time for infiltration of water into a soil with and without cracks and/or macropores (afterNovák et al. (2000)): (1) irrigation, (2) infiltration into a soil without cracks with water accumulating on the soil surface, (3) infiltration throughthe soil surface with surface runoff or with cracks present, (4) infiltration from the cracks into the soil matrix, (5) flow into cracks or surface runoff, and (6) total infiltration into the soil matrix when cracks are present.

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

at the expense of requiring more parameters that may
also be physically poorly defined. For example, Gwo
et al. (1995) developed the MURF and MURT models
for multi-region flow and transport, respectively,
assuming three overlapping pore regions: macropores,
mesopores, and micropores (or using alternative
terms: primary fractures, secondary fractures andsoil matrix). Although different functions for the soil
hydraulic properties of each region were used, they
assumed applicability of the Richards and convection –dispersion equations in all three regions. Hutson
and Wagenet (1995) developed the TRANSMIT
model that considers n overlapping regions, withflow and transport in each region again describedusing the Richards and convection – dispersion
equations. Similar to Gwo et al. (1995), they allowed
for water and solute to exchange between all regions.
4.6. Functional approaches
A large number of functional models, generally
capacity type approaches, has also been developed toaccount for preferential water flow and/or solutetransport (Addiscott, 1977; Addiscott et al., 1986;
Corwin et al., 1991; Emermann, 1995). These
approaches are usually motivated by a desire towork with relatively simple preferential models that
require only a few parameters. One disadvantage
inherent to these functional models is that the rate ofmacropore flow may depend implicitly on the adopted
time step in the model. Emermann (1995) showed thata simple tipping bucket model for macropore flow
could be reformulated in terms of a differential
equation, thereby avoiding this limitation. The resulting model resembles the kinematic wave equation
discussed earlier, with the difference that Emermann(1995) assumed that macropore flow was a linear
function of the macropore water content.

5. Mass transfer
Critical in all dual-porosity/permeability models
are the coupling terms Gw and Gs describing masstransfer between the fracture and matrix regions.
Some descriptions are physically based, thoughapproximate, while others are entirely empirical.

27

Below we list some of the formulations most often
used.
5.1. Mass transfer driven by fluid saturation
The mass transfer rate, Gw ; for water between thefracture and matrix regions in several dual-porosity
studies (Phillip, 1968; Šimůnek et al., 2001) has been
assumed to be proportional to the difference in
effective water contents of the two regions using thefirst-order rate equation:

Gw ¼

›um
¼ v½Sfe 2 Sm
e 
›t

ð20Þ

where um is the matrix water content, v is a first-order
rate coefficient (T21), and Sfe and Sm
e are effective fluid
saturations of the fracture and matrix regions,
respectively. Compared to assuming a pressure
head-based driving force (to be discussed later), thedual-porosity model based on this mass transfer
equation requires significantly fewer parameters
since one does not need to know the retention function
for the matrix region explicitly, but only its residual
and saturated water contents. Coupling (Eq. (20)) witha dual-porosity non-equilibrium flow model leads tothe usual soil hydraulic parameters needed for theequilibrium model, two additional parameters characterizing the matrix region (i.e. its residual, um
r ; andsaturated, um;
water
contents),
andthefirst-ordermassstransfer coefficient v. By additionally assuming thatthe residual water content of the fracture region isequal to zero (and hence that residual water is present
only in the immobile region), one could further
decrease the number of model parameters. The
resulting model was used for the calculations shownin Fig. 3.
Eq. (20) assumes that the mass transfer rate is
proportional to the difference in effective watercontents, rather than pressure heads (Gerke and van
Genuchten, 1993b), which should provide a more
realistic description of the exchange rate between
the fracture and matrix regions. Eq. (20) thus
inherently assumes that the water retention propertiesof the matrix and the fracture domains are identical.
For this reason, Eq. (20) must be used with somecaution and probably only for dual-porosity models.
Dual-porosity models that use Eq. (20) are in somerespects conceptually similar to the non-equilibrium

28

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

model of Ross and Smettem (2000) discussed earlier.
This is because both models assume that the approachto equilibrium (as characterized by either the masstransfer coefficient v or the equilibration time
constant t ) is independent of the actual saturationlevel, only to its absolute deviation from equilibrium.
In the absence of gravity, the Richards equation
can be recast for a homogeneous medium as adiffusion equation with gradients in water content as
the driving force. The dual-permeability model
MACRO (Jarvis, 1994) uses a formulation for themass transfer term that is based on a first-orderapproximation to the water diffusion equation:›ubDw gw
Gw ¼ m ¼
ð21Þð ub 2 um Þ›td2
where ub is the saturated water content in micropores,
d is an effective ‘diffusion’ pathlength (i.e. half theaggregate width or half the fracture spacing), b is ashape factor that depends on the geometry, and gw is ascaling factor ( ¼ 0.4) obtained by matching theresults of the first-order approach at the half-time levelof the cumulative infiltration curve to the numericalsolution of the horizontal infiltration equation (Gerkeand van Genuchten, 1993b). The effective water
diffusivity Dw given by
Dub þ Dum f
Dw ¼
ð22Þ
Se
2
where Dub and Dum are soil water diffusivities atsaturation (saturation is assumed to exist at thefracture – matrix interface, ub ) and in the matrixðum Þ; respectively, and where the effective saturationof the macropores, Sfe ; is introduced to account for an
incomplete wetted contact area between the two poredomains. The value of b varies with the geometry
assumed for the aggregates (van Genuchten and
Dalton, 1986) and in MACRO this parameter is set
to 3 for rectangular slabs. The model based on Eqs.
(21) and (22) is clearly more physically based thanEq. (20) since the water content driving force pertains
to the same material (i.e. the matrix), being given by
the difference between the saturated water content of
the micropores (i.e. the water content at the interface
between micropores and macropores when macropores are active) and the average water content of the
micropores. It is therefore assumed that during active
water exchange between micropores and macropores,

the water pressures at the interface instantaneously
reach equilibrium at the water-entry pressure of themacropores. It can be noted that this is not the same asassuming that the matrix has to be saturated for bypass
flow of water to occur in the macropores. Eq. (21)only describes flow from macropores to matrix and
not in the reverse direction. It is intended for use inmodels based on cut and join hydraulic functions (i.e.
similar to Eq. (5)) and could not be used without
modification in dual-permeability models based on
the concept of overlapping pore continua (i.e. Eqs. (3)and (4)). In MACRO, the boundary between macropores and matrix is defined by the saturated matrix
water content (defined earlier), the equivalent pressure head ðhb Þ and the hydraulic conductivity at ub andhb (Jarvis, 1994). Flow from matrix to macropores
occurs instantaneously if the matrix water content
exceeds ub (i.e. the matrix pressure head exceeds hb ),following the basic physical principle that governs thefilling of pores when the water-entry pressure isexceeded. This can occur internally within the soil,
and also at the soil surface when rainfall intensity
exceeds the matrix infiltration capacity.
5.2. Mass transfer driven by pressure head gradient
The rate of exchange of water between the fracture
and matrix regions, Gw ; can also be assumed to be
proportional to the difference in pressure heads
between the two pore regions (Gerke and van
Genuchten, 1993a):

Gw ¼ aw ðhf 2 hm Þð23Þ

in which aw is a first-order mass transfer coefficient.
Since pressure heads are now needed for both regions,this approach requires estimating retention curves forboth pore regions. For porous media with well-definedgeometries, the first-order mass transfer coefficient,aw ; can be defined as follows (Gerke and van
Genuchten, 1993b):

aw ¼

b
Ka gw
d2

ð24Þ

Gerke and van Genuchten (1996) evaluated theeffective hydraulic conductivity Ka of the fracture –matrix interface using a simple arithmetic average

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

involving both hf and hm as follows
Ka ðhÞ ¼ 0:5½Ka ðhf Þ þ Ka ðhm Þð25Þ

The use of Eq. (24) implies that the medium contains
geometrically well-defined rectangular or other typesof macropores or fractures (Edwards et al., 1979; van
Genuchten and Dalton, 1986). While geometrically
based models are conceptually attractive, they may be
too difficult to use for field applications, partly
because structured soils and rocks usually containmixtures of aggregates and matrix blocks of various
sizes and shapes, but also because the parameters in
Eq. (24) may not be identifiable. Hence, rather thanusing Eq. (24) directly, one could also lump b, d, andgw into one effective hydraulic conductivity Kap of thefracture –matrix interface to giveaw ¼ Kap ðhÞ

ð26Þ
Kap

in which case
can then be used as calibration
parameter. To overcome the requirement for geometrically well-defined macropores or fractures in
Eq. (24), Gerke and van Genuchten (1996) suggested
a method for deriving values of the shape factor b in
Eq. (24) for non-idealized geometries or mixtures ofshapes in order to obtain a more general ‘macroscopic’ approach.
Because one needs to define retention curves forboth regions, mass transfer based on pressure headdifferences is more complex than mass transfer basedon a water content difference. Nevertheless, it may be
preferred to Eq. (20) since the difference in pressure
heads is the actual driving force for mass transfer. On
the other hand, the mass transfer term given by
Eq. (23) may be inherently more unstable numericallysince it involves calculating the product of two highly
non-linear terms, each of which can become extremely small or large. In contrast, Eq. (21) is morestable since the diffusivity term is less non-linear thanK, while the water content differences never become
extreme. Mass transfer based on Eqs. (23) and (24)was used for the calculations presented in Fig. 5.
In order to demonstrate the difficulty in distinguishing between various dual-porosity/permeability models and between various descriptions of masstransfer, which is practically very difficult, if not
impossible, to measure, we repeated the calculations
presented in Fig. 4 using different definitions for themass transfer term with all other parameters, initial29

and boundary conditions the same as before. Fig. 7
shows distributions versus depth of the water mass
exchange rate, Gw ; calculated with the dual-porosity
model when using the saturation driven exchange
term (Eq. (20)) and the pressure head driven exchange
term (Eq. (23)) assuming a constant (independent ofh ) and variable (dependent on h ) parameter aw(Eq. (24)). The soil hydraulic parameters for thematrix, for the pressure head driven case, were taken
to be the same as for the fracture domain. The firstorder rate coefficients v and aw in the saturation andpressure head driven (assuming a constant parameteraw ) mass transfer terms were set equal to1.0 £ 1025 s21 and 1.0 £ 1027 cm21 s21, respectp
in the pressure head driven
ively. The parameter Kas
mass transfer with a variable parameter aw was set
equal to 4.0 £ 1027 cm21 s21. Although radicallydifferent descriptions of the mass transfer term wereused in these calculations, all three runs resulted into
visually almost identical pressure head and watercontent profiles in the fracture domain, while only
very small differences were noticeable in the matrix
domain (Fig. 4). Fig. 7 reflects different levels of nonlinearity in the various definitions of the mass transfer
term. While effective saturation can range between 0
and 1, pressure heads in these calculations ranged
between 2 150 cm and 0. This is the main reason why
the first-order rate coefficients in the saturation andpressure head driven mass transfer terms vary by afactor of up to two.
5.3. Other approaches to mass transfer
Several other approaches for calculating the masstransfer term have also been used. For example, Ahuja
and Hebson (1992) and Novák et al. (2000) modeled
mass transfer using the Green and Ampt equation,
Chen and Wagenet (1992) used the Philip infiltration
model, while Zimmerman et al. (1993) used the
Warren – Root model, as well as developed their own
formulation.
5.4. Solute mass transfer
The transfer rate, Gs ; for solutes between thefracture and matrix regions is usually given as the sumof diffusive and convective fluxes, and can be written

30

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

as (Gerke and van Genuchten, 1996):

Gs ¼ as ð1 2 wm Þðcf 2 cm Þ þ Gw cp

ð27Þ

where cp is equal to cf for Gw . 0 and cm for Gw , 0;and as is the first-order solute mass transfer coefficient(T21) of the form

as ¼

b
Da
d2

ð28Þ

in which Da is an effective diffusion coefficient(L2T21) which represents the diffusion properties ofthe fracture – matrix interface, as well as otherparameters.

6. Model complexity and parameterization

Fig. 7. The water mass exchange rate profiles calculated with thedual-porosity model with (a) the saturation driven exchange term
(Eq. (20)), and the pressure head driven exchange term (Eq. (23))
with (b) a constant and (c) variable parameter aw (Eq. (24)).

The main disadvantage of dual-porosity and dualpermeability models is that, contrary to models for asingle pore region, they require many more input
parameters to characterize both pore systems. For
example, the dual-permeability model of Gerke andvan Genuchten (1993a), in its full complexity, needs
16 parameters to describe water flow. This is because
the hydraulic properties of the fracture and matrix
regions are each characterized by six soil hydraulicparameters (the same number as the equilibrium
system). Additionally, estimates of structural parp
ameter w, as well as Kas; aa ; and na are needed tocharacterize the hydraulic conductivity of thefracture – matrix interface when using the Mualem –van Genuchten model. A previous assumption(Gerke and van Genuchten, 1993a, 1996) is that therelative hydraulic conductivity functions of theinterface and the matrix regions are the same, thus
requiring only the saturated interface conductivitypKas
to scale the interface relative hydraulic conductivity function. If the pore connectivity parameters, l,
are equal to 0.5 (Mualem, 1976) and, as before, theresidual water content of the mobile region isassumed to be zero, the number of parametersmdecreases to eleven: w, ufs ; umr ; u s ; af ; am ; nf ; nm ;fmp
Ks ; Ks and Kas :The number of input parameters decreases slightly
when the kinematic wave approach is used to describe
flow in the fracture domain. In this case, the matrix
domain is again characterized by six parameters(assuming the Mualem – van Genuchten model), while

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

the fracture domain and the inter-domain interface arecharacterized by the macroporosity, usf ; the fracture
conductivity, b, the exponent a, and an effective
diffusion path length d. Thus, this model needs a totalof 10 parameters to describe the porous medium.
A similar number of parameters (i.e. 11) is needed
for the dual-porosity model based on the Richards
equation, the van Genuchten – Mualem model for thesoil hydraulic properties, and the pressure head-basedmass transfer term. In this case, one needs six
parameters to characterize the fracture domain (ur ;us ; a, n, l, Ks ), four parameters for the matrix domainp(ur ; us ; a, n ), and one parameter for the interface ðKasÞ
(assuming that several other factors can be lumped
into this parameter). As discussed earlier, the numberof parameters decreases to 9 when the pressure headbased mass transfer term is replaced with the watercontent-based mass transfer term, since then only theparameters ur and us are needed to characterize thematrix.
The number of parameters of course increases
significantly for multi-porosity and multi-permeability models for which one needs to determine
or calibrate the soil hydraulic (retention and hydraulicconductivity) parameters for each transport domain, as
well as mass transfer coefficients between all transportdomains. In contrast, the non-equilibrium model ofRoss and Smettem (2000) requires only one additionalparameter as compared to equilibrium flow models.
Little guidance is available as to how to obtainthese parameters, either by direct measurement, apriori estimation, or some calibration technique(Beven, 1991; Clothier et al., 1995; Jaynes et al.,1995). Existing experiments rarely provide enoughinformation to fully calibrate non-equilibrium flow
models, let alone to distinguish between the various
non-equilibrium models. As illustrated in Fig. 7,
different non-equilibrium models can produce relatively similar results (see also Šimůnek et al. (2001)).
Experiments or devices hence must be designed thatprovide (either directly or by inverse modeling)
estimates of the many parameters needed in these
relatively complex models. Examples of new
approaches for analyzing water transfer properties ofaggregates with and without cutans based on tension
imbibition have been presented by Leeds-Harrison
et al. (1994) and Gerke and Köhne (2002). Disc
infiltrometers can greatly facilitate measurements of31

the composite (macropore plus matrix) hydraulicproperties. Disc infiltrometry methods involving
ponded and tension infiltrometers are now increasingly used for in situ measurements of the hydraulicconductivity. Advantages of these methods are thatthey can yield well-defined conductivity functions
near saturation, especially when combined with othermethods to measure the hydraulic properties atintermediate pressure head values in the mesopore/micropore range (Mohanty et al., 1997; Jarvis et al.,1999). These measurements may also be useful fordeveloping predictive (pedotransfer) functions formodel parameters in the absence of data. For example,
Smettem and Bristow (1999) showed that hydraulicconductivity measured by tension infiltrometer at
2 4 cm in 20 Australian soils was well explained byvariation in clay content. On a large dataset for 70 soil
horizons from 13 different countries, Jarvis et al.
(2002) found weaker but still significant relationships(r 2 values up to 49%) between textural properties andthe hydraulic conductivity measured at 2 10 cm bytension infiltrometer.
Very little is currently known about thepossibilities and potential problems of applying
inverse modeling techniques to preferential flow
models. Durner et al. (1999) showed that theparameters of bimodal water retention and hydraulic conductivity functions (Eqs. (3) and (4)) could
be determined by inverse modeling of multi-stepoutflow data. They claimed that the procedure wasrobust, leading to unique solutions with limited
data (water flow only), irrespective of the numberof parameters included, providing the underlying
model accurately represented the true soil hydraulicproperties. By comparison, Schwartz et al. (2000)
recently attempted to estimate the parameters of adual-permeability transport model by inverse modeling using steady-state solute breakthrough datafrom a variably charged tropical soil, in which thebromide ion may behave as a weakly sorbed
reactive solute. They encountered great difficulties
in obtaining physically realistic estimates of two
critical parameters: the dispersion coefficient in
the micropores (Dm in Eq. (19), which became far
too large), and the fraction of sorption sites in
the macropores ( f in Eqs. (13), (18), and (19)),
which tended towards unity (i.e. all sorption sites
presumably are located in the macropores).

32

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

They concluded that inverse procedures are problematic for dual-permeability transport models, evenfor the simple case of steady water flow with only
four unknown parameters. However, the data set
available to Schwartz et al. (2000) was not ideal inthat no conservative tracer was applied (bromide
was reactive in their soil), while only flux
measurements at the base of the column wereavailable, but no resident concentrations inside thecolumn. In contrast, Kätterer et al. (2001) used theinverse modeling package SUFI (Abbaspour et al.,1997) to successfully estimate the parameters of thedual-permeability model MACRO using column
breakthrough experiments for multiple non-reactive
tracers (deuterium, bromide, and chloride) and twosoil indigenous solutes (sulphate and nitrate). Theirstudy made use of both outflow flux concentration
data and resident concentrations inside the columns
at the end of the experiments.
7. Need for inter-code comparison
A relatively large number of models has beendeveloped during the last two decades that consider
non-equilibrium flow and transport. While several
studies have been carried out to compare various
codes describing uniform Darcian flow in thevadose zone (Gee et al., 1999; Ogan et al., 1999;
Scanlon et al., 2002), similar comparisons of codes
and/or approaches simulating preferential and/or
non-equilibrium flow are thus far lacking. This may
in large part be due to the fact that most of thepreferential flow codes are not readily available, arenot easy to use, lack numerical stability, and/or are
based on such significantly different descriptions of
the underlying flow processes that an effective
comparison is neither easy nor straightforward.
To be able to carry out such studies, wedeveloped new versions of the HYDRUS-1Dand HYDRUS-2D codes (Šimůnek et al., 1998,
1999a) that include a hierarchical system ofvarious approaches simulating preferential or
non-equilibrium flow and transport. Versions 2.0
of both codes could deal only with equilibrium
water flow (Eq. (1)), and with equilibrium or nonequilibrium mobile – immobile solute transport(Eqs. (2) and (13)). Recently updated versions

additionally include uniform flow with composite
retention and hydraulic conductivity functions(Eqs. (3) and (4)), the single-porosity non-equilibrium model of Ross and Smettem (2000) (Eq.
(7)), a dual-porosity model based on the Richards
equation (9) with either the water content or thepressure head-based mass transfer terms (Eqs. (20)and (23)), and finally two dual-permeability
models: one using the Richards equation for bothpore regions (Eqs. (14) and (15)), and one using
the kinematic wave approach for flow in themacropores (Eq. (16)). The different examples
discussed in this paper were all obtained with theupdated HYDRUS-1D code.
8. Concluding remarks
While macropore flow has important effects on
subsurface hydrology in general, and on infiltration
rates and unsaturated soil water distributions inparticular, its main implications are likely in theaccelerated movement of surface applied fertilizers or
pollutants through the vadose zone (Nielsen et al.,1986). At present it is still very difficult to use themore complex dual-permeability model involving two
coupled Richards equations to describe preferential
flow and transport under field conditions, partly
because of the large number of parameters involved,
and the current lack of standard experimental
techniques to obtain them. At present no examples
exist of such applications in the soil science andvadose zone hydrology literature. Hence, the use of
these models has so far been restricted to theoreticalapplications and laboratory studies carried out underwell-defined and controlled conditions. The dualpermeability model MACRO, based on the kinematicwave equation for flow in macropores, requires fewer
parameters, has been frequently applied to long-termtransient field experiments (Villholth and Jensen,
1998; Larsson and Jarvis, 1999a,b) and is also being
used for risk assessment for pesticide leaching within
the EU (FOCUS, 2000). The main reasons for thewider use of this model are that it is physically based
and numerically robust, while the difficult problem ofparameterization is not insurmountable.
Accurate coupling of the fracture and matrix
domains still represents the greatest challenge in

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

terms of successfully describing non-equilibrium flow
and transport in the vadose zone. Matrix –macropore(or matrix – fracture) interfaces can have very different
properties than the bulk matrix due to the depositionof organic matter, various types of coatings, finetexture mineral particles, or various oxides andhydroxides on the aggregate exteriors or macropore
walls; these coatings can markedly reduce rates ofdiffusion and mass flow between macropores and thesoil matrix (Thoma et al., 1992).
With some exceptions (Larsson and Jarvis,
1999a,b), the limited availability of comprehensive
data sets has so far restricted the field validation ofpreferential flow models. A need still remains toapply currently available modeling tools to natural
systems where sufficient data are available to firstparameterize the models and then to assess theadequacy of model predictions (Šimůnek and deVos, 1999; Evans et al., 2001). This type of twostage calibration process (sometimes called historical validation) may be difficult to follow in practicedue to experimental difficulties related to inaccessibility and observation scale in fractured rocks or
because of the inherently dynamic macropore
system in soils being subject to physical (swell/shrink, freeze/thaw), biological (variations in soil
faunal activity) and man-made disturbances(ploughing and other tillage practices).
Still, while the problem of modeling preferential
flow obviously will remain a challenge, several usefulapproaches have recently become available to study
and model preferential/non-equilibrium flow in structured and/or fractured media. Continuous advances inboth numerical techniques and computer hardware are
now making it increasingly possible to carry out
comprehensive simulations of non-equilibrium flow
processes in the vadose zone. Such simulations,especially if paired with good field data sets, arevital to better understanding and quantifying theoften-controlling effects of heterogeneities, fracturesand macropores on flow and transport at the field scale(van Genuchten et al., 1999).
Acknowledgements
This project is supported by NSF grant, NSFSAHRA at the University of Arizona and Swedish
Research Council.

33

References
Abbaspour, K.C., van Genuchten, M.T., Schulin, R., Schläppi, E.,
1997. A sequential uncertainty domain inverse procedure for
estimating subsurface flow and transport parameters. Water
Resour. Res. 33, 1879–1892.
Addiscott, T.M., 1977. A simple computer model for leaching instructured soils. J. Soil Sci. 28, 554–563.
Addiscott, T.M., Phillippa, J.H., Whitmore, A.P., 1986. Applicationof simple leaching models in heterogeneous soils. Geoderma 38,
194 –285.
Ahuja, L.R., Hebson, C., 1992. Root Zone Water Quality Model.
GPSR Technical Report No. 2, USDA, ARS, Fort Collins, CO.
Beven, K., 1991. Modeling preferential flow: an uncertain future? .
In: Gish, T.J., Shirmohannadi, A. (Eds.), Preferential Flow,
American Society of Agricultural Engineers, St Joseph, MI, pp.
1 –11.
Chen, C., Wagenet, R.J., 1992. Simulation of water and chemicals
in macropore soils. Part 1. Representation of the equivalentmacropore influence and its effect on soil water flow. J. Hydrol.
130, 105–126.
Clothier, B.E., Heng, L., Magesan, G.N., Vogeler, I., 1995. The
measured mobile-water content of an unsaturated soil as afunction of hydraulic regime. Aust. J. Soil Res. 33, 397–414.
Corwin, D.L., Waggoner, B.L., Rhoades, J.D., 1991. A functionalmodel of solute transport that accounts for bypass. J. Environ.
Qual. 20 (3), 647–658.
van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabat, P.,
van Valsum, P.E.V., Groenendijk, P., van Diepen, C.A., 1997.
Theory of SWAP, Version 2.0. Simulation of Water Flow,
Solute Transport and Plant Growth in the Soil – Water –Atmosphere-Plant Environment. Department of Water
Resources, WAU, Report 71, DLO Winand Staring Centre,
Wageningen, Technical Document 45.
Durner, W., 1994. Hydraulic conductivity estimation for soils
with heterogeneous pore structure. Water Resour. Res. 30,
211 –233.
Durner, W., Priesack, E., Vogel, H.-J., Zurmühl, T., 1999.
Determination of parameters for flexible hydraulic functions
by inverse modeling. In: van Genuchten, M.Th., Leij, F.J., Wu,
L. (Eds.), Characterization and Measurement of the Hydraulic
Properties of Unsaturated Porous Media, University of California, Riverside, CA, pp. 817 –829.
Edwards, W.M., van der Ploeg, R.R., Ehlers, W., 1979. A numerical
study of the effects of noncapillary-sized pores upon infiltration.
Soil Sci. Soc. Am. J. 43, 851 –856.
Emermann, S.H., 1995. The tipping bucket equations as a model for
macropore flow. J. Hydrol. 171, 23–47.
Evans, D.D., Rasmussen, T.C., Nicholson, T.J., 2001. Flow andtransport through unsaturated fractured rock: an overview. In:
Evans, D.D., Nicholson, T.J., Rasmussen, T.C. (Eds.), Flow andTransport Through Unsaturated Fractured Rock, 2nd ed,
Geophysical Monograph, vol. 42. AGU, Washington, DC, pp.
1 –18.
Flint, A.L., Flint, L.E., Bodvarsson, G.S., Kwicklis, E.M., FabrykaMartin, J., 2001. Development of the conceptual model of
unsaturated zone hydrology at Yucca Mountain, Nevada,34

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35

Conceptual Models of Flow and Transport in the Fractured
Vadose Zone, National Research Council, National Academy
Press, Washington, DC, pp. 47–85.
Flühler, H., Durner, W., Flury, M., 1996. Lateral solute mixing
processes—a key for understanding field-scale transport ofwater and solutes. Geoderma 70, 165–183.
FOCUS, 2000. FOCUS groundwater scenarios in the EU plant
protection product review process. Report of the FOCUS
Groundwater Scenarios Workgroup, EC Document Reference
Sanco/321/2000, DG SANCO, EU Commission, Brussels, p.
197.
Gee, G.W., Albright, W.H., Wilson, G.V., Fayer, M.J., Ogan, B.,
1999. Alternative Cover Assessment Project: Phase 1 Report,
Desert Research Institute, Las Vegas, NE.
van Genuchten, M.Th., 1980. A closedform equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc.
Am. J. 44, 892–898.
van Genuchten, M.Th., 1981. Non-equilibrium transport parametersfrom miscible displacement experiments. Research Report 119,
US Salinity Laboratory, USDA, ARS, Riverside, CA.
van Genuchten, M.Th., Dalton, F.N., 1986. Models for simulating
salt movement in aggregated field soils. Geoderma 38,
165–183.
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–481.
van Genuchten, M.Th., Schaap, M.G., Mohanty, B.P., Šimůnek,
J., Leij, F.J., 1999. Modeling transport processes at the local
scale. In: Feyen, J., Wiyo, K. (Eds.), Modelling of Transport
Process in Soils at Various Scales in Time and Space,
Proceedings of the International Workshop of EurAgEng’s
Field of Interest on Soil and Water, Leuven, Belgium, 24– 26
November, Wageningen Pers, Wageningen, The Netherlands,
pp. 23 –45.
Gerke, H.H., van Genuchten, M.Th., 1993a. A dual-porosity
model for simulating the preferential movement of water andsolutes in structured porous media. Water Resour. Res. 29,
305–319.
Gerke, H.H., van Genuchten, M.Th., 1993b. Evaluation of a firstorder water transfer term for variably saturated dual-porosity
flow models. Water Resour. Res. 29, 1225–1238.
Gerke, H.H., van Genuchten, M.Th., 1996. Macroscopic representation of structural geometry for simulating water and solute
movement in dual-porosity media. Adv. Water Resour. 19,
343–357.
Gerke, H.H., Köhne, J.M., 2002. Estimating hydraulic properties ofsoil aggregate skins from sorptivity and water retention. Soil
Sci. Soc. Am. J. 66 (1), 26–36.
Germann, P.F., 1985. Kinematic wave approach to infiltration anddrainage into and from soil macropores. Trans. ASAE 28,
745–749.
Germann, P.F., Beven, K., 1985. Kinematic wave approximation to
infiltration into soils with sorbing macropores. Water Resour.
Res. 21 (7), 990 –996.
Gwo, J.P., Jardine, P.M., Wilson, G.V., Yeh, G.T., 1995. A
multiple-pore-region concept to modeling mass transfer insubsurface media. J. Hydrol. 164, 217– 237.

Hendrickx, J.M.H., Flury, M., 2001. Uniform and preferential flow,mechanisms in the vadose zone, Conceptual Models of Flowand Transport in the Fractured Vadose Zone, National ResearchCouncil, National Academy Press, Washington, DC, pp. 149–
187.
Hendrickx, J.M.H., Dekker, L.W., Boersma, O.H., 1993. Unstable
wetting fronts in water repellent field soils. J. Environ. Qual. 22,
109 –118.
Hutson, J.L., Wagenet, R.J., 1995. A multiregion model describing
water flow and solute transport in heterogeneous soils. Soil Sci.
Soc. Am. J. 59, 743–751.
Jarvis, N.J., 1994. The MACRO Model (Version 3.1). Technical
Description and Sample Simulations. Reports and Dissertations
19. Department of Soil Science, Swedish University ofAgricultural Science, Uppsala, Sweden, p. 51.
Jarvis, N.J., 1998. Modeling the impact of preferential flow on
nonpoint source pollution. In: Selim, H.M., Ma, L. (Eds.),
Physical Nonequilibrium in Soils: Modeling and Application,
Ann Arbor Press, Chelsea, MI, pp. 195 –221.
Jarvis, N.J., Messing, I., Larsson, M.H., Zavattaro, L., 1999.
Measurement and prediction of near-saturated hydraulicconductivity for use in dual-porosity models. In: van Genuchten,
M.Th., Leij, F.J., Wu, L. (Eds.), Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media,
University of California, Riverside, CA, pp. 783 –798.
Jarvis, N.J., Zavattaro, L., Rajkai, K., Reynolds, W.D., Olsen, P.-A.,
McGechan, M., Mecke, M., Mohanty, B., Leeds-Harrison, P.B.,
Jacques, D., 2002. Indirect estimation of near-saturated
hydraulic conductivity from readily available soil information.
Geoderma in press.
Jaynes, D.B., Logsdon, S.D., Horton, R., 1995. Field method formeasuring mobile/immobile water content and solute transfer
rate coefficient. Soil Sci. Soc. Am. J. 59, 352 –356.
Kätterer, T., Schmied, B., Abbaspour, K.C., Schulin, R., 2001.
Single and dual-porosity modelling of multiple tracer transport
through soil columns: effects of initial moisture and mode ofapplication. Eur. J. Soil Sci. 52, 1–12.
Larsson, M.H., Jarvis, N.J., 1999a. Evaluation of a dual-porosity
model to predict field-scale solute transport in a macroporous
soil. J. Hydrol. 215, 153–171.
Larsson, M.H., Jarvis, N.J., 1999b. A dual-porosity model toquantify macropore flow effects on nitrate leaching. J. Environ.
Qual. 28, 1298– 1307.
Leeds-Harrison, P.B., Youngs, E.G., Uddin, B., 1994. A device fordetermining the sorptivity of soil aggregates. Eur. J. Soil Sci. 45,
269 –272.
Leij, F.J., Russell, W.R., Lesch, S.M., 1997. Closed-form
expressions for water retention and conductivity data. Ground
Water 35 (5), 848 –858.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., van Genuchten,
M.Th., 1997. New piecewise-continuous hydraulic functions formodeling preferential flow in an intermittent flood-irrigatedfield. Water Resour. Res. 33, 2049–2063.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., Šimůnek, J., vanGenuchten, M.Th., 1998. Preferential transport of nitrate to atile drain in an intermittent flood-irrigated field: model

J. Šimůnek et al. / Journal of Hydrology 272 (2003) 14–35
development and experimental evaluation. Water Resour. Res.
34, 1061–1076.
Mualem, Y., 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour. Res.
12, 513 –522.
National Research Council, 2001. Conceptual Models of Flow andTransport in the Fractured Vadose Zone, National Academy
Press, Washington, DC, pp. 9–45.
Nielsen, D.R., van Genuchten, M.Th., Biggar, J.W., 1986. Water
flow and solute transport processes in the unsaturated zone.
Water Resour. Res. 22 (9), 89S–108S.
Novák, V., Šimůnek, J., van Genuchten, M.Th., 2000. Infiltration ofwater into soils with cracks. ASCE J. Irrig. Drain. Engng 126(1), 41–47.
Ogan, B.D., Wilson, G.V., Albright, W., Gee, G.W., Fayer, M.J.,
Rock, S., 1999. Sensitivity analysis and validation of numericalmodels used in the design of alternative landfill covers. Annual
Meetings Abstracts, SSSA, Salt Lake City, October 31 –November 4, 1999, p. 36.
Othmer, H., Diekkrüger, B., Kutı́lek, M., 1991. Bimodal porosity
and unsaturated hydraulic conductivity. Soil Sci. 152 (3),139–150.
Peters, R.R., Klavetter, E.A., 1988. A continuum model for water
movement in an unsaturated fractured rock mass. Water Resour.
Res. 24 (3), 416 –430.
Phillip, J.R., 1968. The theory of absorption in aggregated media.
Aust. J. Soil Res. 6, 1–19.
Pruess, K., Wang, J.S.Y., 1987. Numerical modeling of isothermal
and non-isothermal flow in unsaturated fractured rock—areview. In: Evans, D.D., Nicholson, T.J. (Eds.), Flow andTransport through Unsaturated Fractured Rock, Geophysics
Monograph, vol. 42. American Geophysical Union, Washington, DC, pp. 11–22.
Ross, P.J., Smettem, K.R., 2000. A simple treatment of physicalnonequilibrium water flow in soils. Soil Sci. Soc. Am. J. 64,
1926–1930.
Scanlon, B.R., Christman, M., Reedy, R., Porro, I., Šimůnek, J.,
Flerchinger, G., 2002. Intercode Comparisons for simulatingwater balance of surficial sediments in semiarid regions. Water
Resour Res., in press.
Schultze, B., Ippisch, O., Huwe, B., Durner, W., 1999. Dynamic
nonequilibrium during unsaturated water flow. In: van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Characterization andMeasurement of the Hydraulic Properties of Unsaturated PorousMedia, University of California, Riverside, CA, pp. 877 –892.
Schwartz, R.C., Juo, A.S.R., McInnes, K.J., 2000. Estimating
parameters for a dual-porosity model to describe non-equilibrium, reactive transport in a fine-textured soil. J. Hydrol. 229,
149–167.
Šimůnek, J., de Vos, J.A., 1999. Inverse optimization, calibrationand validation of simulation models at the field scale. In: Feyen,J., Wiy, K. (Eds.), Modelling of Transport Process in Soils at
Various Scales in Time and Space, Proceedings of theInternational Workshop of EurAgEng’s Field of Interest on
Soil and Water, Leuven, Belgium, 24–26 November, Wageningen Pers, Wageningen, The Netherlands, pp. 431 –445.

35

Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1998. The
HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably
saturated media. Version 2.0, IGWMC-TPS-70, International
Ground Water Modeling Center, Colorado School of Mines,
Golden, Colorado, p. 202.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1999a. The
HYDRUS-2D software package for simulating two-dimensional
movement of water, heat, and multiple solutes in variably
saturated media. Version 2.0, IGWMC-TPS-53, International
Ground Water Modeling Center, Colorado School of Mines,
Golden, Colorado, p. 251.
Šimůnek, J., Wendroth, O., van Genuchten, M.Th., 1999b.
Estimating unsaturated soil hydraulic properties from laboratory
tension disc infiltrometer experiments. Water Resour. Res. 35(10), 2965–2979.
Šimůnek, J., Wendroth, O., Wypler, N., van Genuchten, M.Th.,
2001. Nonequilibrium water flow characterized from an upward
infiltration experiment. Eur. J. Soil Sci. 52 (1), 13–24.
Skopp, J., 1981. Comment on ‘Micro- meso- and macroporosity ofsoil’. Soil Sci. Soc. Am. J. 45, 1246.
Slawinski, C., Sobczuk, H.A., Walczak, R.T., 1996. Submodel of
bypass flow in cracking soil. Part 1. Theory. Int. Agrophys. 10,
189 –196.
Smettem, K.R.J., Chettleborough, D.J., Fichards, B.G., Leaney,
F.W., 1991. The influence of macropores on runoff generationfrom a hillslope soil with a contrasting textural class. J. Hydrol.
122, 235–252.
Smettem, K.R.J., Bristow, K.L., 1999. Obtaining soil hydraulicproperties for water balance and leaching models from survey
data. 2. Hydraulic conductivity. Aust. J. Agric. Res. 50,
1259–1262.
Thoma, S.G., Gallegos, D.P., Smith, D.M., 1992. Impact of fracture
coatings on fracture/matrix flow interactions in unsaturated,porous media. Water Resour. Res. 28, 1357–1367.
Verburg, K., Ross, P.J., Bristow, K.L., 1996. SWIMv2.1 User
Manual. Divisional Report No. 130, CSIRO, Australia.
Villholth, K.G., Jensen, K.H., 1998. Flow and transport processes ina macroporous subsurface drained glacial till soil. II. Model
analysis. J. Hydrol. 207, 121 –135.
Vogel, T., Cı́slerová, M., 1988. On the reliability of unsaturated
hydraulic conductivity calculated from the moisture retention
curve. Transport Porous Media 3, 1–15.
de Vos, J.A., Šimůnek, J., Raats, P.A.C., 1999. Identification of thehydraulic properties of a layered silt loam. In: van Genuchten,
M.Th., Leij, F.J., Wu, L. (Eds.), Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media,
University of California, Riverside, CA, pp. 783 –798.
Wildenschild, D., Hopmans, J.W., Šimůnek, J., 2001. Flow ratedependence of soil hydraulic characteristics. Soil Sci. Soc. Am.
J. 65 (1), 35–48.
Zimmerman, R.W., Chen, G., Hadgu, T., Bodvarsson, G.S., 1993.
A numerical dual-porosity model with semianalytical treatment of fracture/matrix flow. Water Resour. Res. 29,
2127–2137.