Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

MODELING NONEQUILIBRIUM CONTAMINANT TRANSPORT PROCESSESIN SOILS AND GROUNDWATER
Martinus Th. van Genuchten1,*, Elizabeth M. Pontedeiro2, and Renato M. Cotta1
1
Department of Mechanical Engineering, COPPE/LTTC, Federal University of Rio de Janeiro,
UFRJ, Rio de Janeiro, RJ, 21945-970, Brazil
2
Radioactive Waste Division, DIREJ, Brazilian Nuclear Energy Commission, CNEN, Rua da
Passagem 123, Botafogo, Rio de Janeiro, RJ, 20290-030, Brazil(*Corresponding author: rvangenuchten@yahoo.com)

ABSTRACT. Accurate simulation of contaminant transport processes in the subsurface remains amajor challenge in many science and engineering applications. In this paper we review a broad rangeof dual-porosity and dual-permeability formulations for modeling nonequilibrium or preferential fluid
flow and contaminant transport in macroporous soils or fractured rock. We focus especially ontransport processes in the variably-saturated vadose zone between the soil surface and the groundwatertable. One effective modeling approach for flow is to use composite functions for the unsaturated
hydraulic conductivity to account for the separate effects of macropores and micropores, and tocombine this approach with a mobile-immobile water type nonequilibrium formulation for solutetransport. Two example problems are given to illustrate the potentially important effects of preferential
flow on the simulation results. One application concerns pesticide transport in a tile-drained field forwhich only limited data were available for model calibration. A second application involves the longterm environmental fate of a radionuclide decay chain released from a mining installation in Amazonia
processing ore containing natural occurring radioactive materials.

INTRODUCTION
Despite significant advancements during the past several decades, accurate simulation of fluid flow
and contaminant transport in the subsurface remains a challenge in many cases. This is especiallytrue when trying to predict transport in natural field systems, which can be extremely heterogeneousat a range of scales. These heterogeneities often lead to preferential or nonequilibrium flow
processes that are very difficult to capture macroscopically with standard model formulations. The
standard models generally are based on the classical Richards equation for flow in variablysaturated media and the equilibrium advection dispersion equation for solute transport. A large
number of processes are known to contribute to preferential flow at the field scale. These include
the rapid movement of water and dissolved constituents through soil macropores or rock fractures,
unstable or gravity-dominated flow caused by soil textural changes or water repellency, andfunneling of water along inclined textural boundaries. Recent reviews of preferential flow andtransport processes and related modeling approaches are given by Hendrickx and Flury [2001],
Simunek et al. [2003] and Gerke [2006].
A variety of physical and chemical nonequilibrium models have been formulated over the years inattempts to parameterize preferential flow. In this paper we briefly review the classicalformulations for flow and transport in variably-saturated media, and then provide an overview ofavailable dual-porosity, dual-permeability and related formulations for modeling nonequilibrium or
preferential flow and transport in the subsurface. We also show two example applications. One

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

application concerns the transport of a pesticide (bentazone) in a large-scale tile-drained field forwhich only limited data were available for model calibration [Boivin et al., 2006]. The second
application involves the long-term environmental fate and transport of a radionuclide decay chain(238U Æ 234U Æ 230Th Æ 226Ra Æ 210Pb) released from a mining installation in Amazonia processingore containing natural occurring radioactive materials (NORMs).
CLASSICAL DESCRIPTIONS FOR FLOW AND TRANSPORT
Predictions of water and solute movement in the vadose zone are typically made using the classical
Richards equation for flow in unsaturated porous media and the advection-dispersion equations forsolute movement. For a one-dimensional vertical profile, these equations are as follows∂θ (h) ∂ ⎛∂h⎞= ⎜ K ( h) − K ( h) ⎟
∂t
∂x ⎝
∂x
⎠

(1)

∂θ c ∂ρ s ∂ ⎛∂c ⎞ ∂qc−φ+= ⎜θ D ⎟ ∂t∂t∂x ⎝∂x ⎠ ∂x(2)

respectively, where θ is the volumetric water content, h is the soil water pressure head, t is time, x is
distance from the soil surface downward, K is the unsaturated hydraulic conductivity as a functionof h or θ, c is the solution concentration, ρ is the solid phase bulk density, s is the adsorbed (solid
phase) concentration, φ is a general source-sink term, q is the volumetric fluid flux density given
by Darcy-Buckingham's law as

q = − K ( h)∂h
+ K ( h)∂x

(3)and D is the dispersion coefficient given by
D = Doτ + λ v(4)in which Do is the molecular or ionic diffusion coefficient, τ is the tortuosity factor, λ is thedispersivity, and v (=q/θ) is the average pore-water velocity. For simplicity all of the flow andtransport equations discussed in this paper are formulated for one-dimensional vertical flow throughthe vadose zone. Similar equations can be formulated for two- or three-dimensional problems[Simunek et al., 1999; Simunek and van Genuchten, 2006].
The hydraulic (constitutive) relationships θ(h) and K(h) needed in applications of Eq. (1) aredescribed here using the functions [van Genuchten, 1980]

θ ( h) = θ r +θs − θr[1+ |α h|n]m(m = 1 − 1/ n)
2

(5)

(6)
K (h) = K s Se ⎡⎣1 − (1 − Se1/ m ) m ⎤⎦where θr and θs denote the residual and saturated water contents, respectively, α and n are empiricalshape factors, Ks is the saturated hydraulic conductivity, and Se is effective fluid saturation given by

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

S e ( h) =θ -θr
θs - θr(7)

For all of the solute transport applications discussed in this study we assume linear equilibrium
partitioning between the solution and adsorbed phases (i.e., s = Kd c), in which case the transport
equation reduces to∂θ Rc ∂ ⎛∂c ⎞ ∂qc= ⎜θ D ⎟ ∂t∂x ⎝∂x ⎠ ∂x(8)

where R is the retardation factor:

R = 1+

ρ Kd
θ

(9)

The above flow and transport formulations are relatively standard in that they have been used forthe past several decades in various forms, simplifications or extensions in a large number of
research and engineering applications. The equations typically predict uniform flow and transport
processes in the subsurface that ignore the presence of physical non-equilibrium and preferential
flow. We next briefly review possible formulations that may be used to describe preferential flow
in soils and groundwater.

DUAL-POROSITY/PERMEABILITY FORMULATIONS
Preferential flow in macroporous soils and fractured rocks can be described using a variety of dualporosity or dual-permeability models. Figure 1 shows schematics of possible equilibrium andphysical nonequilibrium models for water flow and solute transport [Simunek and van Genuchten,
2008]. Figure 1a pertains to equilibrium flow and transport. The mobile-immobile nonequilibrium(Fig. 1b), dual-porosity (Fig. 1c) and dual-permeability (Fig. 1d) models typically assume that theporous medium consists of two overlapping but interacting pore regions, one associated with theinter-aggregate, macropore, or fracture pore system, and one comprising the micropores (or intraaggregate pores) inside soil aggregates or the rock matrix. Dual-porosity models assume that waterin the matrix is stagnant, while the more complex dual-permeability models also consider waterflow within the soil or rock matrix.
In the dual-permeability approach, Richards equations are applied to both the macropore or fracture
pore system and the matrix region. The flow equations for the fracture (subscript f) and matrix(subscript m) pore systems in that case are given by∂θ f (h f )
∂t

