Journal of Hydrology 238 (2000) 78–89
www.elsevier.com/locate/jhydrol

Modeling flow and transport in a two-dimensional
dual-permeability system with spatially variable
hydraulic properties
T. Vogel a,*, H.H. Gerke b,1, R. Zhang c,2, M.Th. Van Genuchten d,3baCzech Technical University, Faculty of Civil Engineering, Thakurova 7, 16629 Prague, Czech Republic
Department of Soil Landscape Research, Center for Agricultural Landscape and Land Use Research (ZALF), Eberswalder Str. 84,
D-15374 Müncheberg, GermanycUniversity of Wyoming, Department of Renewable Resources, Laramie, WY 82071-3354, USAdGeorge E. Brown Jr., Salinity Laboratory, USDA-ARS, 450 West Big Springs Road, Riverside, CA 92507-4617, USA

Received 22 December 1999; revised 6 June 2000; accepted 21 August 2000

Abstract
Most field soils exhibit soil spatial variability as well as soil structure. The challenge is to account adequately for both typesof spatial heterogeneity in simulation models. A numerical finite element code was used to compare single- and dual-permeability approaches for modeling variably saturated flow and transport in two-dimensional heterogeneous soil systems. The code
was based on the Richards’ equation for water flow and the advection-dispersion equation for solute transport. Spatial
variability in the soil hydraulic properties was accounted for by randomly generating a hydraulic conductivity field using a
one-dimensional first-order Markov process. Soil structural effects were modeled with a two-domain concept in which a firstorder kinetic expression is used to describe the transfer of water and solute between the two domains. Numerical experiments
were carried out for the case of furrow irrigation, including the breakthrough of a conservative solute to the groundwater table.
We compared five different scenarios: a single domain having uniform hydraulic properties (SU), a single domain with arandomly distributed hydraulic conductivity (SR), a dual-permeability system with uniform hydraulic properties (DU), a dualpermeability system with a randomly distributed fracture hydraulic conductivity (DRF), and a dual-permeability system having
a randomly distributed matrix hydraulic conductivity (DRM). All scenarios started with pressure heads in equilibrium with aconstant groundwater table 150 cm below the soil surface and zero initial solute concentrations. The simulated two-dimensional(2D) vertical concentration profiles showed preferential pathways resulting from both the spatial variability (SR) and soil
structure (DRF) scenarios. As expected, drainage of water from the bottom of the profile occurred significantly earlier for dualthan for single-permeability scenarios. The combination of having spatial variability in the hydraulic properties and invoking
the dual-permeability approach yielded the quickest and largest leaching of solute. The 2D dual-permeability approach should
considerably improve the simulation of water and solute movement in naturally heterogeneous field soils. 䉷 2000 Elsevier
Science B.V. All rights reserved.
Keywords: Dual-permeability media; Preferential flow; Unsaturated flow; Solute transport; Spatial variability; Soil structure* Corresponding author. Fax: ⫹420-2-2431-0782.
E-mail addresses: vogel@fsv.cvut.cz (T. Vogel), hgerke@zalf.de (H.H. Gerke), renduo@uwyo.edu (R. Zhang).
1
Fax: ⫹49-33432-82280.
2
Fax: ⫹307-766-6403.
3
Fax: ⫹909-342-4964.
0022-1694/00/$ - see front matter 䉷 2000 Elsevier Science B.V. All rights reserved.
PII: S0022-169 4(00)00327-9

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

1. Introduction
One-dimensional (1D) models assuming homogeneous soil properties often fail to accurately predict
solute movement at the field scale. Realistic models
must consider the horizontal and vertical variability insoil properties at a variety of spatial scales. Most fieldsoils exhibit different types of spatial heterogeneity,such as soil spatial variability and soil structure,
which often also coexist. Within the concept of continuum models (e.g. Bear, 1972), spatial variability
relates to the spatial distribution of macroscopic
model parameters, such as the hydraulic conductivity,
while in structured soils microscale effects sometimes
become so dominant that they affect macroscopic
scale flow and transport processes. In principle, both
spatial variability in soil hydraulic properties andstructure-induced heterogeneity can contribute to theinitiation of preferential pathways. The challenge is toadequately account for both types of spatial heterogeneity in simulation models.
Structured soils and fractured rock formations
contain a highly-permeable macropore or fracture
pore system through which water and dissolved
solutes can move at considerably higher velocities
than in the porous matrix. Consequently, local nonequilibrium conditions in the transient pressure headand solute concentration distributions may develop(we use the term ‘local’ here to indicate a sub-REV
scale, i.e. a scale smaller than that of the representative elementary volume (REV) in continuum
approaches). Such preferential phenomena severely
limit the prediction of water and solute movement in
structured or macroporous soils, as well as in fractured
rock systems. Preferential flow related to structure has
been widely reported in soils containing wormholes,
root channels, and inter-aggregate fissures (e.g.
Ehlers, 1975; Bouma, 1981; Beven and German,
1982; Brusseau and Rao, 1990; Wang, 1991; Thoma
et al., 1992; Flury et al., 1994; Liu et al., 1998; Pruess
1999; Vervoort et al., 1999; among many others).
Additional types of preferential flow have been
linked to textural differences rather than structuraleffects. Two types of preferential flow phenomenathat belong to this category are fingering (Hill andParlange, 1972; Hillel, 1987) and funneled flow
(Kung, 1990; Steenhuis et al., 1990; Walter et al.,2000). The evolution of finger-type preferential flow

