2Integrated modeling of vadose-zone flow and transport
processes
M.T. van Genuchten and Jirka ŠimĤnek
Abstract
Enormous advances have been made during the past several decades in our
understanding and ability to model flow and transport processes in the vadose zone
between the soil surface and the groundwater table. A large number of conceptual models
are now available to make detailed simulations of transient variably-saturated water flow,heat movement and solute transport in the subsurface. In this paper we highlight four
examples illustrating such advances: (1) coupling physical and chemical processes, (2)simulating colloid and colloid-facilitated transport, (3) integrated modeling of surface
and subsurface flow processes, and (4) modeling of preferential flow in the subsurface.
The examples show that improved understanding of underlying processes, continued
advances in numerical methods, and the introduction of increasingly powerful computers
now permit us to make comprehensive simulations of the most important coupled, nonlinear physical, chemical and biological processes operative in the unsaturated zone.

Introduction
Agricultural practices for supplying food to locally burgeoning populations were
probably first implemented some 10,000 years ago in Eurasia, Africa and the Americas(Lawton and Wilke 1979). Hunting activities and gathering food were slowly
complemented, or even replaced, by efforts to cultivate the Earth. This involved
modifying the landscape to grow agricultural crops, raising domesticated animals and,increasingly between roughly 8,000 and 3,000 years ago, practicing irrigation. Many
parts of the world still show impressive remnants of ancient irrigation systems, such as
chains-of-wells systems in China, the Middle and Near East, and North Africa. Thesesystems from around 2500 years ago and onward, supplied water for domestic
consumption and irrigation through sophisticated subsurface irrigation tunnels, also
known as qanatz, kanats or karez (e.g. Cressey 1958; Lawton and Wilke 1979). Some ofthese chains-of-wells systems are still in use today in the Middle East, China andelsewhere.
The ancient irrigation structures and related soil-management and irrigation
practices, were at the time implemented without the advantage of modern computational
tools, such as computers and sophisticated numerical methods, and hence withoutGeorge E. Brown, Jr. Salinity Laboratory, USDA, ARS, 450 West Big Springs Road, Riverside, CA,
92557, USA. E-mail: rvang@ussl.ars.usda.gov
Department of Environmental Sciences, University of California, Riverside, CA, 92521, USA. E-mail:jiri.simunek@ucr.edu

37

Chapter 2
detailed conceptual-numerical descriptions based on equations governing water flow andchemical transport into and through the unsaturated zone. Against the long human historyof manipulating soils and water for improved agricultural practices, it is therefore
amazing that Darcy’s law for saturated flow was first proposed only some 150 years ago(Darcy 1856), while the Richards equation for unsaturated flow was first formulated
barely 70 years ago (Richards 1931). This shows that, compared to the time periods ofhuman history and the implementation of agricultural practices and irrigation, truly
astonishing advancements have been made during the past 30 years or so in processbased descriptions of subsurface flow and transport processes, numerical analysis andcomputer hardware. We have come to a point in history where truly unique opportunities
exist for developing models that integrate the most pertinent processes affecting water,
heat and solute movement in the subsurface. Much of the research in subsurface flow andtransport has historically progressed along mostly disciplinary lines (e.g., soil physics,
hydrogeology, geochemistry, microbiology, plant physiology). As such, physical,
chemical and microbiological processes were often studied and implemented in relative
isolation. The same is true for studies of surface and subsurface flow processes, and forflow and transport in the vadose zone and in groundwater. The introduction ofincreasingly powerful computers, advanced numerical methods and improved
understanding of subsurface flow and transport processes, now provide tremendous
opportunities for integrating the various processes involved.
In this paper we highlight several exciting approaches for coupling physical,
chemical and other processes into integrated descriptions of the subsurface, while
focusing especially on the vadose zone. We will address (1) relatively standard
approaches for modeling flow and transport, (2) multicomponent geochemical transport,(3) colloid and colloid-facilitated transport, (4) integrated surface/subsurface modeling,and (5) process-based descriptions of preferential flow. We start with a relatively
standard description of water, heat and solute movement.

Standard descriptions for water, heat and solute movement
Predictions of water, heat and solute movement in the vadose zone are traditionallymade using the Richards equation for variably-saturated water flow and advection–
dispersion type equations for heat and solute movement. For a one-dimensional soil
profile these equations are given by e.g., ŠimĤnek, Šejna and Van Genuchten (1998)

wT
wt

w
wh
[ K ( h)- K (h)] -S
wz
wz

w C p (T )T
wt

w
wT
wqT
[ O (T ) ] - C w- C w STwz
wz
wz

w( U s) w(T c) w
wc+
= (T D - qc) -I
wt
wt
wz
wz

(1)

(2)

(3)

respectively, where T is the volumetric water content (-), h is the soil water pressure head(L), t is time (T), z is distance from the soil surface downward (L), K is the hydraulicconductivity (LT-1) as a function of h or T, O is the apparent thermal conductivity of thesoil (MLT-3K-1), Cp and Cw are volumetric heat capacities (ML-1T-2K-1) of the porous38

Van Genuchten and ŠimĤnek

medium and the liquid phase, respectively; T is temperature (K), s is the soluteconcentration associated with the solid phase of the soil (MM-1), c is the soluteconcentration of the liquid phase (ML-3), U is the soil bulk density (ML-3), D is the solute
dispersion coefficient (L2T-1), q is the volumetric fluid flux density (LT-1) given by
3
Darcy’s law, and S (T-1) and I (M L T-1) are sinks or sources for water and solutes,
respectively.
Equations (1) through (3) are relatively standard in that they have been popularly
used for the past 50 years or so in various forms, simplifications or extensions. Examplesof possible extensions are their use in two or three dimensions (ŠimĤnek, Huang and Van
Genuchten 1995; Zyvoloski et al. 1997; Pruess, Oldenburg and Moridis 1999), theinclusion of vapor-phase transport (Scanlon et al. 2003) and volatilization (e.g. Jury,
Spencer and Farmer 1983), extensions to multiphase flow involving separate flow
equations for each fluid involved, such as water, air and/or oil (e.g. Pinder and Abriola
1986; Pruess and Battistelli 2002), multicomponent geochemical transport (Yeh andTripathi 1990; Lichtner 2000), and incorporation of dual-porosity or dual-permeability
formulations for preferential flow in macroporous soils or fractured rock (e.g. Gerke and
Van Genuchten 1993; ŠimĤnek et al. 2003; Bodvarsson, Ho and Robinson 2003). Severalof these and related studies have resulted in detailed codes such as TOUGH2 (Pruess,
Oldenburg and Moridis 1999), the HYDRUS codes (ŠimĤnek, Šejna and Van Genuchten
1998; 1999), FEHM (Zyvoloski et al. 1997), HYDROBIOGEOCHEM (Yeh et al. 1998),
RZWQM (Ahuja et al. 1999) and SWAP (Van Dam et al. 1997), among many others.
Much of the work presented in this paper is based on the HYDRUS codes (ŠimĤnek,Šejna and Van Genuchten 1998; 1999), but further modified to enable simulations of thefour types of problems considered here (i.e., multicomponent geochemical transport,
colloid and colloid-facilitated transport, integrated surface/subsurface modeling, andpreferential/non-equilibrium flow and transport). For this reason we first give a briefsummary of these codes. The HYDRUS-1D and HYDRUS-2D software packages are
finite-element numerical models for simulating the one- or two-dimensional movementof water, heat and multiple solutes in variably saturated media. The programs
numerically solve the Richards equation for variably-saturated water flow and advection–
dispersion equations for heat and solute transport.
The solute transport equations in HYDRUS consider advective–dispersive transport
in the liquid phase, and diffusive transport in the gaseous phase. The transport equations
also include provisions for non-linear and/or non-equilibrium reactions between the solidand liquid phases, linear equilibrium reactions between the liquid and gaseous phases,
zero-order production and two first-order degradation reactions: one which is
independent of other solutes, and one which provides the coupling between solutes
involved in sequential first-order decay reactions. Figure 1 shows a schematic of thesequential first-order decay chain reactions incorporated in the HYDRUS codes. Typicalexamples of sequential first-order decay chains are radionuclides (e.g. Lester, Jansen andBurkholder 1975; Rogers 1978; Van Genuchten 1985), nitrification/denitrification chains(e.g. Misra, Nielsen and Biggar 1974; Wagenet and Hutson 1987; Ahuja et al. 1999),
organic-phosphate transport (Castro and Rolston 1977), pesticide decay chains (e.g.
Bromilow and Leistra 1980; Wagenet and Hutson 1987), chlorinated-hydrocarbon decay
chains (PCE to TCE to ethylene) (Schaerlaekens et al. 1999; Casey and ŠimĤnek 2001),
pharmaceuticals (Casey et al. 2003; 2004), explosives and energetics (Selim, Xue andIskandar 1995; Sheremata et al. 1999) and other contaminants.

39

Chapter 2

 Products
Pg,1

Pw,1

g1

A
c1

Ps,1

Jw,1

Pg,2

Pg,1’
Pw,1’
Ps,1’

s1

kg,1 ks,1

Jg,1

Products

Js,1

 Products

Pw,2

g2

B
c2

kg,2

Jg,2

Ps,2

Jw,2

s2

ks,2

C ...

Pg,2’
Pw,2’
Ps,2’

Js,2

Products