=

⎤ Γ
∂h f
∂ ⎡
− K (h f ) ⎥ - w⎢ K f (h f )
∂x ⎣⎢
∂x j
⎦⎥ w

∂θ m (hm ) ∂ ⎡∂h⎤=K m (hm ) m − K m (hm ) ⎥ + Γ w
⎢
∂t
∂x ⎣
∂x
⎦ 1- w

(10a)

(10b)

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

a. Uniform Flowb. Mobile-Immobile Water

Water

Water

Soluteθc. Dual-Porosity

Water
Immob.

SoluteImmob.

MobileMobileSoluteImmob.

θ = θim + θ mo

Mobileθ = θ im + θ mod. Dual-Permeability

Water
Slow

Fast

SoluteSlowFastθ = θM + θF

Figure 1. Conceptual models for equilibrium and preferential flow of water andsolutes. In the schematics, θ is the water content, θmo and θim in (b) and (c) arewater contents of the mobile and immobile flow regions, respectively, and θF (=wθf)and θM (=(1-w)θm ) in (d) are water contents of the matrix and macropore (fracture)regions, respectively (Simunek and van Genuchten, 2008).
respectively, where w is the ratio of the volumes of the fracture domain and the total porous
medium. Different dual-permeability formulations arise depending upon how water and solute
movement in the fractures and the matrix is modeled, and how water and solutes in the matrix andfracture regions are allowed to interact through the exchange term Γw [Gerke and van Genuchten,
1993a; Jarvis, 1999; Simunek et al., 2003]. In previous work we used a simple linear driving force
between the matrix and fracture regions based on discretization of Darcy-Buckingham’s law forhorizontal flow as follows [Gerke and van Genuchten, 1993a]:Γ w = α w (h f − hm )(11)in which αw is a first-order mass transfer coefficient for flow. An expression for αw can be derived
when the fractured or macroporous porous medium is assumed to have an idealized geometry
consisting of well-defined (e.g., rectangular or spherical) matrix blocks, in which case [Gerke andvan Genuchten, 1993b]:αw =β aγ w K a (h)
a2

(12)

where βa is a dimensionless geometry-dependent coefficient, a is the characteristic length of thematrix structure (e.g., half of the fracture spacing, or the radius of spherical aggregates), γw is adimensionless scaling coefficient, and Ka is the effective hydraulic conductivity of the fracturematrix interface. Gerke and van Genuchten [1993b] evaluated this conductivity asKa =[ K a (hm ) + K a (h f )]2

(13)

Alternative expressions for the exchange term Γw and Ka have been suggested by, among others,
Zimmerman et al. [1993]) and Köhne et al. [2004].

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

Analogous to Eqs. (10a,b), the dual-permeability formulation for solute transport is based on
advection-dispersion type equations for transport in both the fracture and matrix regions as follows(Gerke and van Genuchten, 1993a):∂θ f R f c f∂c f ⎞ ∂q f c fΓ∂ ⎛-φf − s⎜θ f D f⎟∂z ⎝∂z ⎠ ∂zw(14a)∂θ m Rm cm ∂ ⎛∂c ⎞ ∂q cΓ= ⎜ θ m Dm m ⎟ - m m - φm + s
1− w
∂t
∂z ⎝
∂z ⎠ ∂z

(14b)

∂t

=

where the subscripts f and m as before refer to the fracture and matrix pore regions, respectively,and Γs is the mass transfer term for solute exchange between the two regions:Γ s = α s (1 − w)θ m (c f − cm ) + Γ w c *(15)

in which αs is a first-order mass transfer coefficient for solute. Equation (15) shows that soluteexchange between the two liquid regions is modeled as the sum of an apparent first-order diffusion
process and advective transport, with c* being equal to cf when Γw>0 (advective transfer from thefracture to the matrix region) and cm when Γw<0 (advective transfer back to the fractures).
The above dual-permeability formulation presents a potentially very powerful model for simulating
preferential flow and transport processes in the subsurface. Example applications of dualpermeability models are given by Pot et al. [2005], Kodesova et al. [2005] and Köhne et al. [2006],
among others. While relatively complicated mathematically, the complexity is merely a reflectionof the complicated nature in which preferential flow occurs, and the many physical and chemical
processes that are known to contribute to preferential flow. And while much has been learned over
the years about preferential flow, the dual-permeability models still provide a very simplified
mathematical description of the processes involved.
When vertical water flow in matrix is ignored, and hence also no vertical advective transport occurs
in the matrix, the above dual-permeability formulations for water and solute reduce to those fordual-porosity models (Fig. 1c). Models of this type further reduce to the more classical two-region
mobile-immobile water type physical nonequilibrium transport models (Fig. 1b) when no water isallowed to exchange between the fracture and matrix regions (i.e., Γw = 0). Lateral advective soluteexchange between the fracture and matrix regions is then also neglected, but diffuse transfer willcontinue. Physical nonequilibrium models of this type are discussed in more detail in the next
section.
Application of dual-permeability models such as Eq. (10a,b) typically require two water retentionfunctions, θ(h), one for the matrix and one for the fracture pore system, and two or three
conductivity functions in terms of their local pressure heads, i.e., Kf(hf) for the fracture network,
Km(hm) for the matrix, and possibly a separate conductivity function Ka(ha) for the fracture/matrix
interface as embedded in the exchange term, Γw [Gerke and van Genuchten, 1993a]. Of thesefunctions, Kf is determined primarily by the structure of the fracture pore system (i.e., the size,
geometry and continuity of the fractures, and possibly the presence of fracture or macropore
coatings). Similarly, Km is determined by the hydraulic properties of single matrix blocks, and thedegree of hydraulic contact between adjacent matrix blocks during unsaturated flow. Suchinformation, unfortunately, is seldom available in practice. For these reasons a number of

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

simplifications are often invoked, such as instantaneous hydraulic equilibration between the fracture
and matrix pore systems.
TWO-REGION PHYSICAL NONEQUILIBRIUM MODELS

To avoid over-parameterization of the governing equations, one simplifying approach is to assume
instantaneous hydraulic equilibration between the fracture and matrix regions such that hf = hm (=h). In that case the coupling term Γw can be eliminated from Eqs. (10a,b) to recover Eq. (1), but
now with composite hydraulic properties of the form

θ (h) = wθ f (h) + (1 − w)θ m (h)(16a)

K (h) = wK f (h) + (1 − w) K m (h)(16b)