79

paths is associated with gravity-driven flow instability(Raats, 1973; Parlange and Hill, 1976; Glass et al.,1989). Fingering occurs in water repellent soils,
when water percolates from a fine-textured into acoarse-textured layer, or when the air pressure
increases ahead of infiltration front (Hendrickx etal., 1993; Ritsema et al., 1993; Ritsema et al., 1998;
Dekker and Ritsema, 1994; Nieber, 1996; Bauters etal., 1998; DiCarlo et al., 1999).
Small-scale heterogeneities related to soil structure
can be modeled either by using a discrete fracture
network model (e.g. Sudicky and McLaren, 1992;
Therrien and Sudicky, 1996) or a multi-continuum
approach (e.g. Gwo et al., 1995, 1996). Within thediscrete fracture network concept, a map of the structural geometry must be known, while with the multicontinuum approach two or more continua, representing matrix and fracture systems, share the same spacedomain.
Heterogeneity in continuum models can be simulated using both single- and dual-porosity or -permeability approaches. Multi-dimensional single-porosity
models assume that preferential pathways may
develop as a result of spatially variable distributed
soil hydraulic properties. Spatially distributed soil
hydraulic properties, for example generated deterministically or as random functions of the spatial coordinates, lead to characteristic patterns of relatively
low and high flow velocities within the domain (e.g.
Roth, 1995; Roth and Hammel, 1996; Tsang et al.,
1988, 1996; Webb and Anderson, 1996; Birkhölzer
and Tsang, 1997). In this case, local equilibrium isassumed and a single Richards’ equation is solved.
The inclusion of randomly generated preferential
pathways has for instance been used by Hopmans etal. (1988) and Vogel and Hopmans (1992) to improvepredictions of drainage from a furrow-irrigated soil
profile.
The dual-porosity approach (Barenblatt et al.,1960), on the other hand, assumes that the medium
consists of two domains with different hydraulic andtransport properties. The term dual-permeability isused to indicate that flow takes place in both domains
in contrast to dual-porosity approach, which is often
used in the context of mobile-immobile-type solutetransport modeling. Water flow and solute transport
in dual-permeability models are described using separate flow and transport equations for the fracture and80

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

matrix pore systems (e.g. Dykhuizen, 1987; Jarvis etal., 1991; Gerke and van Genuchten, 1993a). Crucialcomponents of these types of models are transfer
terms governing the exchange of water and/or solutes
between the fracture and matrix pore systems. Empirical (Othmer et al., 1991) and semi-empirical (Gerkeand van Genuchten, 1993b) expressions exist that are
applicable to transient unsaturated flow.
Field scale studies of preferential flow fromglacial till agricultural soils to tile drains (Villholth et al., 1998; Villholth and Jensen, 1998)
show that both flow in macropores and spatial
variability may affect measured tile outflow andsolute output. Failure to describe the experimentally determined curves using the MACRO modelof Jarvis et al. (1991) was attributed to the factthat a hydraulic resistance at the aggregate surface
was not included. However, their modelingapproach did not allow testing the hypothesisthat the observed spatial variability may wellhave contributed to the discrepancies.
In addition to domain-internal heterogeneity, thegeometry of field boundaries may be irregular, thesoil surface shaped, and the boundary conditions
spatially distributed. While field-scale transport
processes are essentially 3D, they may, at times, be
simplified and reduced to 2D descriptions.
In this paper we extend the 1D approach of Gerkeand van Genuchten (1993a,b) to a 2D dual-permeability water flow and solute transport model bymodifying the SWM II finite element code ofVogel (1987). The model allows for 2D twodomain simulations involving water and solutetransfer between the fracture and matrix domains.
The most important feature of the model is its
ability to distinguish explicitly between heterogeneities due to macro-scale variability in the soil
hydraulic properties and heterogeneities caused by
micro-scale soil structural effects. Specificobjectives of this paper are to compare theeffects of a single- with those of a dualpermeability approach to simulate spatial heterogeneity, and to demonstrate the usefulness of
combining dual-permeability features with a 2Dmodel that considers spatially distributed hydraulicproperties. We present a simulation example toillustrate the interplay of the two types of subsurface heterogeneity.

2. 2D dual-permeability model
The model of Gerke and van Genuchten (1993a)for water and solute movement in a variably-saturated structured medium assumes that all properties of the bulk porous medium are composed of
two local properties, one associated with the fracture (subscript f) and one with the matrix(subscript m) pore system. Bulk and local properties are related by

1  wf 1f ⫹ wm 1m1a

u  wf uf ⫹ wm um

1b