Figure 1. Schematic of system of solutes subject to zero-order production/decay andsequential first-order degradation chain reactions, where c (ML-3), s (MM-1), and g (ML3
) represent concentrations in the liquid, solid and gaseous phases, respectively; thesubscripts s, w, and g refer to solid, liquid and gaseous phases, respectively; straightarrows represent the different zero-order (Ȗ) and first-order (µ, µ') rate reactions, and kg
and kd indicate equilibrium distribution coefficients between phases
Incorporation of the decay chain reactions of Figure 1 into Eq. (3) leads to thefollowing set of transport equations (ŠimĤnek, Šejna and Van Genuchten 1998):
wT ck wU sk wav g k


wt
wt
wt

w g · wqck
wc ·
w §
w §
T D kw k ¸ 
av D kg k ¸ ¨
¨
wz ©
wz ¹
wz ©
wz ¹ wz-( P w, k  P w' , k )T ck - ( P s , k  P s' , k ) U sk - ( P g ,k  P g' ,k )av g k  P w' ,k 1T ck 1 
 P s' ,k 1 U sk 1  P g' , k 1 av g k 1  J w, kT  J s ,k U  J g , kav - S cr ,kk2,3, ..., ns(4)where Pw, Ps , and Pg are first-order rate constants (T-1) for solutes in the liquid, solid andgas phases, respectively; Pw', Ps', and Pg' are similar first-order rate constants (T-1)
providing connections between individual chain species, Jw , Js , and Jg are zero-order rate
constants for the liquid (ML-3T-1), solid (T-1), and gas (ML-3T-1) phases, respectively; av
is the air content, S is the sink (T-1) term in the water flow equation (1), cr is theconcentration of the sink term, Dw is the dispersion coefficient for the liquid phase (L2T1), and Dg is the diffusion coefficient for the gas phase (L2T-1). As before, the subscripts
w, s, and g correspond with the liquid, solid and gas phases, respectively; the subscript krepresents the kth chain number, and ns is the number of solutes involved in the chain
reaction. While seemingly complicated because of the many terms involved, the basic
concepts of Eq. (4) are relatively straightforward. Also, application of models based on
the above equations is greatly facilitated by using graphical interfaces (GUI), such as
those implemented in the HYDRUS software packages (www.hydrus2d.com).

40

Van Genuchten and ŠimĤnek

Multicomponent geochemical transport
Equation (3) provides a highly simplified approximation for solute transport in many
applications in that the model considers the transport of only a single species. Equation(4) is already more realistic by assuming that the solutes are subject to relatively simple
consecutive decay chain reactions. Still, in general, both models ignore the fact that thesoil liquid phase always contains a mixture of many ions that may interact, create
complex species, precipitate, dissolve and/or compete with each other for sorption sites
on the solid phase. The fate and transport of many naturally occurring elements andcontaminants in the subsurface is affected by a multitude of complex, interactive
physical, chemical, mineralogical and biological processes. For example, changes in thechemical composition or pH of the soil solution may impact the retention of heavy metals
on organic matter or iron oxides. Dissolution and precipitation of minerals generally
buffer the transport of a solution with a different pH through the soil profile. Simulationof these and related processes requires the use of a coupled multicomponent reactive
transport code that integrates the physical processes of water flow and advective–
dispersive solute transport with a range of possible biogeochemical processes. Below wegive first a brief overview of multicomponent transport-modeling approaches, and thenillustrate their use with two related examples, one for steady-state flow and one forvariably-saturated conditions.
Overview of past and current research
Many important environmental problems require analysis of the coupled transportand reaction of multiple chemical species (ŠimĤnek and Valocchi 2002). Examples areacid mine drainage (Walter et al. 1994; Lichtner 1996), radionuclide transport(Viswanathan et al. 1998), the fate and transport of metal-organic mixed waste (Rittmann
and VanBriesen 1996; VanBriesen 1998), analysis of redox zones in organiccontaminated aquifers (Abrams, Loague and Kent 1998; Essaid et al. 1995), and reactive
permeable barriers for aquifer remediation (Fryar and Schwartz 1994). Multicomponent
geochemical transport may need to be considered also in applications dealing withsalinity management and evaluation of irrigation or drainage practices in arid andsemiarid areas, including studies of the suitability of water for irrigation, drainage waterreuse and the reclamation of sodic soils (e.g. Tanji 1990; Rhoades 1997; ŠimĤnek andSuarez 1997).
Multicomponent geochemical transport studies were initially limited mostly to thesaturated zone where changes in the fluid velocity, temperature and pH are generally
much more gradual and hence less important than in the unsaturated zone. Consequently,
most multicomponent transport models assume one- or two-dimensional steady-state
saturated water flow with fixed values for the flow velocity, water content (or porosity),
temperature and/or pH (Valocchi, Street and Roberts 1981; Rubin 1983). Early examples
for more general conditions involving major ion chemistry and cation exchange inunsaturated soils, generally using mixing cell approaches, are given by Dutt, Shaffer andMoore (1972), Robbins, Wagenet and Jurinak (1980), and Wagenet and Hutson (1987).
Only recently have more generalized multicomponent transport models been developed
which also consider variably-saturated flow. These include DYNAMIX (Liu andNarasimhan 1989), HYDROBIOGEOCHEM (Yeh and Tripathi 1990; Yeh et al. 1998),
UNSATCHEM-2D (ŠimĤnek and Suarez 1994), FEHM (Zyvoloski et al. 1997),
MULTIFLO (Lichtner and Seth 1996), and FLOTRAN (Lichtner 2000).
Geochemical transport models can be divided into several major groups, such as

41

Chapter 2models with specific chemistry, and more general, loosely coupled models often
involving sequential iteration (Yeh and Tripathi 1990; Steefel and MacQuarrie 1996; Leijet al. 1999; Carrayrou, Mose and Behra 2004). Models with specific chemistry are
generally limited in the number of species they can handle, while their application isusually restricted to problems having a prescribed chemical system. They are, however,much easier to use and computationally can be much more efficient than more generalmodels. Typical examples of models with specific chemistry are those simulating thetransport of major ions, such as UNSATCHEM (ŠimĤnek and Suarez 1994; 1997), or
various reclamation models. Models with more generalized chemistry generally invoke avariety of sequential iterative or non-iterative operator-splitting approaches in whichtransport and chemistry are solved in separate steps (Steefel and MacQuarrie 1996; Belland Binning 2004; Carrayrou, Mose and Behra 2004). Models of this type include the
DYNAMIX, HYDROGEOCHEM, MULTIFLO, and FLOTRAN models mentioned
earlier, as well as the OS3D/GIMRT code of Steefel and Yabusaki (1996). Such models
provide users with far more freedom in designing their particular chemical system, andhence permit a much broader array of applications.
When the source terms and the decay and biodegradation reactions are neglected,
multicomponent transport during one-dimensional transient variably-saturated flow may
be described with the following set of equations (e.g. ŠimĤnek and Suarez 1994):
wT c jwt

U

w(s j  p j )wt

wc j · wqc j
w§
¨T D
¸
wz ©
wz ¹ wz

j 1, 2,..., ns(5)where cj (ML-3), sj (MM-1) and pj (MM-1) are the total dissolved, sorbed/exchanged, andmineral (precipitated) concentrations, respectively, of aqueous component j, and ns is thenumber of aqueous components. The second term on the left-hand side of (5) is absent
for components that do not undergo ion exchange and precipitation/dissolution reactions.
The total concentration of a component j, defined as the sum of the dissolved, sorbed andmineral concentrations, is influenced only by transport processes which act on thesolution concentration cj, but not by chemical reactions. However, the relative fraction of
a component in each of the three phases (solution, sorbed, mineral) depends strongly on
the specific chemical processes of the system being considered. Therefore, Eq. (5) must
be augmented with equations describing the different equilibrium and non-equilibrium
chemical reactions such as complexation, cation exchange, and precipitation/dissolution(ŠimĤnek and Valocchi 2002).
The more general geochemical transport-modeling approach was implemented in ourrecent work (Jacques et al. 2003) on multicomponent transport in which we coupled the
HYDRUS-1D water flow and solute transport model (ŠimĤnek, Šejna and Van
Genuchten 1998) with the PHREEQC geochemical speciation model (Parkhurst andAppelo 1999). PHREEQC considers a variety of chemical reactions such as aqueous
speciation; gas, aqueous and mineral equilibrium; oxidation–reduction reactions; andsolid-solution, surface-complexation, ion-exchange and kinetic reactions. The HYDRUS1D and PHREEQC models were coupled using a non-iterative operator-splitting
approach (SNIA), which first solves the water flow and solute transport part of the
problem, and then the geochemical part. The combined HYDRUS1D-PHREEQC model(Jacques et al. 2003) permits simultaneous simulations of variably-saturated transientwater flow, multicomponent solute transport, and speciation and other geochemical
processes, including a broad range of mixed equilibrium and kinetic reactions.
42Van Genuchten and ŠimĤnek
The following two examples, taken from Jacques et al. (2003), illustrate the potential
power and versatility of the loosely coupled multicomponent geochemical modelingapproach used in HYDRUS1D-PHREEQC. The first example simulates heavy-metaltransport in a multi-layered soil profile assuming steady-state water flow and pHdependent cation exchange capacities. The second example extends the analysis tovariably-saturated flow by simulating the long-term fate and transport of heavy metals
under transient field conditions.
Example for steady-state flow
This example considers the transport of several major cations (Na, K, Ca, Mg) andthree heavy metals (Cd, Zn, Pb) in a 50-cm deep multi-layered soil profile having
different soil hydraulic properties and pH-dependent cation exchange capacities (CECs).
Assuming that the CEC is associated solely with organic matter, the cation exchange
capacity will increase significantly with increasing pH due to the acid-base properties of
its functional groups. This behavior is represented by a multi-site cation exchange
complex consisting of six sites, each having a different selectivity coefficient for theexchange of protons (Appelo, Verweij and Schafer 1998). The top 28 cm of the soil was
assumed to be contaminated with the three heavy metals (initial pH 8.5), while an acidmetal-free solution (pH 3) infiltrated into the soil. Flow was assumed to be steady-state ata relatively low constant flux density of 0.05 cm day-1, which caused the soil profile to beunsaturated (water content values varied between 0.37 and 0.15). The dispersivity anddiffusion coefficients were taken to be 5 cm and 0.80 cm2 day-1, respectively. Measured
soil hydraulic and other properties of the five soil layers are given elsewhere (Seuntjens
et al. 2002).
Results obtained with the combined HYDRUS1D-PHREEQC model will becompared in this example against results obtained with a different geochemical computer
program CRUNCH. The CRUNCH model is based on the GIMRT/OS3D package of
Steefel and Yabusaki (1996) and Steefel (2000). The geochemical reactions and transport
in CRUNCH are coupled in one of two ways: (1) a global implicit approach (GIMRT)
based on simultaneous solutions of the transport and reaction equations, or (2) an
operator time-splitting approach (OS3D) for transport and reactions that is quite similar
to the SNIA sequential non-iterative approach used in the combined HYDRUS1DPHREEQC model. GIMRT generally leads to smaller numerical errors than SNIA. A
comparison of HYDRUS1D-PHREEQC and CRUNCH-GIMRT allows one to assess
numerical discretization errors of the SNIA coupling as a function of the maximum
finite-element Courant number, Cr, which is defined as Cr= q't/T'x, where 't and 'x are
temporal and spatial discretizations of the invoked numerical solution, respectively.
Figure 2 (see Color pages elsewhere in this book) shows selected results for
simulations with CRUNCH-GIMRT, CRUNCH-OS3D (Cr = 0.5) and HYDRUS1DPHREEQC (Cr = 0.5 and 0.1). Infiltration of the low-pH solution causes an increase in
the number of protonated sites on the cation exchange complex (Figure 2c and d), leading
to leaching of the heavy metals. Cd concentrations in the leachate from the 50-cm deep
profile reach a peak after about 0.3 year (b), with most Cd leached from the profile after 1
year.
Results obtained with CRUNCH-SNIA and HYDRUS1D-PHREEQC using Cr = 0.5
showed very good agreement, especially for the outflow curves. The SNIA approach with
Cr = 0.5 produced slightly more numerical dispersion as compared to simulations using
the global implicit approach (CRUNCH-GIMRT), such as in the pH outflow curve, andthe different profiles after 0.7 year. Reducing the Courant number to 0.1 in HYDRUS1D-

43

Chapter 2
PHREEQC produced excellent agreement with CRUNCH-GIMRT.
Extension to transient variably-saturated flow
The coupled model HYDRUS1D-PHREEQC allows one to simulate reactive
transport during variably-saturated unsteady flow. Using the same soil profile data fromthe first example (but now assuming a fixed CEC), the transport of Na, K, Ca, Mg, Cd,
Zn and Pb was simulated for a period of 8.2 year using meteorological data for theKempen region in the northern part of Belgium from 1972 to 1981. Cumulative potential
and actual net fluxes (precipitation minus potential or actual evaporation) across the soil
surface are shown in Figure 3. The actual cumulative downward flux is higher than thepotential flux since the sandy soil could not deliver enough water to maintain thepotential evaporation rate during dry periods, thus leading to lower actual evaporation
rates. In addition to the atmospheric boundary conditions, steady-state flow simulations
were carried out with a constant surface flux density of 0.107 cm day-1.

Cumulative surface flux (cm)

300

200

100

0
00

01

02

03

04

05

Time (year)

06

07

08

Figure 3. Atmospheric boundary conditions used for the heavy-metal transport problem.
Shown are the cumulative potential (black line with dots) and actual (black line) net
surface water fluxes for the transient flow simulation, and the cumulative surface flux
(gray line) for the steady-state simulation
Figure 4 shows distributions of the total concentrations CT of Na, Ca and Cd in the
top 50 cm of the layered soil profile at four times. For the steady-state simulation, Ca andCd concentrations decreased in the top 20 cm. Cd is almost completely leached from thislayer due to aqueous complexation with Cl initially present in the soil solution. The useof atmospheric boundary conditions caused less leaching of elements from the top layer,in part due to upward flow and transport during dry periods. For example, the total
concentration of Na increased significantly after a long dry period at approximately 4.7
years. Because they are more strongly held on the cation exchange complex, Ca and Cd
show somewhat less pronounced accumulation. The effect of using atmosphericboundary conditions was relatively large; total concentrations were up to one order of44

Van Genuchten and ŠimĤnek
magnitude higher than those obtained assuming steady-state flow. This result is important
when, for example, simulating plant uptake of heavy metals or modeling the degradationof organic contaminants. This example shows that the coupling of HYDRUS-1D and
PHREEQC leads to a potentially very powerful tool for simulating a broad range ofinteracting physical, chemical and biological processes affecting solute transport in soils.
0

10

0

0

10

10

20

30

Depth (cm)Depth (cm)Depth (cm)

Na
20

Ca

30

20
10 days
2 year
4.66 years
8.22 years

30

40

40

40

50

50

50

Cd
1x10

-4

1x10

-3

1x10

CT (mol/l)

-2

1x10

-4

CT (mol/l)

1x10

-3

1x10

-7

1x10

-6

CT (mol/l)

1x10

-5

Figure 4. Distributions of the total concentration, cT (mol/l), of Na, Ca and Cd in the top
50 cm of a physically and chemically layered soil profile at four times during steady-state(gray lines) and atmospheric (precipitation/evaporation, black lines) flow conditions

Colloid and colloid-facilitated transport
Colloids and a variety of micro-organisms (including bacteria and viruses) are
subject to the same subsurface fate and transport processes as chemical compounds,
while additionally being subject to their own unique complexities. For example, many
colloids and microbes are negatively charged so that they are electrostatically repelled bynegatively-charged solid surfaces, which may lead to an anion exclusion process causingslightly enhanced transport relative to fluid flow. Size exclusion may similarly enhance
the advective transport of colloids by limiting their presence and mobility to the larger
pores (e.g. Bradford et al. 2003). At the same time, colloids and biologically reactive
solutes and bacteria may be subject to adsorption–desorption at solid surfaces, andaccumulation at air–water interfaces, although the exact physical-chemical processes ofsorption at air–water interfaces are still being debated (Thompson and Yates 1999; Wan
and Tokunaga 2002; Crist et al. 2004). In addition, the transport of colloids and microorganisms is affected by filtration and straining in the porous matrix, which is a functionof the size of the colloid or micro-organism, the water-filled pore size distribution, andthe pore water velocity (Elimelech and O' Melia 1990; Bradford et al. 2003; McCarthy
and McKay 2004). As such, colloid and microbial biomass accumulation may also reduce
the porosity of soils and alter the unsaturated soil hydraulic properties.
Additionally, inorganic, organic and microbiologically active colloids with diameters
between 0.01 and 10 micrometer may enhance the transport of otherwise immobile
chemicals or microbes by colloid-facilitated transport (Kretzschmar et al. 1999; Totscheand Kögel-Knabner 2004). Much evidence now exists that many contaminants aretransported not only in a dissolved state by water, but also sorbed to moving colloids.
Numerous examples exist in the literature illustrating this colloid-facilitated transport,
including for radionuclides (Von Gunten, Waber and Krähenbühl 1988; Noell et al.

45

Chapter 21998), pesticides (Vinten, Yaron and Nye 1983; Kan and Tomson 1990; Lindqvist andEnfield 1992), heavy metals (Grolimund et al. 1996), viruses, pharmaceuticals (Tolls
2001; Thiele-Bruhn 2003), hormones (Hanselman, Graetz and Wilkie 2003) and othercontaminants (Magee, Lion and Lemley 1991; Mansfeldt et al. 2004). Many of thesecontaminants are normally strongly sorbed to soil and hence should not be very mobile inthe subsurface. However, since they can also be sorbed to colloids which often move at
rates similar as non-sorbing tracers, the potential of enhanced transport of colloidadsorbed contaminants can be very significant.
In view of the many processes involved, numerical simulation of colloid and colloidfacilitated transport is still a major and largely unresolved challenge (De Jonge,
Kjaergaard and Moldrup 2004; DeNovio, Saiers and Ryan 2004). Models not only needto describe transport, attachment, detachment, straining and/or blocking or size exclusionof colloids, but also must consider transport, reactions, and kinetic and instantaneous
sorption of contaminants to both the soil and the colloids (Bradford et al. 2003; Totscheand Kögel-Knabner 2004). In this section we describe mathematical models for both
colloid transport and colloid-facilitated contaminant transport. The models have recently
been incorporated in the HYDRUS software packages.
Colloid transport
Colloid fate and transport models are commonly based on some form of theadvection–dispersion equation, but modified to account for colloid filtration (Harvey andGarabedian 1991; Hornberger, Mills and Herman 1992; Corapcioglu and Choi 1996). In
the absence of colloid inactivation and degradation, the colloid transport equation is thengiven bywT wCc w Aaw* c
wS
wCc · wqcCc
w §U c =¨ T w Dc
¸
wt
wt
wt
wx ©
wx ¹
wx

(6)

where Cc is the colloid concentration in the aqueous phase (nL-3), *c is the colloid
concentration adsorbed to the air–water interface (nL-2), Sc is the solid-phase colloid
concentration (nM-1), Tw is the volumetric water content accessible to colloids (L3L-3)(due to ion or size exclusion, Tw may be smaller than the total volumetric water content
T), Aaw is the air–water interfacial area per unit volume (L2L-3), Dc is the dispersion
coefficient for colloids (L2T-1), and qc is the volumetric water flux density for colloids(LT-1). The colloid mass-transfer term between the aqueous and solid phases, Esw, istraditionally given as:
E sw = UwSc= T wk ac\ s Cc - U k dcSc
wt

(7)

in which kac and kdc are first-order colloid attachment and detachment coefficients (T-1),
respectively, and \s is a dimensionless colloid retention function (-). The attachment
coefficient is generally calculated using filtration theory (Logan et al. 1995), a quasiempirical formulation in terms of the median grain diameter of the porous medium (often
termed the collector), the pore-water velocity, and collector and collision (or sticking)efficiencies accounting for colloid removal due to diffusion, interception andgravitational sedimentation (Rajagopalan and Tien 1976; Logan et al. 1995). The firstorder detachment coefficient in Eq. (7) accounts for colloid mobilization, presumably as
46

Van Genuchten and ŠimĤnek
affected by changes in pore-water chemistry (ionic strength, ionic composition, and pH)
and physical perturbations in flow, including changes in the flow rate and the watercontent. The attachment and detachment coefficients in (7) have been found to depend
strongly upon water content, with attachment significantly increasing as the watercontent decreases.
To simulate reductions in the attachment coefficient due to filling of favorable
sorption sites, \s is sometimes assumed to decrease with increasing colloid massretention. Random sequential adsorption (Johnson and Elimelech 1995) and Langmuirian
dynamics (Adamczyk et al. 1994) equations have been proposed for \s to describe thisblocking phenomenon, with the latter equation given by:

1

\s

Sc
Scmax

(8)

in which Scmax is the maximum solid-phase colloid concentration (nM-1). Alternatively,
blocking could possibly be described also using non-linear Freundlich sorption.
Conversely, enhanced colloid retention during porous-medium ripening can theoreticallybe described using a functional form of \s that increases with increasing mass of retained
colloids. To our knowledge, no functional forms for \s have been proposed to describe
porous-medium ripening. We refer to several recent studies (Ginn et al. 2002; DeNovio,
Saiers and Ryan 2004; Rockhold, Yarwood and Selker 2004) for more detailed
discussions of the attachment and detachment coefficients in Eq. (7).
Notice that Eq. (7) lumps the effects of a variety of physical and chemical processes
into a single attachment-coefficient parameter. Bradford et al. (2002; 2003) hypothesizedthat the influence of straining and attachment mechanisms on colloid retention should be
separated into two distinct components: attachment and detachment per se, and straining(being the entrapment of colloids in pore throats that are too small to allow passage).
They modeled the influence of straining using an irreversible first-order expression. In
that case, Esw becomes:
wSc
E sw = Uwt

wScstrwScatt=U
E +E = U
wt
wt
str
att
T w\ s Cc kstr  (T w\ s Cc k ac - U k dcSc )
str
sw

att
sw

(9)

where Eswstr and Eswatt are the rates of mass exchange for colloid straining and attachment/
detachment (nL-3T-1), respectively, kstr is the first-order straining coefficient (T-1), and
Scstr and Scatt are the solid-phase concentrations of strained and attached colloids (nM-1),
respectively. The first term on the right-hand side of the above equation now accounts forstraining and the second term for attachment–detachment.
Application of the first-order attachment-–detachment model given by Eq. (7)typically leads to exponential colloid distributions versus depth. Bradford et al. (2003)
showed that such exponential distributions are often inconsistent with experimental data.
They obtained much better results using a depth-dependent straining coefficient in Eq.
(9) of the form

\

str
s

§d z·= ¨ 50¸© d 50 ¹

-E

(10)
47

Chapter 2

where d50 is the median grain size of the porous media (L), E is a fitting parameter (-),
and z is distance from the porous-medium inlet (L). Data from Bradford (2002; 2003) for
different colloid diameters (dp) and porous media median grain sizes showed an optimalvalue of 0.43 for E while kstr could be described using a unique increasing function of theratio of dp/d50. Subsequent studies with layered soils (Bradford et al. 2004) showed theimportance of straining at and close to textural interfaces when flow occurs in the
direction of coarse-textured to medium- and fine-textured media. Straining was found to
be a significant mechanism for colloid retention for values of dp/d50 greater than 0.005.
That study also showed the importance of liquid-phase velocity distributions, colloid
retention at soil textural interfaces, and size exclusion.
Finally we note that a model similar to Eq. (7) may be used to describe thepartitioning of colloids to the air–water interface:
E aw =

w Aaw* c= T w k aca\ a Cc - Aawk dca ī c
wt

(11)

where Eaw is the colloid mass-transfer term between the bulk water and the air–waterinterface (nL-3T-1), \a is a dimensionless colloid retention function for the air–waterinterface (-) similarly as used in (7), and kaca and kdca are the first-order colloid
1
attachment and detachment coefficients to/from the air–water interface ( T ),respectively.
Colloid-facilitated solute transport
Colloid-facilitated transport requires knowledge of colloid transport, dissolvedcontaminant transport, and colloid-facilitated contaminant transport. Transport and/or
mass-balance equations must therefore be formulated for the total contaminant, forcontaminant sorbed kinetically or instantaneously to the solid phase, and for contaminantsorbed to mobile colloids, to colloids attached to the soil solid phase, and to colloids
accumulating at the air–water interface. As an illustration of the complexities involved,
we review here the equations for colloid-facilitated transport that we recently
incorporated in the HYDRUS software packages.

Mass-balance equation for the total contaminant. The combined dissolved and colloidfacilitated contaminant transport equation (in one dimension) is given by:
wS
wS wT C S
wS S
w Aaw* c S ac
wT C
wJ U e  U k  w c mc  U c ic 
- T R
wt
wt
wt
wt
wt
wt
wx
wC · wq C S
wJw §wC · wqC w §- T ¨T w S mc Dc c ¸  c c mc
¨T D¸wx wx ©
wx ¹ wx wx ©
wx ¹
wx

(12)

where T is the volumetric water content (L3L-3) (note that we use the entire water content
for the contaminant ), C is the dissolved-contaminant concentration in the aqueous phase(ML-3), Se and Sk are contaminant concentrations sorbed instantaneously and kinetically,
respectively, to the solid phase (MM-1); Smc, Sic, and Sac are contaminant concentrations
sorbed to mobile and immobile (attached to solid and air–water interface) colloids (Mn1), respectively; JT is the total flux of the dissolved contaminant (ML-3T-1), D is thedispersion coefficient for contaminants in solution (L2T-1), q is the volumetric water flux
48

Van Genuchten and ŠimĤnek
density for the contaminant (LT-1), and R represents various chemical and biologicalreactions, such as degradation and production (ML-3T-1), discussed below. Note that theleft-hand side sums the mass of contaminant associated with the different phases(contaminant in the liquid phase, contaminant sorbed instantaneously and kinetically tothe solid phase, and contaminant sorbed to mobile and immobile (attached to solid phaseor air–water interface) colloids), while the right-hand side considers various spatial mass
fluxes (dispersion and advective transport of the dissolved contaminant, and dispersion
and advective transport of contaminant sorbed to mobile colloids).
Mass-balance equation for contaminant sorbed to the solid phase. Equation (13) invokes
the concept of two-site sorption for modeling non-equilibrium adsorption–desorption
reactions (e.g. Van Genuchten and Wagenet 1989). The two-site sorption concept
assumes that total sorption, S, can be divided into two fractions:S

Se  Sk(13)

with sorption Se (MM-1) on one fraction of the sites (type-1 sites) assumed to beinstantaneous, and sorption Sk (MM-1) on the remaining sites (type-2 sites) being timedependent according to

U

wSk
wt

T kas C  U kds Sk  Rs

(14)

where kas is the rate of the contaminant sorption to the solid phase (type 2) (T-1), kds is therate of the contaminant desorption from the solid phase (type 2) (T-1), and Rs represents
various chemical and biological reactions of the kinetically sorbed contaminant (ML-3T1).
Mass-balance equation for contaminant sorbed to mobile colloids. The mass-balance
equation for contaminant sorbed to mobile colloids can be written as
wT wCc S mc
wt

wC · wq C S
w §T w S mc Dc c ¸  c c mc 
¨
wx ©
wx ¹
wx

T kamc\ mC  T w kdmc Cc Smc  T w (kac\ s  kstr\ sstr )Cc Smc  U kdc Scatt Sic  Rmc(15)

where kamc is the adsorption rate to mobile colloids (T-1), kdmc is the desorption rate frommobile colloids (T-1), and Rmc represents various chemical and biological reactions forcontaminant sorbed to mobile colloids (ML-3T-1). The parameter \m adjusts the sorption
rate to the number of mobile colloids present, i.e.,

\m

Cc
Ccref

(16)

where Ccref is the reference concentration of colloids for which the sorption rate kamc isvalid (ML-3). In equation (15), the first two terms on the right-hand side represent
dispersion and advective transport, respectively, of contaminant sorbed to mobile
49

Chapter 2
colloids; the third and fourth terms account for sorption and desorption of contaminants
to/from mobile colloids, respectively; the fifth and sixth terms account for the attachment(including straining) and detachment of mobile colloids containing sorbed contaminants,
respectively; while the seventh term represents degradation or other reactions involving
contaminant sorbed to mobile colloids.
Mass-balance equation for contaminant sorbed to immobile colloids. The mass-balance
equation for contaminant sorbed to immobile colloids can be written as follows

U

wSc Sic
wt

T kaic\ i C  U kdic Sc Sic  T w (kac\ s  kstr\ sstr )Cc Smc  U kdc Scatt Sic  Ric(17)

where kaic is the adsorption rate to immobile colloids (T-1), kdic is the desorption rate from
immobile colloids (T-1), and Ric represents various reactions for contaminant sorbed toimmobile colloids (ML-3T-1). The parameter \i adjusts the sorption rate to the number ofimmobile colloids present:

\i

Sc
Scref

(18)

where Scref is the reference concentration of immobile colloids for which the sorption ratekaic is valid (ML-3). In equation (17) the first two terms on the right-hand side represent
adsorption and desorption of contaminant to/from immobile colloids, respectively; thethird and fourth terms describe the attachment (including straining) and detachment ofimmobile colloids with sorbed contaminant, respectively; and the fifth term represents
reactions of contaminant sorbed to immobile colloids.
Mass-balance equation for contaminant sorbed to colloids attached to the air–waterinterface. The mass-balance equation for contaminant sorbed to colloids attached tothe air–water interface may be written as
w Aaw* c S ac
wt

T kaac\ g C  Aawkdac * c Sac  T wk aca\ a Cc Smc - Aawk dca ī c Sac  Rac(19)

where kaac is the adsorption rate to colloids at the air–water interface (T-1), kdac is thedesorption rate from colloids at the air–water interface (T-1), and Rac represents various
reactions for contaminant sorbed to colloids at the air–water interface (ML-3T-1).
Similarly as in (17), the parameter \g adjusts the sorption rate to the number of colloids
at the air–water interface:\g*c* cref(20)

where *cref is the reference concentration of immobile colloids for which sorption rate kaic
is valid (ML-3). In equation (19) the first two terms on the right-hand side represent thesorption and desorption, respectively, of contaminant to/from colloids at the air–waterinterface; the third and fourth terms account for the attachment and detachment,

50

Van Genuchten and ŠimĤnek
respectively, of colloids with sorbed contaminant to/from the air–water interface;whereas the fifth term represents degradation and other reactions of contaminant sorbed
to colloids accumulated at the air–water interface.
Reaction term. The reaction term R in (12) may be used to account for a variety ofchemical and biological reactions and transformations, including degradation andproduction, not already explicitly incorporated in the main total contaminant masstransport equation. Consistent with current capabilities of the HYDRUS software
packages (ŠimĤnek, Šejna and Van Genuchten 1998; 1999) to simulate sequential firstorder decay chains, R may include provisions for two first-order degradation reactions:
one which is independent of other solutes and one which provides the coupling between
solutes involved in sequential first-order decay reactions. As discussed earlier in Section
2, problems of solute transport involving sequential first-order decay reactions frequently
occur in soil and groundwater systems. The reaction term R in (12) for colloid-facilitatedtransport scenarios is now given by:R-( P w  P w' )T C - ( P s  P s' ) Ub ( Se  Sk ) - ( P c  Pc' )(T w Cc S mc  U Sc Sic 
'*

*

'*

*

*

'*

*

*

*

A aw* c Sac )  P wT C  P s U ( S e  S k )  Pc (T w Cc S mc  U Sc Sic  A aw* c S ac )(21)

where Pw, Ps , and Pc are first-order rate constants (T-1) for solutes in the liquid, solid andcolloid phases (T-1), respectively; Pw', Ps', and Pc' are similar first-order rate constants (T1) providing connections between individual chain species, and terms with the superscript* belong to solutes preceding in the chain reaction. The reaction terms Rs, Rm, Rim, and Ra
for reactions in the kinetically sorbed phase, on mobile colloids and on colloids
associated either with the solid phase or the air–water interface, respectively, are as
follows:
Rs

( P s  P s' ) U S k  P s'* U S k*

Rmc*
( P s  P s' )T wCc S mc  Pc'*T wCc S mc

Ric

( P c  Pc' ) U Sc Sic  Pc'* U Sc Sic*

Rac*
( P c  Pc' ) Aaw* c S ac  Pc'* Aaw* c S ac

(22)

The above mathematical development shows that a complete description of colloidfacilitated transport requires a total of nine coupled partial differential equations
involving nine unknown variables (Cc, Scstr, Scatt,*c, C, Sk, Smc, Sic, Sac).
Example application for colloid-facilitated transport
Typical features of colloid and colloid-facilitated transport are demonstrated here fora hypothetical column experiment. The column was assumed to have a length, L, of 10
cm, while the experiment lasted 600 minutes. Water flowed through the column at fullsaturation (the saturated water content was equal to 0.50) at a flux density, q, of 0.1 cm
min-1. Both colloids and the contaminant were assumed to be applied at unit (relative)
concentrations at the top of the column during a time period of 60 min. The soil bulk
density, Uwas set equal to 1.5 g cm-3, while the dispersivity, O, for both the colloids andthe contaminant was assumed to be 0.1 cm. The colloid attachment and detachment
coefficients, ka and kd, were taken to be 0.01 and 0.005 min-1, respectively. Solute

51

Chapter 2
sorption to soil was considered instantaneous with the distribution coefficient, Kd, equal
to 2 L kg-1. Solute adsorption was assumed to be the same for both the mobile andimmobile colloids, but with different rate coefficients of kamc= kaic=0.1 min-1 and kdmc=
kdic=0.02 min-1, respectively, and with Ccref= Scref=1. For these conditions, one porevolume, T, is equal to 50 minutes [t = TL/v = TLTs/q = 1 * (10 cm) * 0.5 / (0.1 cm/min) =50 min], while the retardation factor for the contaminant equals 7 (R = 1 + UKd/Ts = 1 +1.5*2/0.5 = 7).
Figure 5 shows colloid and total solute concentrations at depths of 5 and 10 cm. The
main concentration fronts for both colloids and solute arrived, as expected, at 1 and 7
pore volumes (i.e., after 50 and 350 minutes), respectively. However, a significant
amount of solute arrived much earlier than at 7 pore volumes. This earlier arrival is dueto the accelerated movement of solute sorbed to mobile colloids.
Figure 6 presents colloid and solute fluxes at the bottom of the column. Notice thatthe solute flux (Figure 6b) has two concentration peaks. The first peak corresponds tocontaminant arriving sorbed to colloids, which have a retardation factor, R, equal to 1 (nosorption or anion exclusion), while the second peak corresponds to solute arriving
dissolved in water (R=7, t=350 min). Also notice that the initial contaminant fluxes(Figure 6b) are smaller relative than the colloid fluxes, which show a plateau of certain
duration. This is because contaminant sorption onto the colloids was assumed to be a
kinetic process which requires a certain time period for sorption to be complete. The
contaminant flux for this reason keeps increasing until it suddenly drops at the end of thecolloid pulse. Hence, the largest contaminant fluxes during the first peak occurred justbefore the end of the colloid pulse (Figure 6b). Colloids keep arriving at the bottom of
the column after the main colloid pulse due to kinetic colloid detachment from the solid
phase (Figure 6a). By comparison, considerable solute fluxes are present between the twosolute peaks as a result of solute both sorbed to colloids and, increasingly, dissolved inwater. The dissolved solute fluxes at the end of the column were initially also due toenhanced transport by the colloids in especially the upper part of the columns, andsubsequent desorption into the liquid phase (Figure 6b). As expected, the bulk of thecontaminant arrived at approximately 350 minutes, consistent with a solute retardation
factor, R, of 7.
0.4

1

a)

b)

5 cm
10 cm

5 cm
10 cm

0.8

Concentration [-]Concentration [-]

0.3

0.6

0.4

0.2

0.1
0.2

0

0
0

100

200

300

Time [min]

400

500

600

0

100

200

300

400

500

600

Time [min]Figure 5. Colloid (a) and contaminant (b) concentrations at a depth of 5 and 10 cm

52

Van Genuchten and ŠimĤnek
0.08

0.03

a)

b)
0.025

0.06

Solute Flux [-/min]Solute Flux [-/min]

0.07

0.05
0.04
0.03

0.02

0.015

0.01

0.02
0.005

0.01
0

0
0

100

200

300

400

500

600

Time [min]

0

100

200

300

400

500

600

Time [min]Figure 6. Calculated colloid (a) and contaminant (b) flux densities at the bottom of the
10-cm-long soil column

Coupled surface – subsurface flow
Much progress has been made over the years also in modeling various componentsof overland flow and variably-saturated subsurface flow. Traditionally, however, surface
and subsurface flow processes have been considered mostly separately, with some of theprocesses being severely simplified when coupled surface and subsurface flow wassimulated. For example, existing surface irrigation models typically still use empiricalinfiltration functions such as the Philip, Kostiakov, modified Kostiakov, and Branch
infiltration equations (e.g. Woolhiser, Smith and Goodrich 1990; Strelkoff, Clemmens
and Schmidt 1998), rather than rigorously simulating subsurface water flow and solutetransport using the Richards equation for variably-saturated flow, and advection–
dispersion equations for solute transport. Several weakly coupled, mostly iterative ortime-lagged approaches using the Richards equation have been pursued for larger-scalesurface runoff or small watershed-scale applications. Early examples are given by Zhang
and Cundy (1989) and Govindaraju and Kavvas (1991). Still, because of a lack of tight(numerically implicit) coupling between the surface and subsurface flow and transport
processes, previous models have not proven reliable for untested management practices
under surface-irrigated conditions, nor have they been applied to linked scenarios
involving surface irrigation, chemical runoff, deep drainage, solute leaching and crop
growth.
Recently, several more integrated approaches have been proposed. The first andprobably still the most comprehensive integrated physically-based model was formulatedby VanderKwaak (1999). This model included most or all known stream-flow generation
mechanisms, variably-saturated subsurface flow, advective–dispersive solute transport inboth surface and subsurface flow, and provisions for preferential flow within vadosezone fractures and macropores. The surface and subsurface components were linked
using first-order flux relationships such that one set of discrete algebraic expressions
resulted from the fully-coupled numerical approach. The model was applied to a small
watershed in Canada (VanderKwaak 1999), as well as to a well-studied catchment in
Oklahoma (VanderKwaak and Loague 2001). Panday and Huyakorn (2004) proposed asimilar fully-coupled spatially-distributed model for conjunctive surface/subsurface flow
involving three-dimensional variably-saturated subsurface flow, two-dimensional
overland flow, and one-dimensional flow through such surface structures as rivers andcanals, among other processes. Another recent example of fully coupled variably-53

Chapter 2
saturated flow and surface runoff is given by Beaugendre et al. (in prep.).
The two-dimensional illustrative example application considered here pertains tohill-slope hydrology. Water flow and solute transport in and along hill slopes is arelatively complex nonlinear problem. Rainfall water will infiltrate into the soil profile at
a rate equal to the rainfall rate until the soil infiltration capacity is reached. Surface
runoff will be generated once the soil infiltration capacity is exceeded. This surface
runoff will redistribute water along the land surface by moving it to lower parts of the hillslope where it can infiltrate if locally enough infiltration capacity is present. As morewater infiltrates in the lower part of the hill slope, more water will be available there forvegetation. More vegetation will generally result in more roots which, in turn, may enrichsoil with organic matter, thus improving soil structure and further increasing thehydraulic conductivity and the infiltration capacity (Mattson et al. 2004). We are notaware of existing vadose-zone flow/transport models that describe this dynamic
interaction between overland flow, subsurface infiltration and plant growth.
Governing equations
Hortonian overland flow is usually described using kinematic wave equations, which
are simplifications of the Saint Venant equations and provide excellent approximations
for most overland flow conditions (Woolhiser and Liggett 1967; Morris and Woolhiser
1980; Woolhiser, Smith and Goodrich 1990):

wh wQ

wt wx

q ( x, t )(23)

where h is the unit storage of water (or mean depth for smooth surfaces) (L), Q is thedischarge per unit width (L2T-1), t is time (T), x is the distance co-ordinate over the soil
surface (L), and q(x,t) is the rate of local input, or lateral inflow (i.e., local precipitation
minus local infiltration) (LT-1).
The discharge Q per unit width can be calculated as follows:
Q D hm(24)

where D (L2-mT-1) and m (-) are parameters related to slope, surface roughness and flow
conditions (laminar or turbulent flow). The parameter D is usually evaluated using theManning hydraulic resistance law

Dk0

Sm1/ 2
n

and

m 5/3

(25)

1
where Sm is the slope [-], k0 is an empirical constant (L1/3T-1) equal to 1.49 cm1/3 s , andn is Manning’s roughness coefficient for overland flow (-).
The rate of local input q(x,t) in (23) integrates various inputs and outputs such as
precipitation, evaporation and infiltration. Models that rigorously evaluate thesecomponents, such as the HYDRUS codes, should be ideal tools for calculating q(x,t). By
substituting (24) into (23) we obtain:

54

Van Genuchten and ŠimĤnek
wh wD h mwtwx

q ( x, t )

(26)

or
wh
wh
 D mh m 1wt
wx

q ( x, t )(27)

The above equations were solved using a numerically stable fully implicit four-point
finite-difference method, similar to the one used in the KINEROS model (Woolhiser,
Smith and Goodrich 1990). This numerical solution of the overland equation was
subsequently coupled to the HYDRUS-2D computation module (ŠimĤnek, Šejna and
Van Genuchten 1999) using transient atmospheric boundary conditions.
Example application
We used the tightly coupled overland and subsurface flow approach as incorporated
into HYDRUS-2D to simulate surface runoff from a rainfall event of limited duration
over a 100-m heterogeneous hill slope. The example considers surface runoff generated
by a 10-minute high-intensity rainstorm (0.00667 cm s-1 or 24 cm h-1) covering the entire
hill slope. The hill slope is assumed to consist of two soil materials. The soil in the
middle third of the hillside transect (between 33 and 66 m) is assumed to have a saturated
hydraulic conductivity (Ks= 0.0289 cm s-1 or 25 m d-1) that is two orders of magnitude
higher than that of the other parts of the hill slope (Ks= 0.000289 cm s-1 or 25 cm d-1).
Because of this, the middle section can accommodate the infiltration of water from bothrainfall itself and from runoff coming from the upgradient part of the hill slope. The soiltransect has a slope of 0.01, while the roughness coefficient n was assumed to be equal to
0.01.
Figure 7 shows the depth of the water layer that develops on top of the soil surface(see Color pages elsewhere in this book). The figure also shows the steady-state waterlayer calculated using the analytical solution of Eq. (26) assuming no infiltration (i.e.,q=q0 = 24 cm h-1):§ nq x ·h( x) ¨ 01/ 2 ¸© k0 S m ¹

3/ 5

(28)

The figure shows both the buildup and recession of the water layer along the hill slope atvarious times. Actual and cumulative surface runoff fluxes from the bottom of the hillslope transect, the cumulative effective rain on the transect (rainfall minus infiltration),and the volume of water in the surface layer for a 10-minute rain storm are shown in
Figure 8 (see Color pages elsewhere in this book). Notice that a nearly steady-state
situation for overland flow had developed at the end of the rainfall event.
Finally, Figure 9 shows contours of the water content in the soil profile 6 min after
initiation of the rainfall event (see Color pages elsewhere in this book). Figure 9, as wellas Figure 7, clearly shows that runoff generated in the upper third of the hillside transect
moves to and infiltrates in the upper part of the middle section containing more
permeable soil, and that runoff is generated again in the lower one third of the hill slope.
The example demonstrates the potential of tightly coupled surface/subsurface flow
models.

55

Chapter 2

Preferential flow
A major challenge in dealing with the vadose zone, both in terms of modeling andexperimentation, is its overwhelming heterogeneity. One manifestation of heterogeneity
at intermediate spatial scales involves the preferential movement of water and chemicals
through soil macropores or rock fractures. In structured or macroporous soils, water may
move preferentially through large interaggregate pores, decayed root channels, andearthworm burrows, as well as through drying cracks in fine-textured soils, thereby
bypassing much of the soil matrix. Preferential flow may also occur in the form of
unstable flow (fingering) induced by soil textural layering, water repellency and/or air
entrapment. In unsaturated fractured rock, water may similarly move preferentiallythrough fractures and fissures, thus bypassing much of the rock matrix. Preferential flow
may be caused also by funneling of water through high-conductivity layers, and/or bybeing redirected by sloping less-permeable layers. We refer to several recent reviews andother papers for detailed discussions of the various processes and conditions leading topreferential flow (Ritsema et al. 1993; De Rooij 2000; Evans, Nicholson and Rasmusson
2001; National Research Council 2001; Bodvarsson, Ho and Robinson 2003; ŠimĤnek et
al. 2003; Wang et al. 2004).
An impressive array of models have been developed over the years in attempts toparameterize preferential-flow processes. While not all-inclusive, Figure 10 provides auseful schematic of increasingly complex models that have been used to simulate
preferential-flow processes. Figure 10a represents the traditional case of uniform(equilibrium) flow and transport, assuming applicability of an equivalent-continuum
approach based on the Richards equation (Eq. 1) for variably-saturated flow and Eq. (3)for advective–dispersive solute transport. The simplest situation for apparent preferential
flow arises when the Richards and advection–dispersion equations are still used in an
equivalent matrix and fracture-continuum approach, but now with composite hydraulicconductivity (permeability) curves, K(h), of the type shown in Figure 10b. While stillleading to uniform flow, models using such composite hydraulic properties do allow for
faster flow and transport during conditions near saturation, and as such may provide more
realistic simulations of field data (e.g. Peters and Klavetter 1988; Mohanty et al. 1997;
1998; Zurmühl and Durner 1996; De Vos et al. 1999).
The two parts of the conductivity curves in Figures 10b and 11a may be associatedwith soil structure (near saturation) and soil texture (at lower negative pressure heads).
This conceptualization of K(h) is consistent with analyses by Schaap and Leij (2000) of
the UNSODA soil hydraulic data base (Leij et al. 1996), which show that the ratio of themeasured (Ks) and extrapolated (Kms) values of the hydraulic conductivity is about oneorder of magnitude (Kms would be the value of Ks when no macropores or fractures were
present). We note that the use of composite hydraulic functions that lump the effects ofmatrix and fracture flow into one equation, in conjunction with the Richards equation,
will still lead to uniform flow, and as such cannot reproduce non-uniform moisture
distributions typical of preferential flow. Still deconvolution of bimodal conductivity
functions like those shown in Figure 11a may provide important information for dualpermeability models that consider separate flow domains for the fractures and matrix.

56

Van Genuchten and ŠimĤnek

Figure 10. Alternative conceptual models for flow through variably-saturated
structured media (after Altman et al. 1996)
Non-uniform flow can be generated numerically only using dual-porosity or dualpermeability models for structured media, or using alternative formulations specificallyderived for unstable flow (De Rooij 2000; Eliassi and Glass 2002; Dautov et al. 2002). Dualporosity and dual-permeability models typically assume that the medium consists of two
interacting pore regions, one associated with the macropore or fracture network, and one
with the micropores inside soil aggregates or rock matrix blocks. Different formulations arise
depending upon how water and solute movement in the micropore region is modeled, andhow water and solutes between the micropore and macropore regions are allowed to interact.

57

Chapter 2

Figure 11. Plots of (a) the field-averaged hydraulic-conductivity function (geometricmeans and standard errors, left) and (b) the field-averaged soil water retention function(arithmetic means and standard errors) for the surface horizon at the Las Nutrias
experimental site (after Mohanty et al. 1997)
Dual-porosity type models are obtained when the liquid phase is partitioned into mobile(fracture) and immobile (matrix) liquid pore regions, with water and/or solutes allowed to
exchange between the two liquid regions (Figure 10c). Such models can be made quite
general by permitting transient variably-saturated flow in the fracture, andsimultaneously allowing water to exchange between the fracture and matrix domains(ŠimĤnek et al. 2003). The latter dual-porosity situation would lead to both advective anddiffusive exchange of solutes between the fracture and matrix regions, but still without
vertical flow in the matrix region. Popular early conceptualizations of the dual-porosity
approach, when applied to solute transport only, are two-region or mobile–immobile
water models (e.g. Van Genuchten and Wierenga 1976), in which solute exchange is
described using first-order mass-transfer equations.
More complex dual-permeability models (Figure 10d) arise when water flow occurs
in both the fracture and matrix domains. Examples of dual-permeability models are given
by Gerke and van Genuchten (1993; 1996), Pruess (1991), Jarvis (1998), Liu, Doughty
and Bodvarsson (1998) and Liu, Zhang and Bodvarsson (2003), among many others.
These models all use different formulations for the exchange of water between thefracture and matrix regions. A summary of various exchange (mass transfer) terms fordual-porosity and dual-permeability models is given by ŠimĤnek et al. (2003). We notethat in some dual-permeability models (e.g. Hutson and Wagenet 1995; Wilson, Jardine
and Gwo 1992) more than two domains are considered, all having their unique hydraulicproperties.
The modeling approach can be further refined by considering transient variablysaturated flow and/or transport in discrete well-defined macropores or fractures, either

58

Van Genuchten and ŠimĤnek
without (Figure 10e) or with (Figure 10f) interactions between the fracture and matrix
domains. The approach typically assumes that the flow and transport equations of themacropore or fracture network of prescribed geometry can be solved simultaneously andin a fully-coupled fashion with the corresponding equations for the porous matrix.
Discrete-fracture models of this type (Figure 10f) are given by Shikaze, Sudicky andMendoza (1994) for two-dimensional gas-phase flow and transport through a network ofvadose-zone fractures embedded in a variably-saturated matrix, and by Therrien and
Sudicky (1996) and VanderKwaak (1999) for more general three-dimensional conditions
in which the Richards equation is applied both along a network of interconnected fracture
planes and in the adjoining porous matrix. Discrete fracture network models are likely
unfeasible for most or practical applications problems because of limited data availability
and computational difficulties (Liu, Zhang and Bodvarsson 2003).
Except for the discrete fracture models, all of the models schematically shown in
Figure 10 were recently incorporated in the HYDRUS-1D and HYDRUS-2D software
packages (ŠimĤnek et al. 2003), including the use of composite soil water retention andhydraulic conductivity functions based on the formulation of Durner (1994).
Additionally, we incorporated into the HYDRUS codes a dual-permeability approach
using the kinematic wave equation for flow in the macropores (Jarvis 1998), as well as arelatively simple empirical but still very effective single-porosity non-equilibrium model
proposed by Ross and Smettem (2000) to account for time-dependent water-content
equilibration at a given pressure head [a related more process-based approach involving
dynamic pressure-head equilibration is given by Hassanizadeh, Celia and Dahle (2002)].
Because of the large number of features now included in the updated HYDRUS models,
we believe that the resulting codes provide great flexibility in addressing a large numberof practical field-scale vadose-zone flow and transport problems involving both uniform
and preferential flow (ŠimĤnek et al. 2003; Zhang, ŠimĤnek and Bowman 2004).

Concluding remarks
The topics reviewed in this paper reflect the tremendous advances that have been
made during the past several decades in our ability to describe water and solute transport
processes in the subsurface mathematically. Within the context of our own work, we
reviewed relatively standard approaches for modeling flow and transport, andadditionally ongoing research on (1) multicomponent geochemical transport, (2) colloid
and colloid-facilitated transport, (3) integrated surface/subsurface modeling, and (4)process-based descriptions of preferential flow.
Continued progress in the above four areas of research requires significant advancesin both numerical modeling and the underlying science. Of the four topics, we believethat the main challenges with multicomponent transport modeling and integrated
surface/subsurface flow/transport modeling are first and foremost numerical in terms ofcoupling processes that have for too long been addressed in separate efforts. Still,
formidable challenges remain here. For example, numerical algorithms and databases for
multicomponent transport models must be extended to higher temperatures and ionicstrengths, complex contaminant mixtures (including especially mixed organic andinorganic wastes), multiphase flow, redox disequilibria for low-temperature systems, andcoupled physico-chemical systems to account for possible changes in the soil water
retention and hydraulic-conductivity functions. Integrated surface/subsurface numericalmodels require further research on such issues as spatial and temporal scaling ofhydrological, chemical and biological processes and properties, linking constitutive (soil-

59

Chapter 2
hydraulic) relationships to the measurements scale, preferential flow, and issues ofparameter and model uncertainty. Many of these and related specific challenges invadose-zone research are discussed as part of a recent Department of Energy National
Roadmap for Vadose Zone Science and Technology (Stephens et al. 2002; U.S.
Department of Energy 2001).
By comparison, we believe that the basic scientific issues related to colloid andcolloid-facilitated transport are still largely unresolved, and that our understanding herelags far behind current numerical capabilities. Much work is needed to better understandthe processes of filtration, straining, size exclusion, and mobilization of colloids andmicro-organisms; accumulation at air–water interfaces, interactions between microorganisms and contaminants (including biodegradation), the effects of both physicalfactors (water content, flow velocity, textural interfaces) and chemical processes (ionicstrength, solution composition, pH) on colloid retention and mobilization, and modelingcolloid-facilitated transport during conditions of transient flow. Addressing preferential
flow phenomena, and the related more general problems of subsurface heterogeneity,
poses equally important challenges, as well as enormous opportunities, since those
problems are at the center of many of the unresolved issues in other areas of vadose-zone
flow and transport research.

Acknowledgments
We like to acknowledge the very significant input of Diederik Jacques and Scott
Bradford to research documented in the sections on multicomponent geochemicaltransport and on colloid and colloid-facilitated transport, respectively. This paper is basedon work supported in part by SAHRA (Sustainability of semi-Arid Hydrology andRiparian Areas) under the STC Program of the National Science Foundation, Agreement
No. EAR-9876800.

References
Abrams, R.H., Loague, K. and Kent, D.B., 1998. Development and testing of acompartmentalized reaction network model for redox zones in contaminatedaquifers. Water Resources Research, 34 (6), 1531-1541.
Adamczyk, Z., Siwek, B., Zembala, M., et al., 1994. Kinetics of localized adsorption ofcolloid particles. Advances in Colloid and Interface Science, 48, 151-280.
Ahuja, L.R., Rojas, K.W., Hanson, J.D., et al., 1999. Root zone water quality mode:
modeling management effects on water quality and crop production. Water
Resources Publications, Highlands Ranch.
Altman, S.J., Arnold, B.W., Ho, C.K., et al., 1996. Flow calculations for Yucca
Mountain groundwater travel time (GWTT-95). Sandia National Laboratories,
Albuquerque. Report no. SAND96-0819.
Appelo, C.A.J., Verweij, E. and Schafer, H., 1998. A hydrogeochemical transport model
for an oxidation experiment with pyrite/calcite/exchangers/organic matter
containing sand. Applied Geochemistry, 13 (2), 257-268.
Beaugendre, H., Ern, E., Escaffer, T., et al., in prep. Numerical investigation of surface
runoff in hillslopes with variably-saturated flows. In: Neittaanmäki, P., Rossi, T.,
Majava, K., et al. eds. Proceedings of the European congress on computational
methods in applied sciences and engineering, ECCOMAS 2004, Jyväskylä, July
24-28, 2004.

60

Van Genuchten and ŠimĤnek
Bell, L.S.J. and Binning, P.J., 2004. A split operator approach to reactive transport withthe forward particle tracking Eulerian-Lagrangian localized adjoint method.
Advances in Water Resources, 27 (4), 323-334.
Bodvarsson, G.S., Ho, C.K. and Robinson, B.A. (eds.), 2003. Yucca Mountain Project.
Journal of Contaminant Hydrology, 62/63(special issue).
Bradford, S.A., Bettahar, M., ŠimĤnek, J., et al., 2004. Straining and attachment of
colloids in physically heterogeneous porous media. Vadose Zone Journal, 3 (2),384-394. [http://vzj.scijournals.org/cgi/content/full/3/2/384]
Bradford, S.A., ŠimĤnek, J., Bettahar, M., et al., 2003. Modeling colloid attachment,straining, and exclusion in saturated porous media. Environmental Science andTechnology, 37 (10), 2242-2250.
Bradford, S.A., Yates, S.R., Bettahar, M., et al., 2002. Physical factors affecting thetransport and fate of colloids in saturated porous media. Water Resources
Research, 38 (12), 631-6312.
Bromilow, R.H. and Leistra, M., 1980. Measured and simulated behaviour of aldicarb
and its oxidation products in fallow soils. Pesticide Science, 11 (4), 389-395.
Carrayrou, J., Mose, R. and Behra, P., 2004. Operator-splitting procedures for reactive
transport and comparison of mass balance errors. Journal of ContaminantHydrology, 68 (3/4), 239-268.
Casey, F.X.M., Hakk, H., ŠimĤnek, J., et al., 2004. Fate and transport of testosterone in
agricultural soils. Environmental Science and Technology, 38 (3), 790-798.
Casey, F.X.M. and ŠimĤnek, J., 2001. Inverse analyses of transport of chlorinatedhydrocarbons subject to sequential transformation reactions. Journal ofEnvironmental Quality, 30 (4), 1354-1360.
Casey, F.X.M., Larsen, G.L., Hakk, H., et al., 2003. Fate and transport of 17[beta]estradiol in soil-water systems. Environmental Science and Technology, 37 (11),2400-2409.
Castro, C.L. and Rolston, D.E., 1977. Organic phosphate transport and hydrolysis in soil:
theoretical and experimental evaluation. Soil Science Society of America Journal,
41 (6), 1085-1092.
Corapcioglu, M.Y. and Choi, H., 1996. Modeling colloid transport in unsaturated porous
media and validation with laboratory column data. Water Resources Research, 32(12), 3437-3449.
Cressey, G.B., 1958. Qanats, karez, and foggaras. The Geographical Review, 48, 27-44.
Crist, J.T., McCarthy, J.F., Zevi, Y., et al., 2004. Pore-scale visualization of colloidtransport and retention in partly saturated porous media. Vadose Zone Journal, 3(2), 444-450. [http://vzj.scijournals.org/cgi/content/full/3/2/444]
Darcy, H., 1856. Les fontaines publiques de la ville de Dyon. Dalmont, Paris.
Dautov, R.Z., Egorov, A.G., Nieber, J.L., et al., 2002. Simulation of two-dimensional
gravity-driven unstable flow. In: Hassanizadeh, S.M., Schotting, R.J., Gray,
W.G., et al. eds. Computational methods in water resources: proceedings of the
XIVth international conference on computational methods in water resources(CMWR XIV), June 23-28, 2002, Delft, The Netherlands. Volume 1. Elsevier,
Amsterdam, 9-16. Developments in Water Science no. 47.
De Jonge, L.W., Kjaergaard, C. and Moldrup, P., 2004. Colloids and colloid-facilitatedtransport of contaminants in soils: an introduction. Vadose Zone Journal, 3 (2),321-325. [http://vzj.scijournals.org/cgi/content/full/3/2/321]
De Rooij, G.H., 2000. Modeling fingered flow of water in soils owing to wetting front
instability: a review. Journal of Hydrology, 231/232, 277-294.

61

Chapter 2
De Vos, J.A., ŠimĤnek, J., Raats, P.A.C., et al., 1999. Identification of the hydrauliccharacteristics of a layered silt loam soil. In: Van Genuchten, M.T., Leij, F.J. andWu, L. eds. Characterisation and measurement of the hydraulic properties of
unsaturated porous media: proceedings of the international workshop. Part 1.
University of California, Riverside, 783-798.
DeNovio, N.M., Saiers, J.E. and Ryan, J.N., 2004. Colloid movement in unsaturated
porous media: recent advances and future directions. Vadose Zone Journal, 3 (2),338-351. [http://vzj.scijournals.org/cgi/content/full/3/2/338]
Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore
structure. Water Resources Research, 30 (2), 211-223.
Dutt, G.R., Shaffer, M.J. and Moore, W.J., 1972. Computer simulation model of dynamic
bio-physicochemical processes in soils. University of Arizona, Tucson. Technicalbulletin / Department of Soils, Water and Engineering, University of Arizona no.
196.
Eliassi, M. and Glass, R.J., 2002. On the porous-continuum modeling of gravity-driven
fingers in unsaturated materials: extension of standard theory with a hold-backpile-up effect. Water Resources Research, 38 (11), 161-1611.
Elimelech, M. and O' Melia, C.R., 1990. Kinetics of deposition of colloidal particles inporous media. Environmental Science and Technology, 24 (10), 1528-1536.
Essaid, H.I., Bekins, B.A., Godsy, E.M., et al., 1995. Simulation of aerobic andanaerobic biodegradation processes at a crude oil spill site. Water Resources
Research, 31 (12), 3309-3327.
Evans, D.D., Nicholson, T.J. and Rasmusson, T.C. (eds.), 2001. Flow and transport
through unsaturated fractured rock. 2nd edn. American Geophysical Union,
Washington. Geophysical Monograph no. 42.
Fryar, A.E. and Schwartz, F.W., 1994. Modeling the removal of metals fromgroundwater by a reactive barrier: experimental results. Water Resources
Research, 30 (12), 3455-3469.
Gerke, H.H. and Van Genuchten, M.T., 1993. A dual-porosity model for simulating thepreferential movement of water and solutes in structured porous media. Water
Resources Research, 29 (2), 305-319.
Gerke, H.H. and Van Genuchten, M.T., 1996. Macroscopic representation of structuralgeometry for simulating water and solute movement in dual-porosity media.
Advances in Water Resources, 19 (6), 343-357.
Ginn, T.R., Wood, B.D., Nelson, K.E., et al., 2002. Processes in microbial transport in
the natural subsurface. Advances in Water Resources, 25 (8/12), 1017-1042.
Govindaraju, R.S. and Kavvas, M.L., 1991. Dynamics of moving boundary overland
flows over infiltrating surfaces at hillslopes. Water Resources Research, 27 (8),1885-1898.
Grolimund, D., Borkovec, M., Barmettler, K., et al., 1996. Colloid-facilitated transport ofstrongly sorbing contaminants in natural porous media: a laboratory column
study. Environmental Science and Technology, 30 (10), 3118-3123.
Hanselman, T.A., Graetz, D.A. and Wilkie, A.C., 2003. Manure-borne estrogens as
potential environmental contaminants: a review. Environmental Science andTechnology, 37 (24), 5471-5478.
Harvey, R.W. and Garabedian, S.P., 1991. Use of colloid filtration theory in modelingmovement of bacteria through a contaminated sandy aquifer. Environmental
Science and Technology, 25 (1), 178-185.
Hassanizadeh, S.M., Celia, M.A. and Dahle, H.K., 2002. Dynamic effect in the capillary

62

Van Genuchten and ŠimĤnek
pressure-saturation relationship and its impacts on unsaturated flow. Vadose Zone
Journal, 1 (1), 38-57. [http://vzj.scijournals.org/cgi/content/full/1/1/38]
Hornberger, G.M., Mills, A.L. and Herman, J.S., 1992. Bacterial transport in porous
media: evaluation of a model using laboratory observations. Water Resources
Research, 28, 915-938.
Hutson, J.L. and Wagenet, R.J., 1995. A multiregion model describing water flow andsolute transport in heterogeneous soils. Soil Science Society of America Journal,
59 (3), 743-751.
Jacques, D., ŠimĤnek, J., Mallants, D., et al., 2003. The HYDRUS-PHREEQC
multicomponent transport model for variably-saturated porous media: code
verification and application. In: Poeter, E., Zheng, C., Hill, M., et al. eds.
MODFLOW and more 2003: understanding through modeling: conference
proceedings, International Ground Water Modeling Center, Colorado School ofMines, September 16-19, 2003. Colorado School of Mines, Golden, 23-27.
Jarvis, N., 1998. Modeling the impact of preferential flow on nonpoint source pollution.
In: Selim, H.M. and Ma, L. eds. Physical nonequilibrium in soils: modeling andapplication. Ann Arbor Press, Chelsea, 195-221.
Johnson, P.R. and Elimelech, M., 1995. Dynamics of colloid deposition in porous media:
blocking based on random sequential adsorption. Langmuir, 11, 801-812.
Jury, W.A., Spencer, W.F. and Farmer, W.J., 1983. Behavior assessment model for traceorganics in soil. I. Model description. Journal of Environmental Quality, 12 (4),558-564.
Kan, A.T. and Tomson, M.B., 1990. Ground water transport of hydrophobic organic
compounds in the presence of dissolved organic matter. Environmental
Toxicology and Chemistry, 9 (3), 253-263.
Kretzschmar, R., Borkovec, M., Grolimund, D., et al., 1999. Mobile subsurface colloids
and their role in contaminant transport. Advances in Agronomy, 66, 121-193.
Lawton, H.W. and Wilke, P.J., 1979. Ancient agricultural systems in dry regions. In:
Hall, A.E., Cannell, G.H. and Lawton, H.W. eds. Agriculture in semi-arid
environments. Springer, Berlin, 1-39.
Leij, F.J., Alves, W.J., Van Genuchten, M.T., et al., 1996. The UNSODA Unsaturated
Soil Hydraulic Database: user's manual version 1.0. National Risk Management
Laboratory, Cincinnati. Technical Report / Environmental Protection Agency no.
EPA/600/R-96/095.
Leij, F.J., ŠimĤnek, J., Suarez, D.L., et al., 1999. Nonequilibrium and multicomponent
transport models. In: Skaggs, R.W. and Van Schilfgaarde, J. eds. Agriculturaldrainage. American Society of Agronomy ASA, Madison, 405-430. Agronomy
Monograph no. 38.
Lester, D.H., Jansen, G. and Burkholder, H.C., 1975. Migration of radionuclide chainsthrough an adsorbing medium. In: Abrams, I.M., Zwiebel, I. and Sweed, N.H.
eds. Adsorption and ion exchange. American Institute of Chemical Engineers,
202-213. AIChE Symposium Series no. 152, vol. 71.
Lichtner, P.C., 1996. Continuum formulation of multicomponent-multiphase reactive
transport. In: Lichtner, P.C., Steefel, C.I. and Oelkers, E.H. eds. Reactive
transport in porous media. Mineralogical Society of America, Washington, 1-81.
Reviews in Mineralogy no. 34.
Lichtner, P.C., 2000. FLOTRAN user manual. Los Alamos National Laboratory, Los
Alamos. Report LA-UR-01-2349.
Lichtner, P.C. and Seth, M.S., 1996. User's manual for MULTIFLO. Part II.

63

Chapter 2
MULTIFLOW 1.0 and GEM 1.0: multicomponent-multiphase reactive transportmodel. Center for Nuclear Waste Regulatory Analyses, San Antonio. Report
CNWRA no. 96-010.
Lindqvist, R. and Enfield, C.G., 1992. Biosorption of dichlorodiphenyltrichloroethaneand hexachlorobenzene in groundwater and its implications for facilitatedtransport. Applied and Environmental Microbiology, 58 (7), 2211-2218.
Liu, C.W. and Narasimhan, T.N., 1989. Redox-controlled multiple-species reactive
chemical transport. 1. Model development. Water Resources Research, 25 (5),869-882.
Liu, H.H., Doughty, C. and Bodvarsson, G.S., 1998. An active fracture model forunsaturated flow and transport in fractured rocks. Water Resources Research, 34(10), 2633-2646.
Liu, H.H., Zhang, G. and Bodvarsson, G.S., 2003. The active fracture model: its relation
to fractal flow patterns and an evaluation using field observations. Vadose Zone
Journal, 2 (2), 259-269. [http://vzj.scijournals.org/cgi/content/full/2/2/259]
Logan, B.E., Jewett, D.G., Arnold, R.G., et al., 1995. Clarification of clean-bed filtration
models. Journal of Environmental Engineering, 121 (12), 869-873.
Magee, B.R., Lion, L.W. and Lemley, A.T., 1991. Transport of dissolved organic
macromolecules and their effect on the transport of phenanthrene in porous
media. Environmental Science and Technology, 25 (2), 323-331.
Mansfeldt, T., Leyer, H., Barmettler, K., et al., 2004. Cyanide leaching from soil
developed from coking plant purifier waste as influenced by citrate. Vadose Zone
Journal, 3 (2), 471-479. [http://vzj.scijournals.org/cgi/content/full/3/2/471]
Mattson, E.D., Ankeny, M.D., ŠimĤnek, J., et al., 2004. Overland flow implications on
surface cover designs, designing, building, & regulating evapotranspiration (ET)landfill covers, March 9-10, 2004, Denver, Colorado. U.S. Environmental
Protection Agency's Remediation Technology Development Forum (RTDF)
Program, Grand Junction.
McCarthy, J.F. and McKay, L.D., 2004. Colloid transport in the subsurface: past, present
and future challenges. Vadose Zone Journal, 3 (2), 326-337.
[http://vzj.scijournals.org/cgi/content/full/3/2/326]
Misra, C., Nielsen, D.R. and Biggar, J.W., 1974. Soil nitrogen transformations duringleaching. I. Theoretical considerations. Soil Science Society of America
Proceedings, 38 (2), 289-293.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., et al., 1998. Preferential transport ofnitrate to a tile drain in an intermittent-flood-irrigated field: Model development
and experimental evaluation. Water Resources Research, 34 (5), 1061-1076.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., et al., 1997. New piecewisecontinuous hydraulic functions for modeling preferential flow in an intermittentflood-irrigated field. Water Resources Research, 33 (9), 2049-2063.
Morris, E.M. and Woolhiser, D.A., 1980. Unsteady one-dimensional flow over a plane:
partial equilibrium and recession hydrographs. Water Resources Research, 16 (2),355-360.
National Research Council, 2001. Conceptual models of flow and transport in thefractured vadose zone. National Academy Press, Washington.
Noell, A.L., Thompson, J.L., Corapcioglu, M.Y., et al., 1998. The role of silica colloids
on facilitated cesium transport through glass bead columns and modeling.
Journal of Contaminant Hydrology, 31 (1/2), 23-56.
Panday, S. and Huyakorn, P.S., 2004. A fully coupled physically-based spatially-64

Van Genuchten and Šimůnek
distributed model for evaluating surface/subsurface flow. Advances in Water
Resources, 27 (4), 361-382.
Parkhurst, D.L. and Appelo, C.A.J., 1999. User's guide to PHREEQC (version 2): a
computer program for speciation, batch-reaction, one-dimensional transport, andinverse geochemical calculations. U.S. Department of the Interior, U.S.
Geological Survey, Denver. Water-Resources Investigations Report no. 99-4259.
Peters, R.R. and Klavetter, E.A., 1988. A continuum model for water movement in an
unsaturated fractured rock mass. Water Resources Research, 24, 416-430.
Pinder, G.F. and Abriola, L.M., 1986. On the simulation of nonaqueous phase organic
compounds in the subsurface. Water Resources Research, 22 (9 suppl.), 109-119.
Pruess, K., 1991. TOUGH2: a general-purpose numerical simulator for multiphase fluidand heat flow. Lawrence Berkeley Laboratory, Berkeley. Report no. LBL-29400.
Pruess, K. and Battistelli, A., 2002. TMVOC: a numerical simulator for three-phase nonisothermal flows of multicomponent hydrocarbon mixtures in saturatedunsaturated heterogeneous media. Lawrence Berkeley National Laboratory,
Berkeley. Report LBNL no. 49375.
Pruess, K., Oldenburg, C. and Moridis, G., 1999. TOUGH2 user's guide: version 2.0.
Lawrence Berkeley National Laboratory, Berkeley. Report LBNL no. 43134.
Rajagopalan, R. and Tien, C., 1976. Trajectory analysis of deep-bed filtration with thesphere-in-cell porous media model. AIChE Journal, 22, 523-533.
Rhoades, J.D., 1997. Sustainability of irrigation: an overview of salinity problems andcontrol strategies. In: Canadian Water Resources Association 50th annual
conference: footprints of humanity: reflections on fifty years of water resource
developments, Lethbridge, Alberta, Canada, June 3-6, 1997. CWRA, 1-42.
Richards, L.A., 1931. Capillary conduction of liquids in porous media. Physics, 1, 318333.
Ritsema, C.J., Dekker, L.W., Hendrickx, J.M.H., et al., 1993. Preferential flowmechanism in a water repellent sandy soil. Water Resources Research, 29 (7),2183-2193.
Rittmann, B.E. and VanBriesen, J.M., 1996. Microbiological processes in reactivemodeling. In: Lichtner, P.C., Steefel, C.I. and Oelkers, E.H. eds. Reactive
transport in porous media. Mineralogical Society of America, Washington, 311334. Reviews in Mineralogy no. 34.
Robbins, C.W., Wagenet, R.J. and Jurinak, J.J., 1980. A combined salt transportchemical equilibrium model for calcareous and gypsiferous soils. Soil Science
Society of America Journal, 44 (6), 1191-1194.
Rockhold, M.L., Yarwood, R.R. and Selker, J.S., 2004. Coupled microbial and transport
processes in soils. Vadose Zone Journal, 3 (2), 368-383.
[http://vzj.scijournals.org/cgi/content/full/3/2/368]
Rogers, V.C., 1978. Migration of radionuclide chains in groundwater. Nuclear
Technology, 40 (3), 315-320.
Ross, P.J. and Smettem, K.R.J., 2000. A simple treatment of physical nonequilibrium
water flow in soils. Soil Science Society of America Journal, 64 (6), 1926-1930.
Rubin, J., 1983. Transport of reacting solutes in porous media: relation between
mathematical nature of problem formulation and chemical nature of reactions.
Water Resources Research, 19 (5), 1231-1252.
Scanlon, B.R., Keese, K., Reedy, R.C., et al., 2003. Variations in flow and transport inthick desert vadose zones in response to paleoclimatic forcing (0-90 kyr): field

65

Chapter 2
measurements, modeling, and uncertainties. Water Resources Research, 39 (7),1179.
Schaap, M.G. and Leij, F.J., 2000. Improved prediction of unsaturated hydraulicconductivity with the Mualem-van Genuchten model. Soil Science Society ofAmerica Journal, 64 (3), 843-851.
Schaerlaekens, J., Mallants, D., Šimůnek, J., et al., 1999. Numerical simulation oftransport and sequential biodegradation of chlorinated aliphatic hydrocarbons
using CHAIN_2D. Hydrological Processes, 13 (17), 2847-2859.
Selim, H.M., Xue, S.K. and Iskandar, I.K., 1995. Transport of 2,4,6-trinitrotoluene andhexahydro-1,3,5-trintro-1,3,5-triazine in soils. Soil Science, 160 (5), 328-339.
Seuntjens, P., Mallants, D., Šimůnek, J., et al., 2002. Sensitivity analysis of physical andchemical properties affecting field-scale cadmium transport in a heterogeneous
soil profile. Journal of Hydrology, 264 (1/4), 185-200.
Sheremata, T.W., Thiboutot, S., Ampleman, G., et al., 1999. Fate of 2,4,6-trinitrotoluene
and its metabolites in natural and model soil systems. Environmental Science andTechnology, 33 (22), 4002-4008.
Shikaze, S.G., Sudicky, E.A. and Mendoza, C.A., 1994. Simulation of dense vapor
migration in discretely fractured geologic media. Water Resources Research, 30(7), 1993-2009.
Šimůnek, J., Šejna, M. and Van Genuchten, M.T, 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. International
Ground Water Modeling Center, Colorado School of Mines, Golden.
Šimůnek, J., Šejna, M. and Van Genuchten, M.T., 1998. The HYDRUS-1D software
package for simulating the one-dimensional movement of water, heat, andmultiple solutes in variably-saturated media: version 2.0. IGWMC-TPS-70.
International Ground Water Modeling Center, Colorado School of Mines,
Golden.
Šimůnek, J., Huang, K. and Van Genuchten, M.T., 1995. The SWMS_3D code forsimulating water flow and solute transport in three-dimensional variably
saturated media: version 1.0. USDA-ARS, U.S. Salinity Laboratory, Riverside,
California. Research Report no. 139.
Šimůnek, J., Jarvis, N.J., Van Genuchten, M.T., et al., 2003. Review and comparison ofmodels for describing non-equilibrium and preferential flow and transport in thevadose zone. Journal of Hydrology, 272 (1/4), 14-35.
Šimůnek, J. and Suarez, D.L., 1994. Two-dimensional transport model for variably
saturated porous media with major ion chemistry. Water Resources Research, 30(4), 1115-1133.
Šimůnek, J. and Suarez, D.L., 1997. Sodic soil reclamation using multicomponent
transport modeling. Journal of Irrigation and Drainage Engineering ASCE, 123(5), 367-376.
Šimůnek, J. and Valocchi, A.J., 2002. Geochemical transport. In: Dane, J.H. and Topp,
G.C. eds. Methods of soil analysis. Part 4. Physical methods. SSSA, Madison,1511-1536. SSSA Book Series no. 5.
Steefel, C.I., 2000. New directions in hydrogeochemical transport modeling:
incorporating multiple kinetic and equilibrium reaction pathways. In: Bentley,
L.R., Sykes, J.F., Brebbia, C.A., et al. eds. Computational methods in water
resources : proceedings 13th international conference. Balkema, Rotterdam,331-338.

66

Van Genuchten and Šimůnek
Steefel, C.I. and MacQuarrie, K.T.B., 1996. Approaches to modeling of reactive
transport in porous media. In: Lichtner, P.C., Steefel, C.I. and Oelkers, E.H. eds.
Reactive transport in porous media. Mineralogical Society of America,
Washington, 82-129. Reviews in Mineralogy no. 34.
Steefel, C.I. and Yabusaki, S.B., 1996. OS3D/GIMRT: software for multicomponentmultidimensional reactive transport: user's manual and programmer's guide.
Pacific Northwest National Laboratory, Richland. Report PNNL no. 11166.
Stephens, D.B., Kowall, S.J., Borns, D., et al., 2002. Letter to the editor on A national
strategy for vadose zone science and technology. Vadose Zone Journal, 1 (1),
197-198. [http://vzj.scijournals.org/cgi/content/full/1/1/197]
Strelkoff, T.S., Clemmens, A.J. and Schmidt, B.V., 1998. SRFR version 3.31: a model
for simulating surface irrigation in borders, basins and furrows. U.S. Water
Conservation Laboratory, USDA-ARS, Phoenix.
Tanji, K.K., 1990. Agricultural salinity assessment and management. ASCE, New York.
ASCE Manuals and Reports on Engineering Practice no. 71.
Therrien, R. and Sudicky, E.A., 1996. Three-dimensional analysis of variably-saturated
flow and solute transport in discretely-fractured porous media. Journal ofContaminant Hydrology, 23 (1/2), 1-44.
Thiele-Bruhn, S., 2003. Pharmaceutical antibiotic compounds in soils: a review. Journalof Plant Nutrition and Soil Science, 166 (2), 145-167.
Thompson, S.S. and Yates, M.V., 1999. Bacteriophage inactivation at the air-water-solidinterface in dynamic batch systems. Applied and Environmental Microbiology,
65, 1186-1190.
Tolls, J., 2001. Sorption of veterinary pharmaceuticals in soils: a review. Environmental
Science and Technology, 35 (17), 3397-3406.
Totsche, K.U. and Kögel-Knabner, I., 2004. Mobile organic sorbent affectedcontaminant transport in soil: numerical case studies for enhanced and reduced
mobility. Vadose Zone Journal, 3 (2), 352-367.
[http://vzj.scijournals.org/cgi/content/full/3/2/352]
U.S. Department of Energy, 2001. A national roadmap for vadose zone science andtechnology: understanding, monitoring, and predicting contaminant fate andtransport in the unsaturated zone. Revision 0.0. U.S. Department of Energy.
Report DOE no. ID-10871.
[http://www.inel.gov/vadosezone/documents/roadmap.pdf]
Valocchi, A.J., Street, R.L. and Roberts, P.V., 1981. Transport of ion-exchanging solutes
in groundwater: chromatographic theory and field simulation. Water Resources
Research, 17 (5), 1517-1527.
Van Dam, J.C., Huygen, J., Wesseling, J.G., et al., 1997. Theory of SWAP version 2.0:
simulation of water flow, solute transport and plant growth in the Soil-WaterAtmosphere-Plant environment. DLO Winand Staring Centre, Wageningen.
Report / Wageningen Agricultural University, Department Water Resources no.
71.
Van Genuchten, M.T., 1985. Convective-dispersive transport of solutes involved insequential first-order decay reactions. Computers and Geosciences, 11 (2), 129147.
Van Genuchten, M.T. and Wagenet, R.J., 1989. Two-site/two-region models forpesticide transport and degradation: theoretical development and analyticalsolutions. Soil Science Society of America Journal, 53 (5), 1303-1310.
Van Genuchten, M.T. and Wierenga, P.J., 1976. Mass transfer studies in sorbing porous

67

Chapter 2
media. I. Analytical Solutions. Soil Science Society of America Journal, 40 (4),473-480.
VanBriesen, J.M., 1998. Modeling coupled biogeochemical processes in mixed waste
systems. PhD Thesis Department of Civil Engineering, Northwestern University,
Evanston.
VanderKwaak, J.E., 1999. Numerical simulation of flow and chemical transport in
integrated surface-subsurface hydrologic systems. PhD Thesis Department ofEarth Sciences, University of Waterloo, Ontario.
VanderKwaak, J.E. and Loague, K., 2001. Hydrologic-response simulations for the R-5
catchment with a comprehensive physics-based model. Water Resources
Research, 37 (4), 999-1013.
Vinten, A.J.A., Yaron, B. and Nye, P.H., 1983. Vertical transport of pesticides into soil
when adsorbed on suspended particles. Journal of Agricultural and Food
Chemistry, 31 (3), 662-664.
Viswanathan, H.S., Robinson, B.A., Valocchi, A.J., et al., 1998. A reactive transportmodel of neptunium migration from the potential repository at Yucca Mountain.
Journal of Hydrology, 209 (1/4), 251-280.
Von Gunten, H.R., Waber, U.E. and Krähenbühl, U., 1988. The reactor accident at
Chernobyl: a possibility to test colloid-controlled transport of radionuclides in ashallow aquifer. Journal of Contaminant Hydrology, 2 (3), 237-247.
Wagenet, R.J. and Hutson, J.L., 1987. LEACHM: leaching, estimation and chemistrymodel: a process-based model of water and solute movement, transformations,plant uptake and chemical reactions in the unsaturated zone. Water Resources
Institute, Cornell University, Ithaca.
Walter, A.L., Frind, E.O., Blowes, D.W., et al., 1994. Modeling of multicomponent
reactive transport in groundwater. 2. Metal mobility in aquifers impacted byacidic mine tailings discharge. Water Resources Research, 30 (11), 3149-3158.
Wan, J. and Tokunaga, T.K., 2002. Partitioning of clay colloids at air-water interfaces.
Journal of Colloid and Interface Science, 247 (1), 54-61.
Wang, Z., Jury, W.A., Tuli, A., et al., 2004. Unstable flow during redistribution:
controlling factors and practical implications. Vadose Zone Journal, 3 (2), 549559. [http://vzj.scijournals.org/cgi/content/full/3/2/549]
Wilson, G.V., Jardine, P.M. and Gwo, J.P., 1992. Modeling the hydraulic properties of amultiregion soil. Soil Science Society of America Journal, 56 (6), 1731-1737.
Woolhiser, D.A. and Liggett, J.A., 1967. Unsteady, one-dimensional flow over a plane:
the rising hydrograph. Water Resources Research, 3 (3), 753-771.
Woolhiser, D.A., Smith, R.E. and Goodrich, D.C., 1990. KINEROS: a kinematic runoff
and erosion model: documentation and user manual. U.S. Department ofAgriculture, Agricultural Research Service. Report no. ARS-77.
Yeh, G.T., Salvage, K.M., Gwo, J.P., et al., 1998. HYDROBIOGEOCHEM: a coupledmodel of hydrologic transport and mixed biochemical kinetic/equilibrium
reactions in saturated-unsaturated media. Oak Ridge National Laboratory, Oak
Ridge. Report no. ORNL/TM-13668.
Yeh, G.T. and Tripathi, V.S., 1990. HYDROGEOCHEM: a coupled model ofhydrological and geochemical equilibrium of multicomponent systems. Oak
Ridge National Laboratory, Oak Ridge. Report no. ORNL-6371.
Zhang, P.F., Šimůnek, J. and Bowman, R.S., 2004. Nonideal transport of solute andcolloidal tracers through reactive zeolite/iron pellets. Water Resources Research,
40 (4), NIL_1-NIL_11.

68

Van Genuchten and Šimůnek
Zhang, W. and Cundy, T.W., 1989. Modeling of two-dimensional overland flow. Water
Resources Research, 25 (9), 2019-2035.
Zurmühl, T. and Durner, W., 1996. Modeling transient water and solute transport in abiporous soil. Water Resources Research, 32 (4), 819-829.
Zyvoloski, G.A., Robinson, B.A., Dash, Z.V., et al., 1997. Summary of the models andmethods for the FEHM application: a finite element heat- and mass-transfer
code. Los Alamos National Laboratory, Los Alamos. Report no. LA-13307-MS.

69