Many studies have used composite hydraulic properties of the type given by Eqs. (16a,b), especiallyfor the hydraulic conductivity function [Peters and Klavetter, 1988; Durner, 1994; Mohanty et al.,1997, Zurmühl and Durner, 1996; Vogel and Cislerova, 1998; de Vos et al., 1999]. While stillleading to uniform flow, models using such composite media properties do allow for faster flow andtransport during conditions near saturation, and as such may provide more realistic simulations offield data than the standard approach using unimodal hydraulic properties. In soils, the two parts ofthe conductivity curves may be associated with soil structure (near saturation) and soil texture (atlower negative pressure heads).
Very much consistent with Eqs. (16a,b) are composite hydraulic functions as described with themodel of Durner [1994] as follows
S e ( h) =θ ( h) − θ r
w1
w2
=
+
n m
θs − θr[1+ |α1h| ][1+ |α 2 h|n ]m1( w S + w S ) {w α [1 - (1 - Sl

K ( Se ) = K s

1 e1

2

e2

1 1

1

2

1 /m1 m1
e1

(17)

2

m

) ] + w2α 2 [1 - (1 - S 1e2/m2 ) 2 ]( w1α1 + w2α 2 )

2

}

2

(18)

where wi are the weighting factors for the two overlapping regions, and αi, ni, and mi are empiricalparameters of the separate hydraulic functions (i =1,2, …). An example of composite retention andhydraulic conductivity functions based on Eqs (17) and (18) is shown in Fig. 2 for the following set
of parameters: θr=0.00, θs=0.50, α1=0.01 cm-1, n1=1.50, l=0.5, Ks=1 cm d-1, w1=0.975, w2=0.025,
α2=1.00 cm-1, and n2=5.00. The fracture domain in this case represents only 2.5% of the entire porespace, but accounts for almost 90% of the hydraulic conductivity close to saturation.
Measurements of the composite (fracture plus matrix) hydraulic properties are greatly facilitated by
the use of tension infiltrometers. An advantage of tension infiltrometry is that negative soil waterpressures at the soil-infiltrometer interface can be maintained very close to zero, and that they canbe decreased in small steps to yield well-defined functions near saturation [e.g., Mohanty et al.,1997].

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

0.6
0.5

Water Content [-]Log(Conductivity [cm/days])0
Total
MatrixFracture

0.4
0.3
0.2
0.1
0

-2
Total
Matrix-4

Fracture
-6

-8

-10
-1

0

1

2

3

4

5

-1

Log(|Pressure Head [cm]|)

0

1

2

3

4

5

Log(|Pressure Head [cm]|)

Figure 2. Bimodal (composite) water retention (left) and hydraulic conductivity (right)
functions as described with the model of Durner [1994].
The use of composite hydraulic functions such as those shown in Fig. 2 is consistent with fieldmeasurements suggesting that the macropore conductivity of soils at saturation is generally about
one to two orders of magnitude larger than the matrix conductivity at saturation, depending upon
texture. These findings were recently confirmed by Schaap and van Genuchten, 2006) using adetailed neural network analysis of the UNSODA unsaturated soil hydraulic database [Leij et al.,
1996; Nemes et al., 2001]. The analysis revealed a relatively sharp decrease in the conductivity
away from saturation and a slower decrease afterwards. Schaap and van Genuchten (2005)
suggested an improved composite function for K(h) to account for the effects of macropores nearsaturation as follows:⎛ Ks ⎞K ( h) = ⎜⎟⎝ K m ( h) ⎠R(h)K m ( h)(19)h < -40 cm−40 ≤ h < -4 cm−4 ≤ h ≤ 0 cm

(20)

where
0
⎧
⎪
R (h) = ⎨0.2778 + 0.00694h⎪ 1 + 0.1875h⎩and where Km(h) is the traditional hydraulic conductivity function for the matrix as given by Eq. (6)but with a different exponent for the tortuosity term (the square root term) in Eq. (6). Equations(19) and (20) were found to produce very small systematic errors between the observed (UNSODA)
and calculated hydraulic conductivities across a wide range of pressure heads between saturationand -150 m. While the macropore contribution was found to be most significant between pressure
heads 0 and -4 cm, its influence on the conductivity function extended to pressure heads as low as 40 cm (Eq. 20).
If Eq. (1) is used in conjunction with composite hydraulic functions such as those given by Eqs(16a,b) or Eqs. (19) and (20), the solute transport equation for nonequilibrium flow will no longer
contain an advection term in the solute exchange term between the matrix and fracture regions. The
transport equations for nonequilibrium transport reduce then to a relatively standard dual-porosity
formulation which assumes that the liquid phase can be partitioned into mobile, θmo = wθf, andInt. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

immobile, θim = (1-w)θf, regions, with advective-dispersive transport being restricted to the mobile
region as follows [van Genuchten and Wagenet, 1989]:∂θ mo Rmo cmo ∂ ⎛∂c ⎞ ∂qcmo= ⎜ θ D mo ⎟ − α s (cmo − cim ) − φmo∂t∂x ⎝∂x ⎠ ∂x(21a)∂θi m Rim cim= α s (cmo − cim ) − φim
∂t

(21b)

in which
Rmo = 1 +f ρ Kdθ moRmo = 1 +(1 − f ) ρ K d

θim

(22)

where the subscripts mo and im refer to the mobile (macropore or fracture) and immobile (matrix)
regions, respectively, f is the dimensionless fraction of sorption sites in contact with mobile water, andαs is the solute mass transfer coefficient between the two regions. Please note that the exchange
term in Eq. (21a,b) does not contain the mobile water content, θmo, such as was the case with Eq.
(15). This was done to keep the same model formulation as used previously by van Genuchten and
Wierenga [1976] and van Genuchten and Wagenet [1989].
The totality of Eqs. (1), (16a,b) and (21a,b) provide a very convenient and parsimonious descriptionof nonequilibrium transport in soils and groundwater. One challenge still is the definition of themass transfer coefficient, αs. In earlier work [van Genuchten and Dalton, 1986; Gerke and van
Genuchten, 1993b] we related this value to geometric properties of the fractures or matrix. Whileuseful in concept, field soils are generally too complicated for simplified geometric analyses.
Moreover, evidence exists that the mass transfer coefficient αs decreases with the scale of thetransport problem. Based on analysis of previously published laboratory and intermediate fieldscale experiments, Maraqa [2001] found that αs is inversely related to the residence time of thecontaminant in the system as follows:

α s = atr− b(23)where αs is given in 1/hours, tr is the residence time in hours (tr=LR/vm in which L is the length of thetransport domain and vm the mobile region pore-water velocity given by q/θmo), a was 0.012 for stonysoils, 0.26 for aggregated media, and 0.95 for non-aggregated media; while b was found to be in therange 0.73 to 0.88. Equation (22) shows that αs is a scale-dependent parameter, roughly inversely
proportional (exactly if b=1) to the scale of an experiment. Scale dependency has long been known forseveral transport parameters such as the longitudinal dispersivity, λ, which is generally assumed toincrease with the length of the transport domain [e.g., Anderson, 1984]. Scale dependency may well
exist also macroscopically for the matrix diffusion coefficient [Liu et al., 2004].
The scale dependency of the mass transfer coefficient also makes sense intuitively if one accepts thenotion that some degree of nonequilibrium always will be present during transport in naturally
heterogeneous media, with the degree of non-equilibrium being controlled by new heterogeneities
continuously being encountered as the transport domain increases (e.g., Pachepsky, 2008). This notion
can be demonstrated by placing Eqs. (21a,b) in dimensionless form for steady-state flow in a
homogeneous medium as follows [van Genuchten and Wagenet, 1989]:

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

βR

∂c1
∂c
1 ∂ 2 c1 ∂c1*+ (1 − β ) R 2 =− ω (c1 − c2 ) − φmo∂T∂T P ∂z ∂x(1 − β ) R∂c2= ω (c1 − c2 ) − φim*
∂T

(24a)

(24b)

where
T=

qt
θL

β=

θ mo + ρ fK d θ mo Rmo=θ + ρ Kd
θR

X=

x
L

c1 =

c − CiCo − Ci

ω=

c2 =

αL
q

c − CiCo − Ci(25a,b,c)qLθ mo D(25d,e,f)

P=

in which L is the length of the transport domain, T and z are dimensionless time (often referred to as
pore volume) and distance, respectively, Ci and Co are the initial and inlet (at x=0) concentrations(assumed to be constants here), β accounts for the fraction of the porous medium containing the mobile
liquid and the sorption sites in equilibrium with mobile liquid (0 ≤ β ≤ l), ω is the dimensionless mass*transfer coefficient, and φmoand φim* are dimensional source-sink terms.
If the degree of nonequilibrium of the medium would not change with the spatial scale of the transport
problem, then the dimensionless mass transfer ω in Eq. (24b) would become a constant independent ofscale. Constancy of ω independent of scale (or the resident time) would imply an inverse relationshipbetween αs and L, consistent with Eq. (23) if b = -1. Assuming b = -1 and using the inventoried data of
Maraqa [2001] for both aggregated and intermediate field-scale data suggests the following
approximate relationship for αs:αs =q
8 Lθβ R

(26)

which is now independent of the units used for αs, q and L. The applicability of Eq. (26) definitely
needs further testing in future studies, especially the constant 8 in the denominator of Eq. (26).
APPLICATIONS

We now briefly summarize two applications of the physical nonequilibrium models. One application(Boivin et al., 2006) concerns pesticide transport in a large-scale tile-drained field for which only
limited data were available for model calibration [Boivin et al., 2006]. A second application[Pontedeiro et al., 2009] involves the long-term transport of a radionuclide decay chain released from a
conventional mining installation in Amazonia processing ore containing natural occurring radioactive
materials.
Example 1. Pesticide Transport in a Tile-Drained Field This example concerns the subsurface
transport of the herbicide bentazone in a tile-drained agricultural field in northeastern France. The
field consisted of two separate experimental sites (Fig. 3): Bouzule-1 made up of mostly silt loam,
and Bouzule-2 consisting of relatively fine-textured silty clay. Here we consider only transport at
the Bouzule-2 site, which was equipped with a subsurface tile drain system (0.05 m diameter), withtiles at an average depth of 0.9 m and a drain spacing of 8 m. Precipitation rates and othermeteorological variables were recorded at a weather station approximately 2 km from the field site.

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

Bouzule-1
Silt loam soil
(2.83 ha)Tile drain system(Ø = 5 cm)

Bouzule-2
Silty clay soil

Collector drain(Ø = 10 cm)(1.85 ha)Automatic samplers

Figure 3. Layout of the tile drainage systems at the Bouzule field sites [after Boivin et al., 2006].

Bentazone was sprayed onto the Bouzule-2 site on March 11, 2002. We used the HYDRUS-2D
software package [Simunek et al., 1999] to simulate the subsurface transport of this pesticide for aperiod of 100 days after application. The governing equations were exactly the same as before,
except for their extension to two dimensions [Boivin et al., 2006]. In the study we assumed isotropic
media, included provisions for root water uptake, but neglected pesticide degradation since laboratory
experiments indicated little or no degradation over time periods pertinent to the field experiments.
We refer to Boivin et al. [2006] for a detailed discussion of the field experiments, the invoked flow
and transport models, and the data used for the simulations.
Figure 4 shows the unstructured finite element mesh used in the simulations. The grid was stratified
in accordance with observed soil profile layering. Boundary conditions included no-flow
boundaries along the sides, atmospheric conditions at the soil surface (i.e., daily precipitation and

Depth (cm)Soil layers

0

1

20

2
Drain

3

120
0

100

200

400

Figure 4. Unstructured finite element grid used for the layered transport domain
representing half the drain spacing of the simulated Bouzule-2 silty clay site. The
grid consisted of 2298 triangular elements and 1237 nodes [Boivin et al., 2005].

250

Bouzule-2 Site

a

200
150

Observed
Simulated VGM

100

Simulated MVGM
50
0

Cumulative Drain Discharge(m3/ha)Drain Discharge Rate (m3/ha/d)Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

600

b

500
400

Observed
Simulated VGM

300

Simulated MVGM

200
100

3/10 3/15 3/20 3/25 3/30 4/4 4/9 4/14 4/19

0
3/10 3/15 3/20 3/25 3/30 4/4 4/9 4/14 4/19

Concentration of Bentazone(µg/L)

300
250
200
150

a
Observed
Simulated, ADE
Simulated, MIM

100
50
0
3/10 3/15 3/20 3/25 3/30 4/4 4/9 4/14 4/19

Exported Cumulative Amountof Bentazone (% applied)Figure 5. Observed and simulated instantaneous (a) and cumulative (b) tile draindischarge rates. Flow simulations were carried out with the standard (VGM)
hydraulic functions given by Eqs. (5) and (6), and the modified MVGM) functions
given by Eqs. (27) and (28) [Boivin et al., 2005].

4b3

Observed
Simulated, ADE

2

Simulated, MIM

1
0
3/10 3/15 3/20 3/25 3/30 4/4

4/9 4/14 4/19

Figure 6. Observed and simulated bentazone concentrations in the drainage water (a)and cumulative amount of bentazone exported (b) for the Bouzule-2 silty clay site.
Solute transport simulations were carried out with both the equilibrium advectiondispersion equation (ADE) and the mobile-immobile (MIM) physical nonequilibrium transport model [Boivin et al., 2005].