q  wf qf ⫹ wm qm

1c

uc  wf uf cf ⫹ wm um cm

1d
⫺3

where 1 is porosity (L L ), u the water content(L 3L ⫺3), q the fluid flux density (LT ⫺1), c the soluteconcentration (ML ⫺3), and wf is the relative volumetric proportion of the fracture pore system, wm 
1 ⫺ wf : Flow of water in the dual-permeability
medium is described by means of two coupled
Richards’ equations as follows
3

Cf

2hf
G
 7· K f 7h f  ⫹ 7· K f 7z ⫺ w
2t
wf

Cm

2

2hm
G
 7· K m 7h m  ⫹ 7· K m 7z ⫹ w
2twm

where h is the pressure head (L), K the hydraulicconductivity tensor (LT ⫺1), C the specific water capacity (L ⫺1), z the vertical coordinate taken positive
upward (L), t time (T), and G w is the transfer term(T ⫺1) for water exchange between the two poresystems. The exchange of water between the matrix
and the fracture pore systems is assumed to be proportional to the pressure head difference between bothpore systems

G w  aw h f ⫺ h m 3

where a w is the first-order water transfer coefficient (L ⫺1T ⫺1). The local transient exchange ofwater causes a decrease in pressure in the drained
domain and a corresponding increase of pressurein the receiving domain consistent with therespective water retention curves and the

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

10 min

81

-90
-60
-30
0

Pressure head (cm)

6 hrs

12 hrsFig. 1. Pressure head profiles at 10 min, 6 and 12 h after the beginning of furrow irrigation in a single domain system with uniform hydraulicproperties (Scenario SU). The size of the flow domain is 2.0 × 1.5 m 2.

immediate flow conditions in the domains asprescribed by Eq. (2). Gerke and van Genuchten(1993b) obtained the following general expression
for the water transfer coefficient:aw  gwbKa
a2

equations as follows
2
u R c   7· u f Df 7c f  ⫺ 7· q f cf  ⫺ uf mf cf
2t f f f⫺4

where a is the characteristic radius or half-widthof the matrix structure (L), b the dimensionlessgeometry coefficient, Ka the hydraulic conductivity(LT ⫺1) of the matrix at or near the surface of thefracture pore system, and g w is a corrective coefficient. A single average value of 0.4 for g w was
found to be applicable for different hydraulicproperties and initial conditions (Gerke and van
Genuchten, 1993b). The geometry coefficient b
was evaluated by Gerke and van Genuchten(1996) to be 3 for rectangular slabs.
Similar to Eq. (2), the transport of solutes in a dualpermeability medium, with the solutes subject also tolinear equilibrium adsorption and first-order decay, isgoverned by two coupled convection–dispersion

Gs
wf

5

2
u R c   7· u m Dm 7c m 
2t m m m⫺7· q m cm  ⫺ um mm cm ⫹Gs
wm

where D is the dispersion coefficient tensor (L 2T ⫺1),
m a first-order decay coefficient (T ⫺1), R the retardation factor (dimensionless), and G s is the solute masstransfer term (ML ⫺3T ⫺1) evaluated as (Gerke and van
Genuchten, 1993b)

G s  ^G w cf ⫹ as wm um cf ⫺ cm 

6

where a s is the first-order solute mass transfer coefficient (T ⫺1). The first term on the right-hand side ofEq. (6) defines the advective contribution to G s, while
the second term gives the diffusive contribution. The
first-order mass transfer coefficient is of the form (Van

82

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

3.1. Soil hydraulic and transport properties

Genuchten and Dalton, 1986)

as 

b
Da
a2

7

where b and a are the same as for water, and Da is theeffective diffusion coefficient (L 2T ⫺1) of the soil
matrix at the fracture/matrix interface.

The soil hydraulic properties were described by theequations (Van Genuchten, 1980)u  ur ⫹ us ⫺ ur 1 ⫹ ⫺ahn ⫺m
1=m m 2
K Se   Ks S1=2e 1 ⫺ ⫺Se   ; m  1 ⫺ 1=n;Se  u ⫺ ur = us ⫺ ur 
3. Numerical experiments
As an illustrative example we consider furrow irrigation involving two furrows (Fig. 1). Water infiltrates into a 1.5-m deep soil profile from two
adjacent furrows over a 12-h period. The water levelin the furrows during the simulation period is kept
10 cm above the soil surface, thus resulting inconstant pressure head boundary conditions. The
water table at the bottom of the soil profile is also
kept constant, leading to a zero pressure head boundary condition. Zero-flux flow conditions are imposedon all remaining boundaries. The profile was assumed
to be initially at equilibrium with the water table(which means that all water fluxes were zero). The
initial condition for solute transport was taken to be
zero concentration, while along the wetted perimeterof the furrows the solute concentration of the infiltrating water was set equal to unity. A zero concentration
gradient was imposed along the bottom boundary,thus enabling solute to exit the profile with draining
water.
We compare five different flow scenarios: a single
domain with uniformly distributed soil hydraulicproperties (SU), a single domain with randomly
generated hydraulic conductivities (SR), a twodomain dual-permeability system with uniformly
distributed hydraulic properties (DU), a two-domaindual-permeability system consisting of a uniform
matrix and a fracture domain having randomly generated hydraulic conductivities (DRF), and a twodomain dual-permeability system with a uniform fracture domain but a randomized matrix domain (DRM).
The soil hydraulic parameters, the finite element grid,and the statistical parameters of the random conductivity field used in our numerical experiments, weretaken from the simulation study presented by Vogeland Hopmans (1992).

8

(9)

where u r and u s are the residual and saturated watercontents, respectively; Ks the saturated hydraulicconductivity (LT ⫺1), and a (L ⫺1) and n the empiricalfitting parameters. The dispersion coefficient forsolute transport was considered to be of the form
D  D0 t ⫹ l

兩q兩
u

10

where D0 is the molecular diffusion coefficient(L 2T ⫺1), l is the dispersivity (L), and t  u7=3 =u2s
(Millington and Quirk, 1961) is a tortuosity factor.
The hydraulic conductivity and dispersion tensors
were both assumed to be isotropic, and hence fully
determined by scalar functions.
The solute transport and soil hydraulic parameters,
except for the saturated hydraulic conductivity, wereassumed to be the same for all domains: D0 
0:01 cm2 min⫺1 ; R  1; m  0; l  5 cm; ur  0:24;
us  0:45; a  0:019 cm⫺1 ; and n  1:83: The average saturated hydraulic conductivities for the single,
matrix and fracture domains were taken to be equal to0.01, 0.0002, and 0.2 cm min ⫺1, respectively. Transfer parameters determining the exchange of water andsolute between the domains were set at wf  0:05;
Ka =Km  0:01; a  1 cm; gw  0:4; b  3; andDa =D0  0:01:
3.2. Random generation of the hydraulic conductivity
field
A log-normally distributed autocorrelated scaling
factor, defined as aK  Ks =K s ; was generated to create
subregions of high and low conductivity in the particular flow domain. A first-order Markov process(Haan, 1982) was used to randomize conductivities
in a horizontal directionyi⫹1  my ⫹ ry yi ⫺ my  ⫹ ti⫹1 s y 1 ⫺ r2y 1=2

11

0.01
0.30
0.60
0.90

Relative concentration

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

12 hrsFig. 2. Concentration profile at 12 h after the beginning of furrow
irrigation in a single domain system with uniform hydraulic properties (Scenario SU).

0.1
1
10

Scaling factor αK

where yi  ln aK is the value of y for a particular
nodal point i, m y and s y are the mean and standard
deviation, respectively, t is N{0,1}, and r y is the autocorrelation coefficient between y-values which are one
lag apart. The lag was assumed to be 3 cm, corresponding with the element size of the finite element
mesh. Parameters of the random process were: m y ⫺1:5; s y  1:7; and ry  0:8: The value of r y,
together with the lag, determines the shape of theMarkov autocorrelation function. The parameter r y
can be determined by fitting the autocorrelation function to a set of autocorrelation data, obtained by geostatistical analysis of a measured hydraulicconductivity field.
For this study we used only one realization of therandom process. The same realization was subsequently used to randomize the single-permeability
domain, the fracture domain, and the matrix domain.
The randomized saturated hydraulic conductivity field
is shown in Fig. 3. Notice that the high conductivity
regions create a pipe-like pattern. This feature allows
water to flow from the top part of the soil profile to theFig. 3. Randomly generated saturated hydraulic conductivity field in
the analyzed soil cross-section of 2.0 m width and 1.5 m depth. The
scaling factor, a K, is equal to the local conductivity divided by theaverage conductivity.

83

groundwater table without any blockage by lowconductivity regions. The topsoil and the bottom
parts of the flow region were not randomized. The
top part of the soil profile was kept homogeneous(e.g. due to cultivation) to prevent a direct connectionof the high-conductivity regions with water in the
furrows. The homogeneous bottom part assuresnatural horizontal equilibration of pressure heads insoil profile with the water table.
3.3. Simulation results
Figs. 1 and 2 illustrate how water and solute movein the single domain with uniformly distributed
hydraulic parameters (Scenario SU). Fig. 1 shows
the more or less classical development of the pressure
head profile during furrow irrigation. Water flow at
the end of simulation (12 h) is close to steady state,while the concentration front (Fig. 2) has reached adepth approximately half-way to the water table.
Results for the second case (Scenario SR) involving
a randomized conductivity field (Fig. 3), are presented
in Fig. 4. The hydraulic connection between the
furrows and the water table is now established much
faster. As a consequence, the solute concentration
front now reaches nearly the bottom of the profile atthe end of the simulation (12 h). The pipe-like hydraulic conductivity distribution leads to macro-scalewater flow channeling, with solute transport occurringin the form of a finger-type front.
Different results are obtained with the two-domain
simulation using a uniform distribution of soil hydraulic properties (Fig. 5). In this case (Scenario DU),
water flow in the fracture domain reaches steadystate relatively soon, while the water content in thematrix domain is increasing only relatively slowly,
predominantly through the absorption of water fromthe fractures. Also, as expected, solute transport rates
are significantly higher in the fracture than in thematrix pore system.
Fig. 6 present simulation results for the case of a
two-domain system in which the fracture domain has
a randomized hydraulic conductivity field (Scenario
DRF). This case shows the most heterogeneous distributions, especially for the solute concentration. As
compared to the previous example (DU), the watercontent in the fracture domain increased much faster.
The randomization apparently provides high-flux

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