evaporation rates), a free drainage bottom boundary conditions, and a seepage face at the tile drain(which implies that the drains function only when the surrounding soil is saturated). Initial
simulations started on January 10, 2002, about 2 months before the pesticide application, with anequilibrium pressure head profile so as to obtain realistic soil water contents in the profile at thetime of the pesticide application.
Observed and measured tile drainage discharge rates and concentrations, as well as cumulative
amounts, are shown in Figs. 5 and 6. Observed data (Fig. 6a) indicated very high pesticideconcentrations in the drainage water (a maximum concentration of 223 µg L-1), substantially higher
than the limit of 0.1 µg L-1 set by the European Union for pesticide concentrations in groundwater[Focus Leaching Modelling Workgroup, 1995], and also above 3 µg L-1 used by the U.S.
Environmental Protection Agency for the maximum contaminant level for bentazone [USEPA,Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

2004]. Evidence of preferential flow at this site, especially during the normally relatively wet
spring and summer seasons, was demonstrated previously by several other experiments conductedwith a bromide tracer and the pesticide metolachlor (Novak et al., 2003)
HYDRUS-2D simulations using the traditional (VGM) hydraulic functions given by Eqs. (5) and(6) were found to severely underestimate the observed drain discharge rates (Fig. 5). Ks values
obtained from laboratory experiments did not account for the influence of drying cracks at the site,
and thus substantially underestimated actual Ks values in the field, especially during spring andsummer when drying developed in especially the soil surface horizon. To account for these cracks,
the modified (MVGM) soil hydraulic functions of Vogel et al. [2000] were used to mimic theeffects of preferential. These functions provide yet another way to account for the effects ofmacroporosity on variably-saturated flow, similarly as was shown in Fig. 2. The hydraulic functions ofVogel et al. [2000] are given by⎧⎪⎪⎪θ ( h) = ⎨
⎪
⎪
⎪⎩

θr+

θ m −θ rm(1+ |α h|n)h < hs(27)θsh ≥ hs⎧ ⎛ θ − θ ⎞l ⎡ 1 − F (θ ) ⎤ 2r⎪ Kk ⎜h ≤ hk
⎟ ⎢
⎥
θ
θ
−
⎪ ⎝ k
r ⎠ ⎣1 − F (θ k ) ⎦⎪⎪(h − hk )( K s − K k )⎪K ( h) = ⎨ K k +hk < h < hs
hs - hk⎪⎪⎪Ksh ≥ hs
⎪
⎪
⎪⎩

(28)

where l is a tortuosity factor in the conductivity term (assumed to be 0.5 in Eq. (6)), and⎡ ⎛ θ − θ ⎞1/ m ⎤rF (θ ) = ⎢1- ⎜
⎟ ⎥
θ
θ
−
⎢⎣ ⎝ m r ⎠ ⎥⎦

m

(29)

Equation (27) introduces a very small air entry value in θ(h) near saturation by replacing θs in Eq. (5)with a fictitious parameter θm slightly larger than θs. While this modification has little effect on theretention curve, the effect on the shape and value of the hydraulic conductivity function can beconsiderable, especially for fine-textured soils [Vogel et al., 2000; Schaap and van Genuchten, 2005].
Equation (28) assumes that the predicted hydraulic conductivity function is matched to a measured
value of the conductivity, Kk=K(θk), at some water content θk≤θs where Kk≤Ks. The conductivity
function is then assumed to increase linearly between hk and saturation to account for macropore flow.
Only a few of the hydraulic parameters in Eqs (27) and (28) were obtained from direct measurements(notably θs, θs and Kk). Others (α and n) were estimated from soil texture using the Rosetta
pedotransfer functions [Schaap et al., 1996], while Ks and hk were calibrated against the observed
drainage data.

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

Using Eqs. (27) and (28) with Ks values up to about 20 times the measured Kk values for some of thelayers, did lead to much more accurate simulations of the drain discharge rate during the spring of
2002. However, we acknowledge that it is very unlikely that the selected combination of Kk values
in the simulations represented a single unique description of the macropore flow process at the site.
Figure 6 shows that the concentrations and amounts of pesticide discharged from the tile-drainage
field were severely underestimated using the equilibrium (ADE) model. An excellent fit could beobtained with the mobile-immobile (MIM) physical nonequilibrium model upon calibration of themass transfer coefficient, αs, and equating the immobile water content, θr, to the residual watercontent, θr for all three layers. A single αs value of 0.027 d-1 for all three layers produced excellent
agreement with the measured concentrations (Fig. 5a) as well as with the total amount of bentazone
exported with the drainage water. The fitted value of 0.027 d-1 for as was well within the range ofvalues used by Maraqa [2001] for intermediate-scale field studies, but slightly lower than mostlaboratory studies analyzed by Maraqa [2001], including αs values predicted with Eq. (26).
The pesticide study shows that a limited amount of input data could be used to successfully
simulate drain discharge rates and chemical concentrations using an equilibrium flow model withcomposite hydraulic conductivity function, and the mobile-immobile physical nonequilibrium
model. The modified soil hydraulic functions in this example apparently provided sufficient
flexibility to permit accurate simulations at the site.
Example 2. Radionuclide Transport from a NORM Mining Waste Site In this example we
investigate the environmental fate of a radionuclides being released from a still operating
conventional mining installation in Amazonia processing ore containing natural occurring
radioactive materials (NORMs). The study was initially motivated by Brazilian regulations
requiring a performance assessment of the disposal facility using a leaching and off-site transportscenario. We used for our analysis the HYDRUS-1D software package [2005] to predict the longterm transport of radionuclides moving vertically through both the landfill and then laterally ingroundwater until a well 100 m downgradient from the site. The analysis was limited to theradionuclide decay chain 238U Æ 234U Æ 230Th Æ 226Ra Æ 210Pb. We refer to Pontedeiro et al. [2009]
for a more detailed description of the site, the hydraulic and transport parameters used in the
numerical analyses, and the modeling results.

For the calculations we assumed a 6-m thick homogeneous waste layer (overlain by a negligibly thinlayer of local soil), and a 5-m thick unsaturated zone. Radionuclides leaching vertically from the vadosezone were assumed to mix with laterally moving groundwater. The long-term average recharge rate at
the site was estimated to be 657 mm/y. This recharge rate was used as a surface flux boundary
condition for steady-state simulations of water and dissolved radionuclides flowing through the site andinto groundwater. The mixing process between vadose zone recharge water and laterally moving
groundwater was described using the EPA mixing zone model [USEPA, 1996]. Initial radionuclide
concentrations of the waste layer were estimated to be 71 Bq/g for 238U and 234U, 67 Bq/g for 230Th, 63
Bq/g for 226Ra and 4.8 Bq/g for 210Pb.
Transport calculations for this problem were previously carried out by Pontedeiro et al. [2009] for botha best-case scenario assuming equilibrium transport in a fine-textured (clay) subsurface, and a worstcase scenario with preferential flow in a very coarse-textured (sand) subsurface using the mobile
immobile physical nonequilibrium model. For those calculations we first used measured or estimatedhydraulic and transport properties of the slags and the local soil, a red saprolite known as Belterra clay
[Truckenbrodt and Kotschoubey, 1981, Belk et al., 2008]. The calculations involved very high Kd

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