Pressure head (cm)

0.90

0.30

0.01

12 hrs

0

-30

-60

-90

12 hrs0.60

84

Relative concentration

Fig. 4. Pressure head and concentration profiles in a single domain system with the randomly generated conductivity field (Scenario SR).

pathways for both water and the dissolved solute.
Notice the extreme finger-type concentration front inboth the fracture and matrix pore systems.
The DRM scenario illustrates the effect of employing a randomized conductivity field only for thematrix domain, while assuming the fracture system
to remain uniform (Fig. 7). The results are quite similar to those obtained for the case with two uniform
domains (DU). The pressure head profile in the matrix
system still exhibits a somewhat distinctive preferential flow pattern. However, the matrix pressure head
Fracture

distributions remain relatively isolated and ‘passive’in that they do not greatly affect solute displacementin the matrix domain.
Fig. 8a shows the calculated drainage rates during
the irrigation for the five scenarios. The water fluxes
for this purpose were integrated along the bottom
boundary. In case of two-domain simulations, thefluxes in the fracture and matrix domains werecombined using Eq. (1c) to yield one overall infiltration or drainage flux for the soil profile. Total cumulative amounts of water and solute entering and

-90
-60
-30
0
12 hrs

12 hrsMatrix
0.01
0.30
0.60
0.90

12 hrsRelative concentration

Fracture

Pressure head (cm)Matrix

12 hrsFig. 5. Pressure head and concentration profiles in a dual-permeability system with uniform hydraulic properties (Scenario DU).

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

Fracture

85

-90
-60
-30
0
12 hrsFracture

Pressure head (cm)Matrix

12 hrs

0.01
0.30
0.60
0.90
12 hrsRelative concentration

Matrix

12 hrsFig. 6. Pressure head and concentration profiles in a dual-permeability system, in which the saturated hydraulic conductivity of the fracturedomain was generated randomly (Scenario DRF).

Fracture

-90
-60
-30
0
12 hrs

12 hrsMatrix
0.01
0.30
0.60
0.90

12 hrsRelative concentration

Fracture

Pressure head (cm)Matrix

12 hrsFig. 7. Pressure head and concentration profiles in a dual-permeability system, in which the saturated hydraulic conductivity of the matrix
domain was generated randomly (Scenario DRM).

86

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

0.5
Solute flux (cm2/min)

Water flux (cm2/min)

3
DRF
2
DU
1
DRM
SR
0

0

4

8

SU

0.4

0.2
DRM
0.1
DU
0

12

DRF

0.3

0

Time (hrs)48Time (hrs)12
SU, SR

Fig. 8. Water and solute outflow from the lower boundary of the soil profile as computed for the respective flow scenarios: single uniform (SU),single with randomly generated conductivity (SR), dual uniform (DU), dual with randomly generated fracture conductivity (DRF), and dualwith randomly generated matrix conductivity (DRM).

leaving the particular domains of the profile during the
12 h of irrigation are listed in Table 1. The values forthe matrix and the fracture systems are provided separately after multiplication with their relative volumes,
i.e. wf and wm, respectively. This guarantees additivityof the matrix and fracture contributions to the overall
inflow/outflow amounts given in Table 1. The results
indicate that the inflow rates are quite similar for all
systems studied, but that the outflow rates varyconsiderably. The dual-permeability systems all
exhibited a faster response in terms of drainage ascompared to the single-permeability systems (Fig.
8a). The two-domain system with a randomized fracture domain (DRF) deviated the most from the homogeneous case by showing the earliest and largest
outflow rates. Randomizing the hydraulic conductivity field of the single-porosity system (SR) also
promotes early drainage, albeit to a lesser extent
than use of a dual-permeability system.
Finally, Fig. 8b shows solute leaching rates to thegroundwater table, obtained by integrating the solute
flux across the lower boundary. The results are consis-

tent with those in Fig. 8a in that the randomized fracture case (DRF) yields by far the largest amount ofsolute and the earliest breakthrough.

4. Discussion and conclusions
The dual-permeability approach used in this paper
is based on two separate but interacting 2D flow
domains. The model provides a useful tool for studying preferential flow in structured soils or unsaturated
fractured rock. Both domains of the dual-permeability
system can exhibit different spatial variabilities intheir hydraulic properties, as well as different degrees
of anisotropy. The further the two domains are
removed from mutual local equilibrium, the more
distinct the simulation results will be as compared tothe single-porosity approach. The degree of local
disequilibrium is controlled by the invoked values oftransfer term parameters.
The effects of randomization on the one hand, anddecomposition of the system into two domains on the