values, taken from ISAM [1999], and hence retardation factors for the different radionuclides, leading
to relatively low transport rates over the assumed 100,000 year simulation period. We subsequently
contrasted those simulations with results for a coarse-textured subsurface. Hydraulic properties of thevadose zone below the waste layer were again described with Eqs. (5) and (6), but now usingparameters for a sand as predicted with the Rosetta pedotransfer functions [Schaap et al., 1998].
The transport of each radionuclide (subscript i) within the radionuclide decay chain was initiallymodeled using the equilibrium advection-dispersion equation:∂θ R1c1 ∂ ⎛∂c ⎞ ∂qc1= ⎜θ D 1 ⎟ − µ1θ R1c1∂t∂x ⎝∂x ⎠ ∂x∂θ Ri ci ∂ ⎛∂c ⎞ ∂qci= ⎜θ D i ⎟ + µi −1θ Ri −1ci −1 − µiθ Ri ci∂t∂x ⎝∂x ⎠ ∂x(30a)(i = 2,5)(30b)where the µi‘s represents first-order source-sink coupling between consecutive radionuclides within the
decay chain, with decay occurring in both the solution and adsorbed phases (accounted for by using
retardation factors in the decay terms).
Equations (29a,b) assume equilibrium transport. We also simulated the transport problem assuming
mobile/immobile type physical nonequilibrium transport. The governing transport equations for the i-th
radionuclide in the decay chain are then of the form [see also van Genuchten and Wagenet, 1989]∂θ Rmo ,i cmo ,i∂t=∂cmo,i ⎞ ∂qcmo ,i∂ ⎛− α s (cmo ,i − cim ,i ) + µi −1θ Rmo,i −1cmo,i −1 − µiθ Rmo ,i cmo ,i
⎜θ D
⎟∂x ⎝
∂x ⎠
∂x
∂θ Rim,i cim,i∂t= α s (cim,i − cim ,i ) + µi −1θ Rim,i −1cim ,i −1 − µiθ Rim ,i cim,i

(30a)

(30b)

where as before the subscripts mo and im refer to the mobile and immobile liquid phases. Equations(30a,b) allow simulation of nonequilibrium transport with a minimum of additional transportparameters. Following Nkedi-Kizza et al. [1983], f was set equal to the relative fraction of mobile
water (i.e., f = θm/θ), which makes these two parameters numerically equal to the dimensionalparameter β given by Eq. (25d). For the analysis we used αs values recommended by Maraqa [2001]
based on Eq. (26) assuming a and b values for non-aggregated media. Estimated values of αs for theBelterra clay calculations varied between approximately 0.0001 and 0.005 year-1, while those for thesand were approximately 10 to 100 times larger depending upon the soil layer (vadose zone or aquifer)and radionuclide involved. Dispersivities were assumed to be one tenth of the length of each
subdomain, while radionuclide half-life data were taken from ICRP [1983].
We first considered equilibrium transport in each of the three media (waste, vadose zone, andgroundwater). Calculated concentrations for Belterra clay are shown in Fig. 7. Please note the different
scales for the concentrations in these and subsequent figures. Results show that 238U and 234U have by
far the highest concentrations at both the bottom of the vadose zone and the well. However, some
contribution of 230Th is also clearly present in the leachate from the vadose zone, but not at the well 100
m downgradient from the site.

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

238

U Æ 234U Æ 230Th Æ 226Ra Æ 210Pb

40000

16000

35000

14000

30000

12000

25000

10000

20000

8000

15000

6000

10000

4000

5000

2000

0

0
0

50000

100000

0time (y)50000

100000time (y)Figure 7. Calculated radionuclide concentrations versus time leaving the vadose zone
(left) and at the well (right) assuming equilibrium transport (Belterra clay).

40000

14000

35000

12000

30000

10000

25000

8000

20000

6000
15000

4000

10000

2000

5000
0

0
0

50000
time (y)100000050000

100000time (y)Figure 8. Calculated radionuclide concentrations versus time leaving the vadose zone
(left) and at the well (right) assuming nonequilibrium transport (Belterra clay).
Figure 8 presents results when accounting for preferential flow through the waste and the Belterra clay
vadose zone and aquifer. Comparing Figs. (7) and (8) shows that the peak concentrations for uranium(having the highest concentrations) are only marginally affected by preferential flow (they are slightly
higher for preferential flow simulations). However, the peak concentrations arrive now much earlier atthe bottom of the vadose zone and especially at the well (after 43,000 years for preferential flow versus
83,000 years for equilibrium transport). The first radionuclides also arrived much earlier at the well
(12,500 years for preferential flow versus 22,000 years for equilibrium transport).
All of the simulations thus far were carried out for Belterra clay, leading to very conservative values forboth the hydraulic conductivity (low) and the estimated Kd‘s (high), thus leading to very low
contaminant transport rates. A complementary set of simulations was carried out to show the impact ofsoil texture on the long-term concentrations. For these simulations we assumed the presence of a very

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia
238

U Æ 234U Æ 230Th Æ 226Ra Æ 210Pb
5.00E+06

6e+0065e+006

4.00E+06

4e+0063.00E+06

3e+0062.00E+06
2e+006

1.00E+06

1e+006

0.00E+01

0
0

50000

0

100000

50000

100000time (y)time (y)Figure 9. Calculated radionuclide concentrations versus time leaving the vadose zone
(left) and at the well (right) assuming equilibrium transport (sand subsurface).

5.00E+06

4.00E+06
3.50E+06

4.00E+06

3.00E+06
2.50E+06

3.00E+06

2.00E+06

2.00E+06

1.50E+06
1.00E+06

1.00E+06

5.00E+05

0.00E+01

0.00E+01

0

50000
time (y)100000050000

100000time (y)Figure 10. Calculated radionuclide concentrations versus time leaving the vadose zone
(left) and at the well (right) assuming nonequilibrium transport (sand subsurface).
coarse-textured subsurface (sand) for both the vadose zone below the waste layer and groundwater,while leaving the physical and chemical properties of the waste layer itself unaltered. Kd were again
taken from ISAM (1999), while for the unsaturated hydraulic parameters in Eqs. (1) and (2) we selected
typical values for sand using Rosetta.
Results for equilibrium transport from the waste layer through the now sandy subsurface (vadose zone
and aquifer) are shown in Fig. 9. As compared to Belterra clay (Fig. 7), the 238U and 234U peaks
reached the bottom of the vadose zone much earlier (17,000 versus 24,000 years) with peakconcentrations of 5.95 106 Bq/m3, more than 100 times that for the clay soil (3.72 104 Bq/m3). This
reflects the much faster transport through the sandy vadose zone. Lateral transport rates through theaquifer were similarly very fast, with the 238U/234U peaks reaching the well 100 m downstream of thesite now after only 18,000 years with a peak concentration of 4.92 106 Bq/m3, compared to 85,000
years and a peak concentration of 1.4 104 Bq/m3 for the Belterra clay subsurface.

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

Finally, Figure 10 shows the effects of preferential flow on the equilibrium results for the sandy
subsurface. Comparing the results with those in Fig. 9 for equilibrium transport shows that, similarly as
for the Belterra clay subsurface, the peak concentrations were again only marginally higher for the
preferential flow scenario (e.g., 4.6 106 Bq/m3 for the sandy vadose zone compared to 3.7 104 Bq/m3 forthe clay). The peak concentrations however occurred much earlier, especially at the well (11,000 years
for the sand versus 43,000 years for the clay). This shows that both soil texture and preferential flow
have a major effect on the radionuclide transport rates.
CONCLUSIONSIn this paper we reviewed a broad range of dual-porosity and dual-permeability formulations formodeling nonequilibrium or preferential fluid flow and contaminant transport in macroporous soils or
fractured rock. We focused especially on transport processes in the vadose zone between the soil
surface and the groundwater table. One effective modeling approach requiring far fewer hydraulicparameters is to employ composite (dual-porosity type) functions for the unsaturated hydraulicconductivity to account for the separate effects of macropores and micropores, and then to combine thisapproach with a mobile-immobile water type nonequilibrium formulation for solute transport. The
approach was used successfully to simulate pesticide transport in a large-scale tile-drained field forwhich only limited data were available for model calibration. A second example application involved
the long-term environmental the radionuclide decay chain 238U Æ 234U Æ 230Th Æ 226Ra Æ 210Pb being
released from a still operating mining installation in Amazonia processing ore containing natural
occurring radioactive materials. The two examples show that preferential flow can have a major effecton the simulation results. The physical nonequilibrium formulation used in this study require an
estimate of the mass transfer coefficient (αs) governing diffusive exchange between the fracture andmatrix regions of structured media. While some guidance is provided by literature data from laboratory
and intermediate scale experiments, more research is needed to determine the exact scale dependencyof this parameter, especially for larger scale field applications.
REFERENCES

Anderson, M. P. [1984], Movement of Contaminants in Groundwater: Groundwater Transport Advection and Dispersion, In Groundwater Contamination, Studies in Geophysics, National
Academy Press, Washington, DC., pp. 37-45.
Belk, E. L., Markewitz, D., Radmussen, T. C., Carvaho, E. J. M., Nepstad, D. C., and Davidson, E. A.
[2007], Modeling the Effects of Throughfall Reduction on Soil Water Content in a Brazilian
Oxisol Under a Moist Tropical Forest, Water Resour. Res. Vol. 43, W08432, doi:
10.1029/2006WR005493.
Boivin, A., Simunek, J., Schiavon, M., and van Genuchten, M. Th. [2006], Comparison of Pesticide
Transport Processes in Three Tile-Drained Field Soils Using HYDRUS-2D, Vadose Zone J.,
Vol. 5, No. 3, pp. 838-849.
de Vos, J. A., Simunek, J., Raats, P. A. C., and Feddes, R. A. [1999], Identification of the HydraulicCharacteristics of a Layered Silt loam. In M. Th. van Genuchten, F. J. Leij, and L. Wu(eds.), Characterization and Measurement of the Hydraulic Properties of Unsaturated
Porous Media. Univ. of California, Riverside, CA. Part 1, pp. 783-798.
Durner, W. [1994], Hydraulic Conductivity Estimation for Soils with Heterogeneous Pore
Structure. Water Resour. Res. 30:211-223.
Focus Leaching Modelling Workgroup [1995], Leaching Models and EU Registration. GuidanceDocument 4952/VI/95. Commission of the European Communities, Directorate-General for
Agriculture VI B II-I, Brussels, Belgium.

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

Gerke, H. H. [2006], Preferential Flow Descriptions for Structured Soils, J. Plant Nutr. & Soil Sci.,
Vol. 169, pp. 382-400.
Gerke, H. H., and Köhne, J. M. [2004], Dual-Permeability Modeling of Preferential Bromide
Leaching from a Tile-Drained Glacial Till Agricultural Field, J. Hydrol.,Vol. 289, No. 1-4,
pp. 239-257.
Gerke, H. H., and van Genuchten, M. Th. [1993a], A Dual-Porosity Model for Simulating The
Preferential Movement of Water and Solutes in Structured Porous Media, Water Resour.
Res., Vol. 29, No. 2, pp. 305-319.
Gerke, H. H., and van Genuchten, M. Th. [1993b], Evaluation of a First-Order Water Transfer Term
for Variably-Saturated Dual-Porosity Models, Water Resour. Res., Vol. 29, No. 4, pp. 12251238.
Hendrickx, J. M. H., and Flury, M. [2001], Uniform and Preferential Flow Mechanisms in the Vadose
Zone, In Conceptual models of Flow and Transport in the Fractured Vadose Zone, National
Academy Press, Washington, DC, pp. 149-187.
ICRP [1983], Radionuclide Transformations: Energy and Intensity of Emissions, Int. Commission on
Radiological Protection, ICRP Publication 38, Pergamon Press, Oxford.
ISAM [1999], Derivation of Quantitative Acceptance Criteria for Disposal of Radioactive Waste to
Near Surface Facilities: Development and Implementation of an Approach, Safety Report,
Version 3.0, March 1999, IAEA, Vienna, Austria.
Jarvis, N. J. [1994], The MACRO Model (Version 3.1), Technical Description and Sample
Simulations, Reports and Dissertations 19, Dept of Soil Science, Swedish Univ. ofAgricultural Sciences, Uppsala, Sweden, 51 p.
Kodesova, R., Kozak, J., Simunek, J., and Vacek, O. [2005], Field and Numerical Study ofChlorotoluron Transport in the Soil Profile: Comparison of Single and Dual-Permeability
Model, Plant, Soil and Environ., Vol. 51, No. 7, pp. 310-315.
Köhne, J. M., Mohanty, B., Simunek, J., and Gerke, H. H. [2004], Numerical Evaluation of aSecond-Order Water Transfer Term for Variably Saturated Dual-Permeability Models,
Water Resour. Res., 40, doi:10.1029/2004WR00385.
Köhne, J. M., Köhne, S., and Simunek, J. [2006], Multi-Process Herbicide Transport in Structured
Soil Columns: Experiment and Model Analysis, J. Contam. Hydrol. Vol. 85, pp. 1-32.
Leij, F. J., Alves, W. J., van Genuchten, M. Th., and Williams, J. R. [1996], The UNSODA
Unsaturated Soil Hydraulic Database; User's Manual, Version 1.0. EPA/600/R-96/095,
National Risk Management Laboratory, Office of Research and Development, U.S.
Environmental Protection Agency, Cincinnati, OH, 103 p.
Liu, H. H., Bodvarsson, G. S., and Zhang, G. [2004], Scale Dependency of the Effective Matrix
Diffusion Coefficient, Vol. 3, pp. 312-315.
Maraqa, M. A. [2001], Prediction of Mass-Transfer Coefficient for Solute Transport in Porous Media.
J. Contam. Hydrol., Vol. 53, pp. 153-171.
Mohanty, B. P., Bowman, R. S., Hendrickx, J. M. H., and van Genuchten, M. Th. [1997], New
Piecewise-Continuous Hydraulic Functions for Modeling Preferential Flow in an
Intermittent Flood-Irrigated Field, Water Resour. Res., Vol. 33, No. 9, pp. 2049-2063.
Nemes, A., Schaap, M. G., Leij, F. J., and Wösten, J. H. M., [2001]. Description of the Unsaturated
Soil Hydraulic Database UNSODA, Version 2.0. J. Hydrol. Vol. 251, pp. 151-162.
Nkedi-Kizza, P., Biggar, J. W., van Genuchten, M. Th., Wierenga, P. J., Selim, H. M., Davidson, J. M.,
and Nielsen, D. R. [1983], Modeling Tritium and Chloride-36 Transport Through an
Aggregated Oxisol, Water Resour. Res., Vol. 19, No. 3, pp. 691-700.
Novak, S. M., Banton, O., and Schiavon, M. [2003], Modeling Metolachlor Exports in Subsurface
Drainage Water from Two Structured Soils under Maize (Eastern France), J. Hydrol., Vol.
270, No. 3-4, pp. 295-308.

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