Table 1
Cumulative inflow and outflow amounts of water and solutes after 12 h of irrigation. The numbers in parenthesis indicate the value for thematrix multiplied by wm  0:95 and the numbers in brackets indicate the value for the fracture pore system multiplied by wf  0:05:
Flow system

Water/solute inflow (cm 2)

Water outflow (cm 2)

Solute outflow (cm 2)

SU
SR
DU
DRF
DRM

1896
1956
1717 (107) [1610]
2123 (118) [2005]
1823 (118) [1705]37500
762 (7) [755]1483 (13) [1470]
465 (4) [461]

0.000
0.009
12 (0.01) [12]207 (0.22) [207]
4 (0.00) [4]T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

other hand, are to a certain extent equivalent. The
decision of which approach is more appropriate
depends on the natural soil system being simulated,
the preferred modeling concept, and our ability torelate the soil hydraulic properties to either macroscale spatial variability or micro-scale structuraleffects. For example, fractured soil systems may bebetter modeled by using discrete fracture networkmodels, and finely aggregated or biopores containing
soils with the two-domain model. Unstructured soils
with relatively small-scale variabilities in the hydraulic properties may be described either by the distributed parameters approach or by the two-domain
approach. In addition to these conceptual criteria,
selection of the particular approach should also be
influenced by criteria of numerical efficiency. Macropore type heterogeneities at the scale of millimeters
are difficult to model by means of a single-porosity or-permeability approach. Randomization in that case is
not practical since the size of elements used for discretizing the flow domain would have to be significantly
smaller than the autocorrelation length of therandomly generated hydraulic properties. This
would require an excessive number of elements tocover the entire flow region.
Still, crucial to application of a multi-dimensional
dual-permeability approach is the need for independent determination of the dual set of soil hydraulic andtransport parameters, as well as of water and solutetransfer parameters characterizing soil structure andthe fracture/matrix interface. Since the parameters are
random variables, additional data on spatial distributions and autocorrelations are required. In many
natural systems structure-related preferential pathways are predominantly active at saturation and in anarrow range of negative pressure near saturation.
This makes it possible to apply relatively straightforward procedures to determine the separate hydraulicproperties of the fracture and matrix domains, forexample by executing a simple two-stage transient
flow experiment, in which during the first stage thepreferential pathways remain inactive (dry) and theprocess is dominated by the matrix, while during thesecond stage (involving saturation) the flow process isstrongly influenced by the preferential flow domain.
The first stage is then used to estimate the matrix
properties, while the second stage provides information necessary for the identification of the preferential

87

flow domain parameters. The parameters can then be
identified by solving the associated inverse problem
for the dual set of Richards and convection-dispersion
equations. Although simple in principle, the abovedescribed inverse procedure may suffer from lack of
uniqueness and likely needs to be complemented withadditional (independent) information on soil structure,
its spatial variability, and the fracture/matrix interface
properties. The problem of identifying the dual-porosity and dual-permeability soil hydraulic parameters
has been addressed by Durner (1994), Mohanty etal. (1997), Köhne et al. (1999), Vogel et al. (1999,
2000) and Schwartz et al. (2000), among others. A
variety of experimental techniques for dual-permeability parameter identification are currently being
developed by the authors.
Spatial variability in the soil hydraulic properties ofa multi-dimensional dual-permeability system can betreated in a deterministic (e.g. by means of soil stratification) or stochastic manner (e.g. using autocorrelated random fields). The dual-permeability model
described in this paper permits one to formulate andassess the usefulness of various alternative hypotheses
about the spatial distribution of hydraulic properties
by comparing model responses with real systemobservations.
Finally, we note that the 2D dual-permeability
model proposed in this paper is capable of generating
preferential flow patterns that qualitatively are in agood agreement with frequently observed distributions under field conditions (e.g. dye tracer studies
by Cislerová et al., 1990; Villholth et al., 1998; andLarsson et al., 1999). The extension of the dualpermeability approach to a 2D (or even 3D) systems
makes it possible to respect the natural spatial variability of hydraulic properties of the matrix as well as
the preferential flow (fracture) domain, which in turnfacilitates more realistic simulation of preferential
flow processes at the field scale. The approach should
help one to more effectively analyze the effects of soil
structure and soil spatial variability on preferential
flow, and design improved experiments for studying
preferential flow processes.

Acknowledgements
The research was partly conducted within the

88

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89

framework of the CTU research project CEZ J04/98
211100002. Additional support was provided byNATO Collaborative Linkage Grant EST 974796.