Pachepsky, Y., Gimenez, D., Lilly, A., and Nemes, A. [2008], Promises of Hydropedology, CAB
Reviews: Perspectives in Agriculture, Veterinary Science, Nutrition and Natural Resources,
2008 3, No. 040.
Peters, R. R., and Klavetter, E. A. [1988], A Continuum Model for Water Movement in an
Unsaturated Fractured Rock Mass, Water Resour. Res., Vol. 24, pp. 416-430.
Pontedeiro, E. M., van Genuchten, M. Th., Cotta, R. M., and Simunek, J. [2009], The Effects ofPreferential Flow and Soil Texture on Risk Assessment of a NORM Waste Disposal Site, J.
Hazadous Materials (submitted).
Pot, V., Simunek, J., Benoit, P., Coquet, Y., Yra, A., and Martínez-Cordón, M.-J. [2005], Impact ofRainfall Intensity on the Transport of Two Herbicides in Undisturbed Grassed Filter Strip
Soil Cores. J. Contam. Hydrol., Vol. 81, pp. 63-88.
Schaap, M. G., Leij, F. J. and van Genuchten, M. Th. [1998], Neural Network Analysis for Hierarchical
Prediction of Soil Hydraulic Properties, Soil Sci. Soc. Am. J.,Vol. 62, No. 4, pp.847-855.
Schaap, M. G., and van Genuchten, M. Th. [2006], A Modified van Genuchten-Mualem Formulation
for Improved Prediction of the Soil Hydraulic Conductivity Function Near Saturation.
Vadose Zone J. Vol. 5, No. 1, pp. 27-34.
Simunek, J., and van Genuchten, M. Th. [2006], Contaminant Transport in the Unsaturated zone;
Theory and Modeling, In J.W. Delleur (ed.), Handbook of Groundwater Engineering, 2nded., CRC Press, Boca Raton, FL, pp. 22.1-46.
Simunek, J., and van Genuchten, M. Th. [2008]. Modeling Nonequilibrium Flow and Transport
Processes using HYDRUS, Vadose Zone J., Vol 7, pp. 782-797.
Simunek, J., Huang, K., Sejna, M., and van Genuchten, M. Th. [1999], 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, Int. Ground Water ModelingCenter, Colorado School of Mines, Golden, CO., 251 p.
Simunek, J., Jarvis, N. J., van Genuchten, M. Th., and Gärdenäs, A. [2003], Review andComparison of Models for Describing Non-Equilibrium and Preferential Flow andTransport in the Vadose Zone. J. Hydrol. Vol. 272, No. 1-4, pp. 14-35.
Simunek, J., van Genuchten, M. Th., and Sejna, M. [2005], The HYDRUS-1D Software Package forSimulating the One-Dimensional Movement of Water, Heat and Multiple Solutes in VariablySaturated Media. Version 3.0. HYDRUS Software Series 1, Dept. of Environmental Sciences,
Univ. of California, Riverside, 240 p.
Simunek, J., van Genuchten, M. Th., and Sejna, M. [2008], Development and Applications of the
HYDRUS and STANMOD Software Packages and Related Codes, Vadose Zone J., Vol. 7, pp.
587-600.
Truckenbrodt, W., and Kotschoubey, B. [1981], Cobertura Terciária das Bauxitas Amazônicas., Revista
Braileira Geociências, Vol. 11, pp. 203-208.
USEPA [1996], Soil Screening Guidance: Technical Background Document, EPA/540/R-95/128,
PB96–963502, Office of Solid Waste and Emergency Response, U.S. Environmental
Protection Agency, Washington, DC.
USEPA [2004], Drinking Water Standards and Health Advisories, Office of Water, U.S.
Environmental Protection Agency, Washington, DC, 20 p.
van Genuchten, M. Th. [1980]. A Closed-Form Equation for predicting the Hydraulic Conductivity of
Unsaturated Soils, Soil Sci. Soc. Am., Vol. 44, pp. 892-898.
van Genuchten, M. Th., and Dalton, F. N. [1986], Models for Simulating Salt Movement in AggregatedField Soils, Geoderma, Vol. 38, pp. 165-183.
van Genuchten, M. Th., and Wagenet, R. J. [1989], Two-Site/Two-Region Models for Pesticide
Transport and Degradation: Theoretical Development and Analytical Solutions. Soil Sci. Soc.
Am. J., Vol. 53, No. 5, pp. 1303-1310.

Int. Symp. on Convective Heat and Mass Transfer in Sustainable Energy
April 26 – May 1, 2009, Tunisia

van Genuchten, M. Th., and Wierenga, P. J. [1976], Mass Transfer Studies in Sorbing Porous Media. I.
Analytical Solutions. Soil Sci. Soc. Am. J., Vol. 40, pp. 473-480.
Vogel, T., and Císlerova, M. [1988], On the Reliability of Unsaturated Hydraulic Conductivity
Calculated from the Moisture Retention Curve, Transp. Porous Media, Vol. 3, pp. 1-15.
Vogel, T. van Genuchten, M. Th., and Cislerova, M. [2000], Effect of the Shape of the Soil Hydraulic
Functions near Saturation on Variably-Saturated Flow Predictions. Adv. Water Resour. Vol. 24,
No. 2, pp. 133-144.
Zimmerman, R. W., Chen, G., Hadgu, T., Bodvarsson, G. S. [1993], A Numerical Dual-Porosity
Model with Semi-Analytical Treatment of Fracture/Matrix Flow, Water Resour. Res., Vol.
29, pp. 2127–2137.
Zurmühl, T., Durner, W. [1996], Modeling Transient Water and Solute Transport in a Biporous
Soil, Water Resour. Res., Vol. 32, pp. 819–829.