References
Barenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts
in the theory of seepage of homogeneous liquids in fissured
rocks. J. Appl. Math. Mech. 24 (5), 1286–1303.
Bauters, T.W.J., DiCarlo, D.A., Steenhuis, T.S., Parlange, J.-Y.,
1998. Preferential flow in water repellent sands. Soil Sci. Soc.
Am. J. 62, 1185–1190.
Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, NewYork.
Beven, K.J., German, P., 1982. Macropores and water flow in soils.
Water Resour. Res. 18, 1311–1325.
Birkhölzer, J., Tsang, C.-F., 1997. Solute channeling in unsaturated
heterogeneous porous media. Water Resour. Res. 33, 2221–
2238.
Bouma, J., 1981. Soil morphology and preferential flow along
macropores. Agric. Water Manage. 3, 235–250.
Brusseau, M.L., Rao, P.S.C., 1990. Modeling solute transport instructured soils: a review. Geoderma 46, 169–192.
Cislerová, M., Vogel, T., Šimunek, J., 1990. The infiltration–
outflow experiment used to detect flow deviations. In: Roth,K., Fluhler, H., Jury, W.A., Parker, J.C. (Eds.). Field-Scale
Water and Solute Flux in Soils. Monte Verita, Birkhauser
Verlag, Basel, pp. 109–117.
Dekker, L., Ritsema, C., 1994. How water moves in a water repellent sandy soil: 1. Potential and actual water repellency. Water
Resour. Res. 30, 2507–2517.
DiCarlo, D.A., Bauters, T.W.J., Darnault, C.J.G., Steenhuis, T.S.,
Parlange, J.-Y., 1999. Lateral expansion of preferential flow
path in sands. Water Resour. Res. 35 (2), 427–434.
Durner, W., 1994. Hydraulic conductivity estimation for soils withheterogeneous pore structure. Water Resour. Res. 30, 211–233.
Dykhuizen, R.C., 1987. Transport of solutes through unsaturated
fractured media. Water Res 21 (12), 1531–1539.
Ehlers, W., 1975. Observations on earthworm channels and infiltration on tilled and untilled loess soils. Soil Sci. 119, 242–249.
Flury, M., Flühler, M., Jury, W.A., Leuenberger, J., 1994. Susceptibility of soils to preferential flow of water: a field study. Water
Resour. Res. 30, 1945–1954.
Gerke, H.H., van Genuchten, M.Th., 1993a. A dual-porosity model
for simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res. 29 (2), 305–319.
Gerke, H.H., van Genuchten, M.Th., 1993b. Evaluation of a firstorder water transfer term for variably saturated dual-porosity
models. Water Resour. Res. 29 (4), 1225–1238.
Gerke, H.H., van Genuchten, M.T., 1996. Macroscopic representation of structural geometry for simulating water and solute
movement in dual-porosity media. Adv. Water Resour. 19 (6),343–357.
Glass, R.J., Parlange, J.-Y., Steenhuis, T.S., 1989. Mechanisms forfinger persistence in homogeneous, unsaturated, porous media:
theory and verification. Soil. Sci. 148, 60–70.
Gwo, J.P., Jardine, P.M., Wilson, G.V., Yeh, G.T., 1995. A multiple-pore-region concept to modeling mass transfer in subsurface
media. J. Hydrol. 164, 217–234.
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.
Haan, C.T., 1982. Statistical Methods in Hydrology. The Iowa State
University Press, Ames, IA.
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.
Hill, D.E., Parlange, J.-Y., 1972. Wetting front instability in layeredsoils. Soil Sci. Soc. Am. J. 36, 697–702.
Hillel, D., 1987. Unstable flow in layered soils: a review. Hydrol.
Processes 1, 143–147.
Hopmans, J.W., Schukking, H., Torfs, P.J.J.F., 1988. Two-dimensional steady unsaturated water flow in heterogeneous soils with
autocorelated soil hydraulic properties. Water Resour. Res. 24(12), 2005–2017.
Jarvis, N.J., Jansson, P-E., Dik, P.E., Messing, I., 1991. Modelling
water and solute transport in macroporous soil. I. Modeldescription and sensitivity analysis. J. Soil Sci. 42, 59–70.
Köhne, M., Meyer-Windel, S., Gerke, H.H., 1999. Parametrizing adual retention and conductivity model. Geophysical Research
Abstracts, EGS 24th General Assembly, The Hague, vol. 1, no.
2, p. 329.
Kung, K.-J.S., 1990. Preferential flow in a sandy vadose zone. 2.
Mechanisms and implications. Geoderma 46, 59–71.
Larsson, M.H., Jarvis, N.J., Torstensson, G., Kasteel, R., 1999.
Quantifying the impact of preferential flow on solute transport
to tile drains in a sandy field soils. J. Hydrol. 215, 116–134.
Liu, H.H., Doughty, C., Bodvarsson, G.S., 1998. An active fracture
model for unsaturated flow and transport in fractured rocks.
Water Resour. Res. 34 (10), 2633–2646.
Millington, R.J., Quirk, J.M., 1961. Permeability of porous solids.
Trans. Faraday Soc. 57, 1200–1207.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., van Genuchten,
M.Th., 1997. New piecewise- continuous hydraulic functions
for modeling preferential flow in an intermittent flood-irrigatedfield. Water Resour. Res. 33 (9), 2049–2063.
Nieber, J.L., 1996. Modeling finger development and persistence ininitially dry porous media. Geoderma 70, 209–229.
Othmer, H., Diekkrüger, B., Kutilek, M., 1991. Bimodal porosity andunsaturated hydraulic conductivity. Soil Sci. 152 (3), 139–150.
Parlange, J.-Y., Hill, D.E., 1976. Theoretical analysis of wetting
front instability in soils. Soil Sci. 122, 236–239.
Pruess, K., 1999. A mechanistic model for water seepage throughthick unsaturated zones in fractured rocks of low matrix permeability. Water Resour. Res. 35 (4), 1039–1051.
Raats, P.A.C., 1973. Unstable wetting fronts in uniform and nonuniform soils. Soil Sci. Soc. Am. Proc. 37, 681–685.
Ritsema, C.J., Dekker, L.W., Hendrickx, J.M.H., Hamminga, W.,
1993. Preferential flow mechanism in a water repellent sandy
soil. Water Resour. Res. 29, 2183–2193.

T. Vogel et al. / Journal of Hydrology 238 (2000) 78–89
Ritsema, C.J., Dekker, L.W., Nieber, J.L., Steenhuis, T.S., 1998.
Modeling and field evidence of finger recurrence in a water
repellent sandy soil. Water Resour. Res. 34 (4), 555–567.
Roth, K., 1995. Steady state flow in an unsaturated, two-dimensional, macroscopically homogeneous, Miller-similar medium.
Water Resour. Res. 31, 2127–2140.
Roth, K., Hammel, K., 1996. Transport of conservative chemical
through an unsaturated two-dimensional Miller-similar medium
with steady state flow. Water Resour. Res. 32, 1653–1663.
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.
Steenhuis, T.S., Kung, K.-J.S., Parlange, J.-Y., Selker, J.S., Chen,
X.X., 1990. Flow regimes in sandy soils with inclined layers. In:
Morel-Seytoux, H. (Ed.). Proceedings of the 10th Annual
Hydrology Days. Hydrology Days, Atherton, California, pp.
78–94.
Sudicky, E.A., McLaren, R.G., 1992. The Laplace Transform Galerkin technique for large-scale simulation of mass transport indiscretely fractured porous formations. Water Resour. Res. 28(2), 499–514.
Therrien, R., Sudicky, E.A., 1996. Three-dimensional analysis ofvariably-saturated flow and solute transport in discretely fractured porous media. J. Contam. Hydrol. 23, 1–44.
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 (5), 1357–1367.
Tsang, Y.W., Tsang, C.F., Neretnieks, I., 1988. Flow and tracer transport in fractured media: a variable-aperture channel model and its
properties. Water Resour. Res. 24, 2049–2060.
Tsang, Y.W., Tsang, C.F., Hale, F.V., Dverstorp, B., 1996. Tracer
transport in a stochastic continuum model of fractured media.
Water Resour. Res. 32, 3077–3092.
Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc.
Am. J. 44, 892–898.

89

Van Genuchten, M.Th., Dalton, F.N., 1986. Models for simulating
salt movement in aggregated field soils. Geoderma 38, 165–
183.
Vervoort, R.W., Radcliffe, D.E., West, L.T., 1999. Soil structuredevelopment and preferential solute flow. Water Resour. Res.
35 (4), 913–928.
Villholth, K.G., Jensen, K.H., Fredericia, J., 1998. Flow and transport processes in a macroporous subsurface-drained glacial till
soil. I. Field investigations. J. Hydrol. 207, 98–120.
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., 1987. SWM II-Numerical model of two-dimensional
flow in avariably saturated porous medium. Research report
no. 87, Department of Hydraulics and Catchment Hydrology,
Agriculture Univversity Wageningen, The Netherlands, 1987.
Vogel, T., Hopmans, J.W., 1992. Two-dimensional analysis offurrow infiltration. J. Irrig. Drain. Engng — ASCE 118 (5),791–806.
Vogel, T., Cislerová, M., Stodolovská, A., 1999. Hydraulic parameters of acid brown soil by inverting infiltration-outflow
experiment. Geophysical Research Abstracts, EGS 24th General
Assembly, The Hague, vol. 1, no. 2, p. 305.
Vogel, T., Šanda, M., CÍslerová, M., 2000. Modeling structure
related preferential flow component of subsurface runoff at thehillslope scale. GRA, EGS 25th General Assembly, Nice, vol. 2,
p. 272.
Walter, M.T., Kim, J.-S., Steenhuis, T.S., Parlange, J.-Y., Heilig,A., Braddock, R.D., Selker, J.S., Boll, J., 2000. Funneled flowmechanisms in a sloping layered soil: laboratory investigation.
Water Resour. Res. 36 (4), 841–849.
Wang, J.S.Y., 1991. Flow and transport in fractured rocks. Rev.
Geophys. (Suppl), 254–262.
Webb, E.K., Anderson, M.P., 1996. Simulation of preferential
flow in three-dimensional, heterogeneous conductivity fields
with realistic internal architecture. Water Resour. Res. 32,
533–545.