JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 1 — #1

22

Contaminant

Transport in the

Unsaturated Zone

Theory and Modeling

Ji ̆rí Šim  ̊unek

University of California Riverside

Martinus Th. van

Genuchten

George E. Brown, Jr, Salinity

Laboratory, USDA-ARS

22.1 Introduction.............................................. 22-1

22.2 Variably Saturated Water Flow .......................... 22-3

Mass Balance Equation • Uniform Flow • Preferential Flow

22.3 Solute Transport ......................................... 22-8

Transport Processes • Advection–Dispersion Equations •

Nonequilibrium Transport • Stochastic Models •

Multicomponent Reactive Solute Transport • Multiphase

Flow and Transport • Initial and Boundary Conditions

22.4 Analytical Models ........................................ 22-22

Analytical Approaches • Existing Models

22.5 Numerical Models ....................................... 22-27

Numerical Approaches • Existing Models

22.6 Concluding Remarks .................................... 22-36

References ....................................................... 22-38

22.1 Introduction

Human society during the past several centuries has created a large number of chemical substances that

often find their way into the environment, either intentionally applied during agricultural practices,

unintentionally released from leaking industrial and municipal waste disposal sites, or stemming from

research or weapons production related activities. As many of these chemicals represent a significant

health risk when they enter the food chain, contamination of both surface and subsurface water supplies

has become a major issue. Modern agriculture uses an unprecedented number of chemicals, both in

plant and animal production. A broad range of fertilizers, pesticides and fumigants are now routinely

applied to agricultural lands, thus making agriculture one of the most important sources for non-point

source pollution. The same is true for salts and toxic trace elements, which are often an unintended

consequence of irrigation in arid and semiarid regions. While many agricultural chemicals are generally

beneficial in surface soils, their leaching into the deeper vadose zone and groundwater may pose serious

problems. Thus, management processes are being sought to keep fertilizers and pesticides in the root

22-1
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 2 — #2

22-2 The Handbook of Groundwater Engineering

zone and prevent their transport into underlying or down-gradient water resources. Agriculture also

increasingly uses a variety of pharmaceuticals and hormones in animal production many of which, along

with pathogenic microorganisms, are being released to the environment through animal waste. Animal

waste and wash water effluent, in turn, is frequently applied to agricultural lands. Potential concerns about

the presence of pharmaceuticals and hormones in the environment include: (1) abnormal physiological

processes and reproductive impairment; (2) increased incidences of cancer; (3) development of antibiotic

resistant bacteria; and (4) increased toxicity of chemical mixtures. While the emphasis above is mostly on

non-point source pollution by agricultural chemicals, similar problems arise with point-source pollution

from industrial and municipal waste disposal sites, leaking underground storage tanks, chemicals spills,

nuclear waste repositories, and mine tailings, among other sources.

Mathematical models should be critical components of any effort to optimally understand and quantify

site-specific subsurface water flow and solute transport processes. For example, models can be helpful tools

for designing, testing and implementing soil, water, and crop management practices that minimize soil

and water pollution. Models are equally needed for designing or remediating industrial waste disposal sites

and landfills, or for long-term stewardship of nuclear waste repositories. A large number of specialized

numerical models now exist to simulate the different processes at various levels of approximation and for

different applications.

Increasing attention is being paid recently to the unsaturated or vadose zone where much of the

subsurface contamination originates, passes through, or can be eliminated before it contaminates surface

and subsurface water resources. Sources of contamination often can be more easily remediated in the

vadose zone, before contaminants reach the underlying groundwater. Other chapters in this Handbook

deal with water flow (Chapters xx) and solute transport (Chapters xx) in fully saturated (groundwater)

systems and with water flow in the unsaturated zone (Chapter 5). The focus of this chapter thus will be

on mathematical descriptions of transport processes in predominantly variably saturated media.

Soils are generally defined as the biologically active layer at the surface of the earth’s crust that is made

up of a heterogeneous mixture of solid, liquid, and gaseous material, as well as containing a diverse

community of living organisms (Jury and Horton, 2004). The vadose (unsaturated) zone is defined as

the layer between the land surface and the permanent (seasonal) groundwater table. While pores between

solid grains are fully filled with water in the saturated zone (groundwater), pores in the unsaturated zone

are only partially filled with water, with the remaining part of the pore space occupied by the gaseous

phase. The vadose zone is usually only partially saturated, although saturated regions may exist, such as

when perched water is present above a low-permeable fine-textured (clay) layer or a saturated zone behind

the infiltration front during or after a high-intensity rainfall event.

As the transport of contaminants is closely linked with the water flux in soils and rocks making up the

vadose zone, any quantitative analysis of contaminant transport must first evaluate water fluxes into and

through the vadose zone. Water typically enters the vadose zone in the form of precipitation or irrigation

(Figure 22.1), or by means of industrial and municipal spills. Some of the rainfall or irrigation water

may be intercepted on the leaves of vegetation. If the rainfall or irrigation intensity is larger than the

infiltration capacity of the soil, water will be removed by surface runoff, or will accumulate at the soil

surface until it evaporates back to the atmosphere or infiltrates into the soil. Part of this water is returned

to the atmosphere by evaporation. Some of the water that infiltrates into the soil profile may be taken

up by plant roots and eventually returned to the atmosphere by plant transpiration. The processes of

evaporation and transpiration are often combined into the single process of evapotranspiration. Only

water that is not returned to the atmosphere by evapotranspiration may percolate to the deeper vadose

zone and eventually reach the groundwater table. If the water table is close enough to the soil surface, the

process of capillary rise may move water from the groundwater table through the capillary fringe toward

the root zone and the soil surface.

Because of the close linkage between water flow and solute transport, we will first briefly focus on the

physics and mathematical description of water flow in the vadose zone (Section 22.2). An overview is given

of the governing equations for water flow in both uniform (Section 22.2.2) and structured (Section 22.2.3)

media. This section is followed by a discussion of the governing solute transport equations (Section 22.3),
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 3 — #3

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-3

Transpiration

Rainfall

Interception

Evaporation

Irrigation

Runoff

Root water

uptake

Deep drainage

Capillary rise

Groundwater recharge

Water table

Unsaturated zone

zone

Saturated

Root zone

Storage

FIGURE 22.1 Schematic of water fluxes and various hydrologic components in the vadose zone.

again for both uniform (Section 22.3.2) and structured (fractured) (Section 22.3.3) media. We also

briefly discuss alternative formulations for colloid (Section 22.3.3.2.2) and colloid-facilitated transport

(Section 22.3.3.3), multicomponent geochemical transport (Section 22.3.5), and stochastic approaches

for solute transport (Section 22.3.4). This is followed by a discussion of analytical (Section 22.4) and

numerical (Section 22.5) approaches for solving the governing flow and transport equations, and an

overview of computer models currently available for simulating vadose zone flow and transport processes

(Sections 22.4.2 and 22.5.2).

22.2 Variably Saturated Water Flow

In this section, we briefly present the equations governing variably saturated water flow in the subsurface. More details about this topic, including the description of the soil hydraulic properties and their

constitutive relationship, are given in Chapter xx. Traditionally, descriptions of variably saturated flow

in soils are based on the Richards (1931) equation, which combines the Darcy–Buckingham equation

for the fluid flux with a mass balance equation. The Richards equation typically predicts a uniform flow

process in the vadose zone, although possibly modified macroscopically by spatially variable soil hydraulic

properties (e.g., as dictated by the presence of different soil horizons, but possibly also varying laterally).
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 4 — #4

22-4 The Handbook of Groundwater Engineering

Unfortunately, the vadose zone can be extremely heterogeneous at a range of scales, from the microscopic

(e.g., pore scale) to the macroscopic (e.g., field or larger scale). Some of these heterogeneities can lead

to a preferential flow process that macroscopically is very difficult to capture with the standard Richards

equation. One obvious example of preferential flow is the rapid movement of water and dissolved solutes

through macropores (e.g., between soil aggregates, or created by earthworms or decayed root channels)

or rock fractures, with much of the water bypassing (short-circuiting) the soil or rock matrix. However,

many other causes of preferential flow exist, such as flow instabilities caused by soil textural changes or

water repellency (Hendrickx and Flury, 2001; Šim  ̊unek et al. 2003; Ritsema and Dekker 2005), and lateral

funneling of water due to inclined or other textural boundaries (e.g., Kung 1990). Alternative ways of

modeling preferential flow are discussed in a later section. Here we first focus on the traditional approach

for uniform flow as described with the Richards equation.

22.2.1 Mass Balance Equation

Water flow in variably saturated rigid porous media (soils) is usually formulated in terms of a mass balance

equation of the form:

∂θ

∂t = −∂qi∂xi

− S (22.1)

where θ is the volumetric water content [L3L−3], t is time [T], xi is the spatial coordinate [L], qi is the

volumetric flux density [LT−1], and S is a general sink/source term [L3L−3T−1], for example, to account

for root water uptake (transpiration). Equation 22.1 is often referred to as the mass conservation equation

or the continuity equation. The mass balance equation in general states that the change in the water

content (storage) in a given volume is due to spatial changes in the water flux (i.e., fluxes in and out of

some unit volume of soil ) and possible sinks or sources within that volume. The mass balance equation

must be combined with one or several equations describing the volumetric flux density (q) to produce the

governing equation for variably saturated flow. The formulations of the governing equations for different

types of flow (uniform and preferential flow) are all based on this continuity equation.

22.2.2 Uniform Flow

Uniform flow in soils is described using the Darcy–Buckingham equation:

qi = −K(h)



KA

ij

∂h

∂xj

+ KA

iz 

(22.2)

where K is the unsaturated hydraulic conductivity [LT−1], and KA

ij are components of a dimensionless

anisotropy tensor KA (which reduces to the unit matrix when the medium is isotropic). The Darcy–

Buckingham equation is formally similar to Darcy’s equation, except that the proportionality constant

(i.e., the unsaturated hydraulic conductivity) in the Darcy–Buckingham equation is a nonlinear function

of the pressure head (or water content), while K(h)in Darcy’s equation is a constant equal to the saturated

hydraulic conductivity, Ks (e.g., see discussion by Narasimhan [2005]).

Combining the mass balance Equation 22.1 with the Darcy–Buckingham Equation22.2 leads to the

general Richards equation (Richards, 1931)

∂θ(h)

∂t = ∂∂xi



K(h)



KA

ij

∂h

∂xj

+ KA

iz  − S(h) (22.3)

This partial differential equation is the equation governing variably saturated flow in the vadose zone.

Because of its strongly nonlinear makeup, only a relatively few simplified analytical solutions can be

derived. Most practical applications of Equation 22.3 require a numerical solution, which can be obtained

using a variety of numerical methods such as finite differences or finite elements (Section 22.5a).
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 5 — #5

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-5

Equation 22.3 is generally referred to as the mixed form of the Richards equation as it contains two

dependent variables, that is, the water content and the pressure head. Various other formulations of the

Richards equation are possible.

22.2.3 Preferential Flow

Increasing evidence exists that variably saturated flow in many field soils is not consistent with the

uniform flow pattern typically predicted with the Richards equations (Flury et al., 1994; Hendrickx and

Flury, 2001). This is due to the presence of macropores, fractures, or other structural voids or biological

channels through which water and solutes may move preferentially, while bypassing a large part of the

matrix pore-space. Preferential flow and transport processes are probably the most frustrating in terms

of hampering accurate predictions of contaminant transport in soils and fractured rocks. Contrary to

uniform flow, preferential flow results in irregular wetting of the soil profile as a direct consequence of

water moving faster in certain parts of the soil profile than in others. Hendrickx and Flury (2001) defined

preferential flow as constituting all phenomena where water and solutes move along certain pathways,

while bypassing a fraction of the porous matrix. Water and solutes for these reasons can propagate quickly

to far greater depths, and much faster, than would be predicted with the Richards equation describing

uniform flow.

The most important causes of preferential flow are the presence of macropores and other structural

features, development of flow instabilities (i.e., fingering) caused by profile heterogeneities or water

repellency (Hendrickx et al., 1993), and funneling of flow due to the presence of sloping soil layers that

redirect downward water flow. While the latter two processes (i.e., flow instability and funneling) are

usually caused by textural differences and other factors at scales significantly larger than the pore scale,

macropore flow and transport are usually generated at the pore or slightly larger scales, including scales

where soil structure first manifests itself (i.e., the pedon scale) (Šim  ̊unek et al., 2003).

Uniform flow in granular soils and preferential flow in structured media (both macroporous soils and

fractured rocks) can be described using a variety of single-porosity, dual-porosity, dual-permeability,

multi-porosity, and multi-permeability models (Richards, 1931; Pruess and Wang, 1987; Gerke and van

Genuchten, 1993a; Gwo et al., 1995; Jarvis, 1998; Šim  ̊unek et al., 2003, 2005). While single-porosity

models assume that a single pore system exists that is fully accessible to both water and solute, dualporosity and dual-permeability models both assume that the porous medium consists of two interacting

pore regions, one associated with the inter-aggregate, macropore, or fracture system, and one comprising

the micropores (or intra-aggregate pores) inside soil aggregates or the rock matrix. Whereas dual-porosity

models assume that water in the matrix is stagnant, dual-permeability models allow also for water flow

within the soil or rock matrix.

Figure 22.2 illustrates a hierarchy of conceptual formulations that can be used to model variably

saturated water flow and solute transport in soils. The simplest formulation (Figure 22.2a) is a singleporosity (equivalent porous medium) model applicable to uniform flow in soils. The other models

apply in some form or another to preferential flow or transport. Of these, the dual-porosity model

of Figure 22.2c assumes the presence of two pore regions, with water in one region being immobile

and in the other region mobile. This model allows the exchange of both water and solute between the

two regions (Šim  ̊unek et al., 2003). Conceptually, this formulation views the soil as consisting of a soil

matrix containing grains/aggregates with a certain internal microporosity (intra-aggregate porosity) and

a macropore or fracture domain containing the larger pores (inter-aggregate porosity). While water and

solutes are allowed to move through the larger pores and fractures, they can also flow in and out of

aggregates. By comparison, the intra-aggregate pores represent immobile pockets that can exchange,

retain, and store water and solutes, but do not contribute to advective (or convective) flow. Models that

assume mobile–immobile flow regions (Figure 22.2b) are conceptually somewhere in between the singleand dual-porosity models. While these models assume that water will move similarly as in the uniform

flow models, the liquid phase for purposes of modeling solute transport is divided in terms of mobile and

immobile fractions, with solutes allowed to move by advection and dispersion only in the mobile fraction
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 6 — #6

22-6 The Handbook of Groundwater Engineering

(a) Uniform flow

Water Water Immob

Immob. Immob.

Water Water

Solute Solute Solute

Solute

Mobile

Mobile Mobile

Slow

Slow

Fast

Fast

(b) Mobile-immobile water (c) Dual-porosity (d) Dual-permeability

 = im + mo   = im + mo  = m + f

FIGURE 22.2 Conceptual models of water flow and solute transport (θ is the water content, θmo and θim in (b) and

(c) are water contents in the mobile and immobile flow regions, respectively, and θm and θf in (d) are water contents

in the matrix and macropore (fracture) regions, respectively).

and between the two pore regions. This model has long been applied to solute transport studies (e.g., van

Genuchten and Wierenga, 1976).

Finally, dual-permeability models (Figure 22.2d) are those in which water can move in both the interand intra-aggregate pore regions (and matrix and fracture domains). These models in various forms are

now also becoming increasingly popular (Pruess and Wang, 1987; Gerke and van Genuchten, 1993a; Jarvis,

1994; Pruess, 2004). Available dual-permeability models differ mainly in how they implement water flow

in and between the two pore regions (Šim  ̊unek et al., 2003). Approaches to calculating water flow in the

macropores or inter-aggregate pores range from those invoking Poiseuille’s equation (Ahuja and Hebson,

1992), the Green and Ampt or Philip infiltration models (Ahuja and Hebson, 1992), the kinematic wave

equation (Germann, 1985; Germann and Beven, 1985; Jarvis, 1994), and the Richards equation (Gerke

and van Genuchten, 1993a). Multi-porosity and multi-permeability models (not shown in Figure 22.2)

are based on the same concept as dual-porosity and dual-permeability models, but include additional

interacting pore regions (e.g., Gwo et al., 1995; Hutson and Wagenet, 1995). These models can be readily

simplified to the dual-porosity/permeability approaches. Recent reviews of preferential flow processes and

available mathematical models are provided by Hendrickx and Flury (2001) and Šim  ̊unek et al. (2003),

respectively.

22.2.3.1 Dual-Porosity Models

Dual-porosity models assume that water flow is restricted to macropores (or inter-aggregate pores and

fractures), and that water in the matrix (intra-aggregate pores or the rock matrix) does not move at all.

This conceptualization leads to two-region type flow and transport models (van Genuchten and Wierenga,

1976) that partition the liquid phase into mobile (flowing, inter-aggregate), θmo, and immobile (stagnant,

intra-aggregate), θim, regions [L3L−3]:

θ = θmo + θim (22.4)

The dual-porosity formulation for water flow can be based on a mixed formulation of the Richards

Equation 22.3 to describe water flow in the macropores (the preferential flow pathways) and a mass

balance equation to describe moisture dynamics in the matrix as follows (Šim  ̊unek et al., 2003):

∂θmo(hmo)

∂t = ∂∂xi



K(hmo)



KA

ij

∂hmo

∂xj

+ KA

iz  − Smo(hmo) − w

∂θim(him)

∂t = −Sim(him) + w

(22.5)
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 7 — #7

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-7

where Sim and Smo are sink terms for both regions [T−1], and w is the transfer rate for water from the

inter- to the intra-aggregate pores [T−1].

Several of the above dual-porosity features were recently included in the HYDRUS software packages

(Šim  ̊unek et al., 2003, 2005). Examples of their application to a range of laboratory and field data involving

transient flow and solute transport are given by Šim  ̊unek et al. (2001), Zhang et al. (2004), Köhne et al.

(2004a, 2005), Kodešová et al. (2005), and Haws et al. (2005).

22.2.3.2 Dual-Permeability Models

Different types of dual-permeability approaches may be used to describe flow and transport in structured

media. While several models invoke similar governing equations for flow in the fracture and matrix

regions, others use different formulations for the two regions. A typical example of the first approach is

the work of Gerke and van Genuchten (1993a, 1996) who applied Richards equations to each of two pore

regions. The flow equations for the macropore (fracture) (subscript f) and matrix (subscript m) pore

systems in their approach are given by:

∂θf(hf)

∂t = ∂∂xi



Kf(hf)



KA

ij

∂hf

∂xj

+ KA

iz  − Sf(hf) − w

w (22.6)

and

∂θm(hm)

∂t = ∂∂xi



Km(hm)



KA

ij

∂hm

∂xj

+ KA

iz  − Sm(hm) +

w

1 − w (22.7)

respectively, where w is the ratio of the volumes of the macropore (or fracture or inter-aggregrate)

domain and the total soil system [−]. This approach is relatively complicated in that the model requires

characterization of water retention and hydraulic conductivity functions (potentially of different form)

for both pore regions, as well as the hydraulic conductivity function of the fracture–matrix interface.

Note that the water contents θf and θm in Equation 22.6 and Equation 22.7 have different meanings than

in Equation 22.5 where they represented water contents of the total pore space (i.e., θ = θmo + θim),

while here they refer to water contents of the two separate (fracture or matrix) pore domains such that

θ = wθf + (1 − w)θm.

22.2.3.3 Mass Transfer

The rate of exchange of water between the macropore and matrix regions, w, is a critical term in both the

dual-porosity model (Equation 22.5) and the dual-permeability approach given by (Equation 22.6) and

(Equation 22.7). Gerke and van Genuchten (1993a) assumed that the rate of exchange is proportional to

the difference in pressure heads between the two pore regions:

w = αw(hf − hm) (22.8)

in which αw is a first-order mass transfer coefficient [T−1]. For porous media with well-defined geometries,

the first-order mass transfer coefficient, αw, can be defined as follows (Gerke and van Genuchten, 1993b):

αw = β

d2 Kaγw (22.9)

where d is an effective diffusion path length [L] (i.e., half the aggregate width or half the fracture spacing),

β is a shape factor that depends on the geometry [−], and γw (= 0.4) is a scaling factor [−] obtained

by matching the results of the first-order approach at the half-time level of the cumulative infiltration

curve to the numerical solution of the horizontal infiltration equation (Gerke and van Genuchten, 1993b).

Several other approaches based on water content of relative saturation differences have also been used

(Šim  ̊unek et al., 2003).
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 8 — #8

22-8 The Handbook of Groundwater Engineering

22.3 Solute Transport

Similarly, as shown in Equation 22.1 for water flow, mathematical formulations for solute transport are

usually based on a mass balance equation of the form:

∂CT

∂t = −∂JTi∂xi

− φ (22.10)

where CT is the total concentration of chemical in all forms [ML−3], JTi is the total chemical mass flux

density (mass flux per unit area per unit time) [ML−2T−1], and φ is the rate of change of mass per unit

volume by reactions or other sources (negative) or sinks (positive) such as plant uptake [ML−3T−1]. In

its most general interpretation, Equation 22.10 allows the chemical to reside in all three phases of the

soil (i.e., gaseous, liquid, and solid), permits a broad range of transport mechanisms (including advective

transport, diffusion, and hydrodynamic dispersion in both the liquid and gaseous phases), and facilitates

any type of chemical reaction that leads to losses or gains in the total concentration.

While the majority of chemicals are present only in the liquid and solid phases, and as such are

transported in the vadose zone mostly only in and by water, some chemicals such as many organic

contaminants, ammonium, and all fumigants, can have a significant portion of their mass in the gaseous

phase and are hence subject to transport in the gaseous phase as well. The total chemical concentration

can thus be defined as:

CT = ρbs + θc + ag (22.11)

where ρb is the bulk density [ML−3], θ is the volumetric water content [L3L−3], a is the volumetric air

content [L3L−3], and s[MM−1], c[ML−3], and g [ML−3] are concentrations in the solid, liquid, and

gaseous phases, respectively. The solid phase concentration represents solutes sorbed onto sorption sites

of the solid phase, but can include solutes sorbed onto colloids attached to the solid phase or strained by

the porous system, and solutes precipitated onto or into the solid phase.

The reaction term φ of Equation 22.10 may represent various chemical or biological reactions that lead

to a loss or gain of chemical in the soil system, such as radionuclide decay, biological degradation, and

dissolution. In analytical and numerical models these reactions are most commonly expressed using zeroand first-order reaction rates as follows:

φ = ρbsμs + θcμw + agμg − ρbγs − θγw − aγg (22.12)

where μs, μw, and μg are first-order degradation constants in the solid, liquid, and gaseous phases [T−1],

respectively, and γs [T−1], γw [ML−3T−1], and γg [ML−3T−1] are zero-order production constants in the

solid, liquid, and gaseous phases, respectively.

22.3.1 Transport Processes

When a solute is present in both the liquid and gaseous phase, then various transport processes in both of

these phases may contribute to the total chemical mass flux:

JT = Jl + Jg (22.13)

where Jl and Jg represent solute fluxes in the liquid and gaseous phases [ML−2T−1], respectively. Note

that in Equation 22.13, and further below, we omitted the subscript i accounting for the direction of flow.

The three main processes that can be active in both the liquid and gaseous phase are molecular diffusion,

hydrodynamic dispersion, and advection (often also called convection). The solute fluxes in the two phases
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 9 — #9

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-9

are then the sum of fluxes due to these different processes:

Jl = Jlc + Jld + Jlh

Jg = Jgc + Jgd + Jgh

(22.14)

where the additional subscripts c, d, and h denote convection (or advection), molecular diffusion, and

hydrodynamic dispersion, respectively.

22.3.1.1 Diffusion

Diffusion is a result of the random motion of chemical molecules. This process causes a solute to move

from a location with a higher concentration to a location with a lower concentration. Diffusive transport

can be described using Fick’s law:

Jld = −θξl(θ )Dw

l

∂c

∂z = −θDs

l

∂c

∂z

Jgd = −aξg(θ )Dw

g

∂g

∂z = −aDs

g

∂g

∂z

(22.15)

where Dw

l and Dw

g are binary diffusion coefficients of the solute in water and gas [L2T−1], respectively; Ds

l

and Ds

g are the effective diffusion coefficients in soil water and soil gas [L2T−1], respectively; and ξl and

ξg are tortuosity factors that account for the increased path lengths and decreased cross-sectional areas

of the diffusing solute in both phases (Jury and Horton, 2004). As solute diffusion in soil water (air) is

severely hampered by both air (water) and solid particles, the tortuosity factor increases strongly with

water content (air content). Many empirical models have been suggested in the literature to account for

the tortuosity (e.g., Moldrup et al., 1998). Among these, the most widely used model for the tortuosity

factor is probably the equation of Millington and Quirk (1961) given by:

ξl(θ ) = θ 7/3

θ 2

s

(22.16)

where θs is the saturated water content (porosity) [L3L−3]. A similar equation may be used for the

tortuosity factor of the gaseous phase by replacing the water content with the air content.

22.3.1.2 Dispersion

Dispersive transport of solutes results from the uneven distribution of water flow velocities within and

between different soil pores (Figure 22.3). Dispersion can be derived from Newton’s law of viscosity which

states that velocities within a single capillary tube follow a parabolic distribution, with the largest velocity

in the middle of the pore and zero velocities at the walls (Figure 22.3a). Solutes in the middle of a pore,

for this reason, will travel faster than solutes that are farther from the center. As the distribution of solute

ions within a pore depends on their charge, as well as on the charge of pore walls, some solutes may move

significantly faster than others. In some situations (i.e., for negatively charged anions in fine-textured

soils, leading to anion exclusion), the solute may even travel faster than the average velocity of water (e.g.,

Nielsen et al., 1986). Using Poiseuille’s law, one can further show that velocities in a capillary tube depend

strongly on the radius of the tube, and that the average velocity increases with the radius to the second

power. As soils consist of pores of many different radii, solute fluxes in pores of different radii will be

significantly different, with some solutes again traveling faster than others (Figure 22.3b).

The above pore-scale dispersion processes lead to an overall (macroscopic) hydrodynamic dispersion

process that mathematically can be described using Fick’s law in the same way as molecular diffusion,

that is,

Jlh = −θDlh

∂c

∂z = −θλv

∂c

∂z = −λq

∂c

∂z (22.17)
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 10 — #10

22-10 The Handbook of Groundwater Engineering

R

0

–R

C = C0 C = 0

Velocity

(a) (b)

FIGURE 22.3 Distribution of velocities in a single capillary (a) and distribution of velocities in a more complex pore

system (b).

where Dlh is the hydrodynamic dispersion coefficient [L2T−1], v is the average pore-water velocity [LT−1],

and λ is the dispersivity [L]. The dispersion coefficient in one-dimensional systems has been found to be

approximately proportional to the average pore-water velocity, with the proportionality constant generally

referred as the (longitudinal) dispersivity (Biggar and Nielsen, 1967). The discussion above holds for onedimensional transport; multi-dimensional applications require a more complicated dispersion tensor

involving longitudinal and transverse dispersivities (e.g., Bear, 1972).

Dispersivity is a transport property that is relatively difficult to measure experimentally. Estimates are

usually obtained by fitting measured breakthrough curves with analytical solutions of the advection–

dispersion equation (discussed further below). The dispersivity often changes with the distance over

which solute travels. Values of the longitudinal dispersivity typically range from about 1 cm for packed

laboratory columns, to about 5 or 10 cm for field soils. Longitudinal dispersivities can be significantly

larger (on the order of hundreds of meters) for regional groundwater transport problems (Gelhar et al.,

1985). If no other information is available, a good first approximation is to use a value of one-tenth of the

transport distance for the longitudinal dispersivity (e.g., Anderson, 1984), and a value of one-hundreds

of the transport distance for the transverse dispersivity.

22.3.1.3 Advection

Advective transport refers to solute being transported with the moving fluid, either in the liquid phase

(Jlc) or the gas phase (Jgc), that is,

Jlc = qc

Jgc = Jgg

(22.18)

where Jg is the gaseous flux density [LT−1]. Advective transport in the gaseous phase is often neglected as

its contribution in many applications is negligible compared to gaseous diffusion.

The total solute flux density in both the liquid and gaseous phases is obtained by incorporating

contributions from the various transport processes into Equation 22.14 to obtain

Jl = qc − θDs

l

∂c

∂z − θDlh

∂c

∂z = qc − θDe

∂c

∂z

Jg = −aDs

g

∂g

∂z

(22.19)

where De is the effective dispersion coefficient [L2T−1] that accounts for both diffusion and hydrodynamic

dispersion. Dispersion in most subsurface transport problems dominates molecular diffusion in the liquid

phase, except when the fluid velocity becomes relatively small or negligible. Notice that Equation 22.19

neglects advective and dispersive transport in the gaseous phase.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 11 — #11

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-11

22.3.2 Advection–Dispersion Equations

22.3.2.1 Transport Equations

The equation governing transport of dissolved solutes in the vadose zone is obtained by combining the

solute mass balance (Equation 22.10) with equations defining the total concentration of the chemical

(Equation 22.11) and the solute flux density (Equation 22.19) to give

∂(ρbs + θc + ag )

∂t = ∂∂xi



θDeij

∂c

∂xj



+

∂

∂xi



aDs

gij

∂g

∂xj



− ∂



qic



∂xi

− φ (22.20)

Notice that this equation is again written for multidimensional transport, and that Deij and Ds

ij are thus

components of the effective dispersion tensor in the liquid phase and a diffusion tensor in the gaseous

phase [L2T−1], respectively.

Many different variants of Equation 22.20 can be found in the literature. For example, for

one-dimensional transport of nonvolatile solutes, the equation simplifies to

∂(ρbs + θc)

∂t = ∂(θRc)∂t = ∂∂z



θDe

∂c

∂z



− ∂(qc)

∂z − φ (22.21)

where q is the vertical water flux density [LT−1] and R is the retardation factor [−]

R = 1 + ρb

θ

ds(c)

dc (22.22)

For transport of inert, nonadsorbing solutes during steady-state water flow we obtain

∂c

∂t = De

∂2c

∂z2 − v

∂c

∂z (22.23)

The above equations are usually referred to as advection–dispersion equations (ADEs).

22.3.2.2 Linear and Nonlinear Sorption

The ADE given by Equation 22.20 contains three unknown concentrations (those for the liquid, solid,

and gaseous phases), while Equation 22.21 contains two unknowns. To be able to solve these equations,

additional information is needed that somehow relates these concentrations to each other. The most

common way is to assume instantaneous sorption and to use adsorption isotherms to relate the liquid and

adsorbed concentrations. The simplest form of the adsorption isotherm is the linear isotherm given by

s = Kdc (22.24)

where Kd is the distribution coefficient [L3M−1]. One may verify that substitution of this equation into

Equation 22.22 leads to a constant value for the retardation factor (i.e., R = 1 + ρbKd/θ).

While the use of a linear isotherm greatly simplifies the mathematical description of solute transport,

sorption and exchange are generally nonlinear and most often depend also on the presence of competing

species in the soil solution. The solute retardation factor for nonlinear adsorption is not constant, as is

the case for linear adsorption, but changes as a function of concentration. Many models have been used

in the past to describe nonlinear sorption. The most commonly used nonlinear sorption models are those

by Freundlich (1909) and Langmuir (1918) given by:

s = Kfcβ (22.25)

s = Kdc

1 + ηc (22.26)
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 12 — #12

22-12 The Handbook of Groundwater Engineering

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Dissolved concentration [M/L3]

0 0.5 1 1.5 0 0.5 1 1.5

Sorbed concentration [M/M]

0.5

0.75

1

1.5

2

2.5

3

5

Increasing

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Dissolved concentration [M/L3]

Sorbed concentration [M/M]

0

0.5

0.75

1

1.5

2

2.5

3

5

Increasing

(a) (b)

FIGURE 22.4 Plots of the Freundlich adsorption isotherm given by (22.25), with Kd = 1 and β given in the

caption (a), and the Langmuir adsorption isotherm given by (22.26), with Kd = 1 and η given in the caption (b).

TABLE 22.1 Equilibrium Adsorption Equations (adapted from van Genuchten and Šim  ̊unek (1996))

Equation Model Reference

s = k1c + k2 Linear Lapidus and Amundson (1952)

Lindstrom et al. (1967)

s = k1ck3

1 + k2ck3 Freundlich–Langmuir Sips (1950), Šim  ̊unek et al. (1994, 2005)

s = k1c

1 + k2c +

k3c

1 + k4c Double Langmuir Shapiro and Fried (1959)

s = k1cck2/k3 Extended Freundlich Sibbesen (1981)

s = k1c

1 + k2c + k3

√c Gunary Gunary (1970)

s = k1ck2 − k3 Fitter–Sutton Fitter and Sutton (1975)

s = k1{1 − [1 + k2ck3 ]

k4 } Barry Barry (1992)

s = RT

k1

ln(k2c) Temkin Bache and Williams (1971)

s = k1c exp(−2k2s) Lindstrom et al. (1971)

van Genuchten et al. (1974) s

sT

= c[c + k1(cT − c) exp{k2(cT − 2c)}]−1 Modified Kielland Lai and Jurinak (1971)

k1, k2, k3, k4: empirical constants; R: universal gas constant; T: absolute temperature; cT: maximum solute

concentration; sT : maximum adsorbed concentration.

Source: Adapted from van Genuchten, M.Th. and Šimunek, J. In P.E. Rijtema and V. Eliáš (Eds.), Regional

Approaches to Water Pollution in the Environment, NATO ASI Series: 2. Environment. Kluwer, Dordrecht, The

Netherlands, pp. 139–172, 1996.

respectively, where Kf [M−β L−3β] and β [−] are coefficients in the Freundlich isotherm, and η [L3M−1]

is a coefficient in the Langmuir isotherm. Examples of linear, Freundlich and Langmuir adsorption

isotherms are given in Figure 22.4. Table 22.1 lists a range of linear and other sorption models frequently

used in solute transport studies.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 13 — #13

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-13

22.3.2.3 Volatilization

Volatilization is increasingly recognized as an important process affecting the fate of many organic chemicals, including pesticides, fumigants, and explosives in field soils (Jury et al., 1983, 1984; Glotfelty and

Schomburg, 1989). While many organic pollutants dissipate by means of chemical and microbiological

degradation, volatilization may be equally important for volatile substances, such as certain pesticides.

The volatility of pesticides is influenced by many factors, including the physicochemical properties of the

chemical itself as well-such environmental variables as temperature and solar energy. Even though only a

small fraction of a pesticide may exist in the gas phase, air-phase diffusion rates can sometimes be comparable to liquid-phase diffusion as gas-phase diffusion is about four orders of magnitude greater than

liquid phase diffusion. The importance of gaseous diffusion relative to other transport processes depends

also on the climate. For example, while transport of MTBE (gasoline oxygenate) is generally dominated

by liquid advection in humid areas, gaseous diffusion may be equally or more important in arid climates;

this even though only about 2% of MTBE may be in the gas phase.

The general transport equation given by Equation 22.20 can be simplified considerably when assuming

linear equilibrium sorption and volatilization such that the adsorbed (s) and gaseous (g ) concentrations

are linearly related to the solution concentration (c) through the distribution coefficients, Kd in (22.24)

and KH, that is,

g = KHc (22.27)

respectively, where KH is the dimensionless Henry’s constant [−]. Equation 22.20 for one-dimensional

transport then has the form:

∂(ρbKd + θ + aKH)c

∂t = ∂∂z



θDe

∂c

∂z



+

∂

∂z



aDs

gKH

∂c

∂z



− ∂



qc

∂x − φ (22.28)

or

∂θRc

∂t = ∂∂z



θDE

∂c

∂z



− ∂



qc

∂x − φ (22.29)

where the retardation factor R [−] and the effective dispersion coefficient DE [L2T−1] are defined as

follows:

R = 1 + ρbKd + aKH

θ

DE = De +

aDs

gKH

θ

(22.30)

Jury et al. (1983, 1984) provided for many organic chemicals their distribution coefficients Kd, Henry’s

constants KH, and calculated percent mass present in each phase.

22.3.3 Nonequilibrium Transport

As equilibrium solute transport models often fail to describe experimental data, a large number of

diffusion-controlled physical nonequilibrium and chemical-kinetic models have been proposed and used

to describe the transport of both non-adsorbing and adsorbing chemicals. Attempts to model nonequilibrium transport usually involve relatively simple first-order rate equations. Nonequilibrium models have

used the assumptions of two-region (dual-porosity) type transport involving solute exchange between

mobile and immobile liquid transport regions, and one-, two- or multi-site sorption formulations (e.g.,

Nielsen et al., 1986; Brusseau, 1999). Models simulating the transport of particle-type pollutants, such

as colloids, viruses, and bacteria, often also use first-order rate equations to describe such processes as

attachment, detachment, and straining. Nonequilibrium models generally have resulted in better descriptions of observed laboratory and field transport data, in part by providing additional degrees of freedom

for fitting observed concentration distributions.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 14 — #14

22-14 The Handbook of Groundwater Engineering

22.3.3.1 Physical Nonequilibrium

22.3.3.1.1 Dual-Porosity and Mobile–Immobile Water Models

The two-region transport model (Figures 22.2b and Figure 22.2c) assumes that the liquid phase can be

partitioned into distinct mobile (flowing) and immobile (stagnant) liquid pore regions, and that solute

exchange between the two liquid regions can be modeled as a first-order exchange process. Using the same

notation as before, the two-region solute transport model is given by (van Genuchten and Wagenet, 1989;

Toride et al., 1993):

∂θmocmo

∂t

+

∂f ρsmo

∂t = ∂

∂z



θmoDmo

∂cmo

∂z



− ∂qcmo

∂z − φmo − s

∂θimcim

∂t

+ ∂(1 − f )ρsim

∂t = −φim + s

(22.31)

for the mobile (macropores, subscript mo) and immobile (matrix, subscript im) domains, respectively,

where f is the dimensionless fraction of sorption sites in contact with the mobile water [−], φmo and

φim are reactions in the mobile and immobile domains [ML3T−1], respectively, and s is the solute

transfer rate between the two regions [ML3T−1]. Notice that the same equations Equation 22.31 can be

used to describe solute transport using both the mobile–immobile and dual-porosity models shown in

Figures 22.2b and Figure 22.2c, respectively.

22.3.3.1.2 Dual-Permeability Model

Analogous to Equations 22.6 and Equation 22.7 for water flow, the dual-permeability formulation for

solute transport can be based on advection–dispersion type equations for transport in both the fracture

and matrix regions as follows (Gerke and van Genuchten, 1993a):

∂θfcf

∂t

+ ∂ρsf

∂t = ∂∂z



θf Df

∂cf

∂z



− ∂qfcf

∂z − φf − sw (22.32)

∂θmcm

∂t

+ ∂ρsm

∂t = ∂∂z



θmDm

∂cm

∂z



− ∂qmcm

∂z − φm − s1 − w (22.33)

where the subscript f and m refer to the macroporous (fracture) and matrix pore systems, respectively; φf

and φm represent sources or sinks in the macroporous and matrix domains [ML3T−1], respectively; and

w is the ratio of the volumes of the macropore-domain (inter-aggregate) and the total soil systems [−].

Equation 22.32 and Equation 22.33 assume complete advective–dispersive type transport descriptions for

both the fractures and the matrix. Several authors simplified transport in the macropore domain, for

example, by considering only piston displacement of solutes (Ahuja and Hebson, 1992; Jarvis, 1994).

22.3.3.1.3 Mass Transfer

The transfer rate, s, in Equation 22.31 for solutes between the mobile and immobile domains in the

dual-porosity models can be given as the sum of diffusive and advective fluxes, and can be written as

s = αs(cmo − cim) + wc∗ (22.34)

where c∗ is equal to cmo for w > 0 and cim for w < 0, and αs is the first-order solute mass transfer

coefficient [T−1]. Notice that the advection term of Equation 22.34 is equal to zero for the mobile–

immobile model (Figure 22.2b) as the immobile water content in this model is assumed to be constant.

However, w may have a nonzero value in the dual-porosity model depicted in Figure 22.2c.

The transfer rate, s, in Equation 22.32 and Equation 22.33 for solutes between the fracture and matrix

regions is also usually given as the sum of diffusive and advective fluxes as follows (e.g., Gerke and van

Genuchten, 1996):

s = αs(1 − wm)(cf − cm) + wc∗ (22.35)
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 15 — #15

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-15

in which the mass transfer coefficient, αs [T−1], is of the form:

αs = β

d2 Da (22.36)

where Da is an effective diffusion coefficient [L2T−1] representing the diffusion properties of the fracture–

matrix interface.

Still more sophisticated models for physical nonequilibrium transport may be formulated. For example,

Pot et al. (2005) and Köhne et al. (2006) considered a dual-permeability model that divides the matrix

domain further into mobile and immobile subregions and used this model successfully to simulate bromide

transport in laboratory soil columns at different flow rates or for transient flow conditions, respectively.

22.3.3.2 Chemical Nonequilibrium

22.3.3.2.1 Kinetic Sorption Models

An alternative to expressing sorption as an instantaneous process using algebraic equations (e.g., Equation 22.24, Equation 22.25 or Equation 22.26) is to describe the kinetics of the reaction using ordinary

differential equations. The most popular and simplest formulation of a chemically controlled kinetic

reaction arises when first-order linear kinetics is assumed:

∂s

∂t = αk(Kdc − s) (22.37)

where αk is a first-order kinetic rate coefficient [T−1]. Several other nonequilibrium adsorption expressions were also used in the past (see Table 2 in van Genuchten and Šim  ̊unek, 1996). Models based on this

and other kinetic expressions are often referred to as one-site sorption models.

As transport models assuming chemically controlled nonequilibrium (one-site sorption) generally did

not result in significant improvements in their predictive capabilities when used to analyze laboratory

column experiments, the one-site first-order kinetic model was further expanded into a two-site sorption

concept that divides the available sorption sites into two fractions (Selim et al., 1976; van Genuchten and

Wagenet, 1989). In this approach, sorption on one fraction (type-1 sites) is assumed to be instantaneous

while sorption on the remaining (type-2) sites is considered to be time-dependent. Assuming a linear

sorption process, the two-site transport model is given by (van Genuchten and Wagenet, 1989)

∂(f ρbKd + θ )c

∂t = ∂∂z



θDe

∂ c

∂z



− ∂



qc

∂z − φe

∂sk

∂t = αk[(1 − f )Kdc − sk] − φk

(22.38)

where f is the fraction of exchange sites assumed to be at equilibrium [−], φe [ML3T−1] and φk

[MM−1T−1] are reactions in the equilibrium and nonequilibrium phases, respectively, and the subscript k refers to kinetic (type-2) sorption sites. Note that if f = 0, the two-site sorption model reduces

to the one-site fully kinetic sorption model (i.e., when only type-2 kinetic sites are present). On the other

hand, if f = 1, the two-site sorption model reduces to the equilibrium sorption model for which only

type-1 equilibrium sites are present.

22.3.3.2.2 Attachment/Detachment Models

Additionally, transport equations may include provisions for kinetic attachment/detachment of solutes

to the solid phase, thus permitting simulations of the transport of colloids, viruses, and bacteria. The

transport of these constituents is generally more complex than that of other solutes in that they are

affected by such additional processes as filtration, straining, sedimentation, adsorption and desorption,

growth, and inactivation. Virus, colloid, and bacteria transport and fate models commonly employ a
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 16 — #16

22-16 The Handbook of Groundwater Engineering

modified form of the ADE, in which the kinetic sorption equations are replaced with equations describing

kinetics of colloid attachment and detachment as follows:

ρ

∂s

∂t = θkaψc − kdρs (22.39)

where c is the (colloid, virus, bacteria) concentration in the aqueous phase [NcL−3], s is the solid phase

(colloid, virus, bacteria) concentration [NcM−1], in which Nc is a number of (colloid) particles, ka is

the first-order deposition (attachment) coefficient [T−1], kd is the first-order entrainment (detachment)

coefficient [T−1], andψ is a dimensionless colloid retention function [−]. The attachment and detachment

coefficients in Equation 22.39 have been found to strongly depend upon water content, with attachment

significantly increasing as the water content decreases.

To simulate reductions in the attachment coefficient due to filling of favorable sorption sites, ψ is sometimes assumed to decrease with increasing colloid mass retention. A Langmuirian dynamics (Adamczyk

et al., 1994) equation has been proposed for ψ to describe this blocking phenomenon:

ψ = smax − s

smax

= 1 − s

smax

(22.40)

in which smax is the maximum solid phase concentration [NcM−1].

A similar equation as Equation 22.39 was used by Bradford et al. (2003, 2004) to simulate the process of

pore straining. Bradford et al. (2003, 2004) hypothesized that the influence of straining and attachment

processes on colloid retention should be separated into two distinct components. They suggested the

following depth-dependent blocking coefficient for the straining process:

ψ =

dc + z − z0

dc

−β

(22.41)

where dc is the diameter of the sand grains [L], z0 is the coordinate of the location where the straining

process starts [L] (the surface of the soil profile, or interface between soil layers), and β is an empirical

factor (Bradford et al., 2003) [−].

The attachment coefficient is often 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, and gravitational sedimentation (Rajagopalan and Tien,

1976; Logan et al., 1995):

ka = 3(1 − θ )

2dc

ηαv (22.42)

where dc is the diameter of the sand grains [L], α is the sticking efficiency (ratio of the rate of particles

that stick to a collector to the rate they strike the collector) [−], v is the pore-water velocity [LT−1], and η

is the single-collector efficiency [−].

In related studies, Schijven and Hassanizadeh (2000) and Schijven and Šim  ̊unek (2002) used a twosite sorption model based on two equations (22.39) to successfully describe virus transport at both the

laboratory and field scale. Their model assumed that the sorption sites on the solid phase can be divided

into two fractions with different properties and various attachment and detachment rate coefficients.

22.3.3.3 Colloid-Facilitated Solute Transport

There is an increasing evidence that many contaminants, including radionuclides (Von Gunten et al.,

1988; Noell et al., 1998), pesticides (Vinten et al., 1983; Kan and Tomson, 1990, Lindqvist and Enfield,

1992), heavy metals (Grolimund et al., 1996), viruses, pharmaceuticals (Tolls, 2001; Thiele-Bruhn, 2003),

hormones (Hanselman et al., 2003), and other contaminants (Magee et al., 1991; Mansfeldt et al., 2004)

are transported in the subsurface not only with moving water, but also sorbed to mobile colloids. As many
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 17 — #17

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-17

colloids and microbes are negatively charged and thus electrostatically repelled by negatively charged solid

surfaces, which may lead to an anion exclusion process, their transport may be slightly enhanced 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). The transport of contaminants sorbed

to mobile colloids can thus significantly accelerate their transport relative to more standard advection–

transport descriptions.

Colloid-facilitated transport is a relatively complicated process that requires knowledge of water flow,

colloid transport, dissolved contaminant transport, and colloid-contaminant interaction. Transport and

mass balance equations, hence, must be formulated not only for water flow and colloid transport, but

also for the total contaminant, for contaminant sorbed kinetically or instantaneously to the solid phase,

and for contaminant sorbed to mobile colloids, to colloids attached to the soil solid phase, and to colloids

accumulating at the air–water interface. Development of such a model is beyond the scope of this chapter.

We refer interested readers to several manuscripts dealing with this topic: Mills et al. (1991), Corapcioglu

and Jiang (1993), Corapcioglu and Kim (1995), Jiang and Corapcioglu (1993), Noell et al. (1998), Saiers

et al. (1996), Saiers and Hornberger (1996), van de Weerd et al. (1998), and van Genuchten and Šim  ̊unek

(2004).

22.3.4 Stochastic Models

Much evidence suggests that solutions of classical solute transport models, no matter how refined to

include the most relevant chemical and microbiological processes and soil properties, often still fail to

accurately describe transport processes in most natural field soils. A major reason for this failure is the fact

that the subsurface environment is overwhelmingly heterogeneous. Heterogeneity occurs at a hierarchy of

spatial and time scales (Wheatcraft and Cushman, 1991), ranging from microscopic scales involving timedependent chemical sorption and precipitation/dissolution reactions, to intermediate scales involving the

preferential movement of water and chemicals through macropores or fractures, and to much larger scales

involving the spatial variability of soils across the landscape. Subsurface heterogeneity can be addressed

in terms of process-based descriptions which attempt to consider the effects of heterogeneity at one or

several scales. It can also be addressed using stochastic approaches that incorporate certain assumptions

about the transport process in the heterogeneous system (e.g., Sposito and Barry, 1987; Dagan, 1989).

In this section we briefly review several stochastic transport approaches, notably those using stream tube

models and the transfer function approach.

22.3.4.1 Stream Tube Models

The downward movement of chemicals from the soil surface to an underlying aquifer may be described

stochastically by viewing the field as a series of independent vertical columns, often referred to as “stream

tubes” (Figure 22.5), while solute mixing between the stream tubes is assumed to be negligible. Transport

V1

V2

V3

Vn–1

Vn

FIGURE 22.5 Schematic illustration of the stream tube model (Toride et al., 1995).
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 18 — #18

22-18 The Handbook of Groundwater Engineering

in each tube may be described deterministically with the standard ADE, or modifications thereof to include

additional geochemical and microbiological processes. Transport at the field scale is then implemented by

considering the column parameters as realizations of a stochastic process, having a random distribution

(Toride et al., 1995). Early examples are by Dagan and Bresler (1979) and Bresler and Dagan (1979) who

assumed that the saturated hydraulic conductivity had a lognormal distribution.

The stream tube model was implemented into the CXTFIT 2.0 program (Toride et al., 1995) for a

variety of transport scenarios in which the pore-water velocity in combination with either the dispersion

coefficient, De, the distribution coefficient for linear adsorption, Kd, or the first-order rate coefficient for

nonequilibrium adsorption, αk, are stochastic variables (Toride et al., 1995).

22.3.4.2 Transfer Function Models

Jury (1982) developed an alternative formulation for solute transport at the field scale, called the transfer function model. This model was developed based on two main assumptions about the soil system

(a) the solute transport is a linear process, and (b) the solute travel time probabilities do not change

over time. These two assumptions lead to the following transfer function equation that relates the solute

concentration at the outflow end of the system with the time-dependent solute input into the system:

cout(t) =

t

0

cin(t − t)f (t) dt (22.43)

The outflow at time t, cout(t) [ML−3], consists of the superposition of solute added at all times less than

t, cin(t − t

) [ML−3], weighted by its travel-time probability density function (pdf) f (t) [T−1] (Jury and

Horton, 2004). One important advantage of the transfer function approach is that it does not require

knowledge of the various transport processes within the flow domain. Different model distribution

functions can be used for the travel-time probability density function f (t) in Equation 22.43. Most

commonly used (Jury and Sposito, 1985) are the Fickian probability density function

f (t) = L

2

√

πDt 3 exp 

−(L − vt)2

4Dt 

(22.44)

and the lognormal distribution

f (t) = 1

√2πσt

exp 

(ln t − μ)2

2σ2



(22.45)

where D[L2T−1], v[LT−1],μ, and σ are model parameters, and L is the distance from the inflow boundary

to the outflow boundary [L]. To accommodate conditions when the water flux through the soil system is

not constant, Jury (1982) expressed the travel-time pdf as a function of the cumulative net applied water I:

I =

t

0

q(t) dt (22.46)

leading to the following transfer function equation:

cout(I) =



I(t)

0

cin(I − I)f (I) dI (22.47)
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 19 — #19

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-19

22.3.5 Multicomponent Reactive Solute Transport

The various mathematical descriptions of solute transport presented thus far all considered solutes that

would move independently of other solutes in the subsurface. In reality, the transport of reactive contaminants is more often than not affected by many often interactive physico-chemical and biological processes.

Simulating these processes requires a more comprehensive approach that couples the physical processes of

water flow and advective–dispersive transport with a range of biogeochemical processes. The soil solution

is always a mixture of many ions that may be involved in mutually dependent chemical processes, such as

complexation reactions, cation exchange, precipitation–dissolution, sorption–desorption, volatilization,

redox reactions, and degradation, among other reactions (Šim  ̊unek and Valocchi, 2002). The transport

and transformation of many chemical contaminants is further mediated by subsurface aerobic or anaerobic bacteria. Bacteria catalyze redox reactions in which organic compounds (e.g., hydrocarbons) act as

the electron donor and inorganic substances (oxygen, nitrate, sulfate, or metal oxides) as the electron

acceptor. By catalyzing such reactions, bacteria gain energy and organic carbon to produce new biomass.

These and related processes can be simulated using integrated reactive transport codes that couple the

physical processes of water flow and advective–dispersive solute transport with a range of biogeochemical

processes. This section reviews various modeling approaches for such multicomponent transport systems.

22.3.5.1 Components and Reversible Chemical Reactions

Multi-species chemical equilibrium systems are generally defined in terms of components. Components

may be defined as a set of linearly independent chemical entities such that every species in the system

can be uniquely represented as the product of a reaction involving only these components (Westall et al.,

1976). As a typical example, the chemical species CaCO0

3

Ca2+ + CO2−

3 ⇔ CaCO0

3 (22.48)

consists of the two components Ca2+ and CO2−

3 .

Reversible chemical reaction processes in equilibrium systems are most often represented using mass

action laws that relate thermodynamic equilibrium constants to activities (the thermodynamic effective

concentration) of the reactants and products (Mangold and Tsang, 1991; Appelo and Postma, 1993;

Bethke, 1996). For example, the reaction

bB + cC ⇔ dD + eE (22.49)

where b and c are the number of moles of substances B and C that react to yield d and e moles of products

D and E, is represented at equilibrium by the law of mass action

K = ad

D ae

E

ab

B acC

(22.50)

where K is a temperature-dependent thermodynamic equilibrium constant, and ai is the ion activity,

being defined as the product of the activity coefficient ( γi) and the ion molality (mi), that is, ai = γimi.

Single-ion activity coefficients may be calculated using either the Davies equation, an extended version of

the Debye–Hückel equation (Truesdell and Jones, 1974), or by means of Pitzer expressions (Pitzer, 1979).

Equation 22.50 can be used to describe all of the major chemical processes, such as aqueous complexation,

sorption, precipitation–dissolution, and acid–base and redox reactions, provided that the local chemical

equilibrium assumption is valid (Šim  ̊unek and Valocchi, 2002).
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 20 — #20

22-20 The Handbook of Groundwater Engineering

22.3.5.2 Complexation

Equations for aqueous complexation reactions can be obtained using the law of mass action as follows

(e.g., Yeh and Tripathi, 1990; Lichtner, 1996):

xi = Kx

i

γ x

i



Na

k=1

(γ a

k ck )

ax

ik i = 1, 2, ... , Mx (22.51)

where xi is the concentration of the ith complexed species, Kx

i is the thermodynamic equilibrium constant

of the ith complexed species, γ x

i is the activity coefficient of the ith complexed species, Na is the number of

aqueous components, ck is the concentration of the kth aqueous component, γ a

k is the activity coefficient

of the kth aqueous component species, ax

ik is the stoichiometric coefficient of the kth aqueous component

in the ith complexed species, Mx is the number of complexed species, and subscripts and superscripts x

and a refer to complexed species and aqueous components, respectively.

22.3.5.3 Precipitation and Dissolution

Equations describing precipitation–dissolution reactions are also obtained using the law of mass action,

but contrary to the other processes are represented by inequalities rather than equalities, as follows

(Šim  ̊unek and Valocchi, 2002):

Kp

i ≥ Qpi = 

Na

k=1

(γ a

k ck )

a

p

ik i = 1, 2, ... , Mp (22.52)

where Mp is the number of precipitated species, Kp

i is the thermodynamic equilibrium constant of the

ith precipitated species, that is, the solubility product equilibrium constant, Qp

i is the ion activity product

of the ith precipitated species, and a

p

ik is the stoichiometric coefficient of the kth aqueous component

in the ith precipitated species. The inequality in Equation 22.52 indicates that a particular precipate is

formed only when the solution is supersaturated with respect to its aqueous components. If the solution

is undersaturated then the precipitated species (if it exists) will dissolve to reach equilibrium conditions.

22.3.5.4 Cation Exchange

Partitioning between the solid exchange phase and the solution phase can be described with the general

exchange equation (White and Zelazny, 1986):

zj · Azi+Xzi + zi · Bzj+ ⇔ zi · Bzj+Xzj + zj · Azi+ (22.53)

where A and B are chemical formulas for particular cation (e.g., Ca2+ or Na+), X refers to an “exchanger”

site on the soil, and zi is the valence of species. The mass action equation resulting from this exchange

reaction is

Kij =



 c

zj+

j

γ a

j c

zj+

j





zi 



γ a

i c

z+

i

i

c

z+

i

i





zj

(22.54)

where Kij is the selectivity coefficient, and ci is the exchanger-phase concentration of the ith component

(expressed in moles per mass of solid).

22.3.5.5 Coupled System of Equations

Once the various chemical reactions are defined, the final system of governing equations usually consists

of several partial differential equations for solute transport (i.e., ADEs for each component) plus a set of

nonlinear algebraic and ordinary differential equations describing the equilibrium and kinetic reactions,

respectively. Each chemical and biological reaction must be represented by the corresponding algebraic

or ordinary differential equations depending upon the rate of the reaction. Since the reaction of one
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 21 — #21

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-21

species depends upon the concentration of many other species, the final sets of equations typically are

tightly coupled. For complex geochemical systems, consisting of many components and multidimensional

transport, numerical solution of these coupled equations is challenging (Šim  ̊unek and Valocchi, 2002). As

an alternative, more general models have recently been developed that also more loosely couple transport

and chemistry using a variety of sequential iterative or non-iterative operator-splitting approaches (e.g.,

Bell and Binning, 2004; Jacques and Simunek, 2005). Models based on these various approaches are

further discussed in Section 22.5b.

22.3.6 Multiphase Flow and Transport

While the transport of solutes in variably saturated media generally involves two phases (i.e., the liquid

and gaseous phases, with advection in the gaseous phase often being neglected), many contamination

problems also increasingly involve nonaqueous phase liquids (NAPLs) that are often only slightly miscible

with water. Nonaqueous phase liquids may consist of single organic compounds such as many industrial

solvents, or of a mixture of organic compounds such as gasoline and diesel fuel. Some of these compounds

can be denser than water (commonly referred to as DNAPLs) or lighter than water (LNAPLs). Their fate

and dynamics in the subsurface is affected by a multitude of compound-specific flow and multicomponent

transport processes, including interphase mass transfer and exchange (also with the solid phase).

Multiphase fluid flow models generally require flow equations for each fluid phase (water, air, NAPL).

Two-phase air–water systems hence could be modeled also using separate equations for air and water.

This shows that the standard Richards Equation 22.3 is a simplification of a more complete multiphase

(air–water) approach in that the air phase is assumed to have a negligible effect on variably saturated

flow, and that the air pressure varies only little in space and time. This assumption appears adequate for

most variably saturated flow problems. Similar assumptions, however, are generally not possible when

NAPLs are present. Hence mathematical descriptions of multiphase flow and transport in general require

flow equations for each of the three fluid phases, mass transport equations for all organic components

(including those associated with the solid phase), and appropriate equations to account for interphase

mass transfer processes. We refer readers to reviews by Abriola et al. (1999) and Rathfelder et al. (2000)

for discussions of the complexities involved in modeling such systems subject to multiphase flow, multicomponent transport, and interphase mass transfer. An excellent overview of a variety of experimental

approaches for measuring the physical and hydraulic properties of multi-fluid systems is given by Lenhard

et al. (2002).

22.3.7 Initial and Boundary Conditions

22.3.7.1 Initial Conditions

The governing equations for solute transport can be solved analytically or numerically provided that the

initial and boundary conditions are specified. Initial conditions need to be specified for each equilibrium

phase concentration, that is,

c(x, y,z, t) = ci(x, y,z, 0) (22.55)

where ci is the initial concentration [ML−3], as well as for all nonequilibrium phases such as concentrations

in the immobile region, sorbed concentrations associated with kinetic sites, and initially attached or

strained colloid concentrations.

22.3.7.2 Boundary Conditions

Complex interactions between the transport domain and its environment often must be considered for

the water flow part of the problem being considered since these interactions determine the magnitude of

water fluxes across the domain boundaries. By comparison, the solute transport part of most analytical

and numerical models usually considers only three types of boundary conditions. When the concentration
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 22 — #22

22-22 The Handbook of Groundwater Engineering

at the boundary is known, one can use a first-type (or Dirichlet type) boundary condition:

c(x, y,z, t) = c0(x, y,z, t) for (x, y,z) ∈ d (22.56)

where c0 is a prescribed concentration [ML−3] at or along the D Dirichlet boundary segments. This

boundary condition is often referred to as a concentration boundary condition. A third-type (Cauchy

type) boundary condition may be used to prescribe the concentration flux at the boundary as follows:

−θDij

∂c

∂xj

ni + qinic = qinic0 for (x,z) ∈ C (22.57)

in which qini represents the outward fluid flux [LT−1], ni is the outward unit normal vector and c0

is the concentration of the incoming fluid [ML−3]. In some cases, for example, when a boundary is

impermeable (q0=0) or when water flow is directed out of the region, the Cauchy boundary condition

reduces to a second-type (Neumann type) boundary condition of the form:

θDij

∂c

∂xj

ni = 0 for (x,z) ∈ N (22.58)

Most applications require a Cauchy boundary condition rather than Dirichlet (or concentration) boundary

condition. Since Cauchy boundary conditions define the solute flux across a boundary, the solute flux

entering the transport domain will be known exactly (as specified). This specified solute flux is then in

the transport domain divided into advective and dispersive components. On the other hand, Dirichlet

boundary condition controls only the concentration on the boundary, and not the solute flux which,

because of its advective and dispersive contributions, will be larger than for a Cauchy boundary condition.

The incorrect use of Dirichlet rather than Cauchy boundary conditions may lead to significant mass

balance errors at an earlier time, especially for relative short transport domains (van Genuchten and

Parker, 1984).

A different type of boundary condition is sometimes used for volatile solutes when they are present in

both the liquid and gas phases. This situation requires a third-type boundary condition, but modified to

include an additional volatilization term accounting for gaseous diffusion through a stagnant boundary

layer of thickness d [L] above the soil surface. The additional solute flux is often assumed to be proportional

to the difference in gas concentrations above and below this boundary layer (e.g., Jury et al., 1983). This

modified boundary condition has the form:

−θDij

∂c

∂xj

ni + qinic = qinic0 +

Dg

d (kgc − gatm) for (x,z) ∈ C (22.59)

where Dg is the molecular diffusion coefficient in the gas phase [L2T−1] and gatm is the gas concentration

above the stagnant boundary layer [ML−3]. We note that Jury et al. (1983) assumed gatm in Equation 22.59

to be zero.

Still other types of boundary conditions can be used. One example is the use of the Bateman equations

(Bateman, 1910) to account for a finite rate of release of multiple solutes that are subject to a first-order

sequential decay (e.g., radionuclides, nitrogen species). These solutes are typically released from a waste

site into the environment as a consequence of decay reactions in the waste site (e.g., van Genuchten, 1985).

22.4 Analytical Models

A large number of computer models using both analytical and numerical solutions of the flow and solute

transport equations are now available for a wide range of applications in research and management of

natural subsurface systems (Šim  ̊unek, 2005). Modeling approaches range from relatively simple analytical
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 23 — #23

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-23

and semi-analytical models, to more complex numerical codes that permit consideration of a large

number of simultaneous nonlinear processes. While for certain conditions (e.g., for linear sorption, a

concentration independent sink term φ , and a steady flow field) the solute transport equations (e.g.,

Equation 22.20, Equation 22.21, Equation 22.23, Equation 22.31, and Equation 22.38) become linear, the

governing flow equations (e.g., Equation 22.3) are generally highly nonlinear because of the nonlinear

dependency of the soil hydraulic properties on the pressure head or water content. Consequently, many

analytical solutions have been derived in the past for solute transport equations, and are now widely used

for analyzing solute transport during steady-state flow. Although a large number of analytical solutions

also exist for the unsaturated flow equation, they generally can be applied only to relatively simple flow

problems. The majority of applications for water flow in the vadose zone require a numerical solution of

the Richards equation.

Analytical methods are representative of the classical mathematical approach for solving differential

equations to produce an exact solution for a particular problem. Analytical models usually lead to an

explicit equation for the concentration (or the pressure head, water content, or temperature) at a particular time and location. Hence, one can evaluate the concentration directly without the time stepping

that is typical of numerical methods. While exceptions exist (e.g., Liu et al., 2000; Cotta et al., 2005),

analytical solutions usually can be derived only for simplified transport systems involving linearized governing equations, homogeneous soils, simplified geometries of the transport domain, and constant or

highly simplified initial and boundary conditions. Unfortunately, analytical solutions for more complex

situations, such as for transient water flow or nonequilibrium solute transport with nonlinear reactions,

are generally not available and cannot be derived, in which case numerical models must be employed

(Šim  ̊unek, 2005).

22.4.1 Analytical Approaches

Analytical solutions are usually obtained by applying various transformations (e.g., Laplace, Fourier or

other transforms) to the governing equations, invoking a separation of variables, and using the Green’s

function approach (e.g., Leij et al., 2000; Cotta et al., 2005).

22.4.1.1 Laplace Transformation

The Laplace transform, £, of the solute concentration with respect to time is defined as follows:

c ̄(x, s) = L[c(x, t)] = ∞

0

c(x, t) exp(−st) dt (22.60)

where s is the transform variable [T−1]. Laplace transforms can greatly simplify the governing solute

transport equations by eliminating one independent variable, usually time. They also transform the

original transport equation from a partial to an ordinary differential equation. Using Laplace transforms,

one thus obtains a governing ADE in the Laplace domain that is much simpler to solve analytically than

the original equation. The analytical solution in the Laplace space is subsequently inverted to the real

space using either a table of Laplace transforms, by applying inversion theorems, or by using a numerical

inversion program.

22.4.1.2 Fourier Transformation

Two- and three-dimensional problems often require not only Laplace transforms with respect to time

(transformation from t tos) and one spatial coordinate (usually x being transformed to r), but often also

a double Fourier transform with respect to two other spatial coordinates (y and z), which is given by (Leij
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 24 — #24

22-24 The Handbook of Groundwater Engineering

and Toride, 1997):

FYZ [ ̄c(r, y,z, s)]= ̃c(r, γ , κ, s) = 1

2π

∞

−∞

∞

−∞

c(r, y,z, s) exp(iγ y + iκz) dy dz (22.61)

where i

2 = −1 and γ and κ are transformation variables pertaining to the y and z coordinates. Similarly

as for Laplace transforms, Fourier transforms lead to an equation that can be solved analytically more

easily. One also needs to transform this analytical solution from the Laplace and Fourier domains back

to the original time and space domains. Several analytical solutions, including the derivation process, for

different initial and boundary conditions, and different transport domain geometries are given by Leij

and Toride (1997) in the N3DADE manual.

22.4.1.3 Method of Moments

Statistical moments are often used to characterize the statistical distribution of the solute concentration

versus time and space. For example, moments can be used to evaluate the average residence time or the

width of the residence time distribution. These attributes can be quantified with the moments of the

distribution (Skaggs and Leij, 2002). The pth time moment, mp, is defined as:

mp =

∞

0

t pc(t) dt (p = 0, 1, 2, . . .) (22.62)

where c(t) is the breakthrough curve measured at a certain location. The zeroth moment is related to

the mass of solute contained in the breakthrough curve, while the first moment is related to the mean

residence time. Normalized and central moments are defined as:

Mp = mp

m0

Mp = 1

m0

∞

0

(1 − M1)

pc(t) dt

(22.63)

The second central moment relates to the degree of solute spreading. Similar expressions can be written

also for the spatial (depth) moments characterizing the spatial solute distribution. Moments for particular

problems or processes can then be obtained by substituting appropriate analytical solutions in the above

equations.

22.4.2 Existing Models

Numerous analytical solutions of the linear advection–dispersion solute transport equations (e.g., Equation 22.21, Equation 22.23, Equation 22.31, and Equation 22.38), or their two- and three-dimensional

equivalents, have been developed during the last 40 years and are now widely used for predictive purposes

and for analyzing laboratory and field observed concentration distributions. The majority of these solutions pertain to solute transport equations assuming constant water content, θ, and flux, q, values (i.e.,

for steady-state water flow conditions in a homogeneous medium).

22.4.2.1 One-Dimensional Models

Some of the more popular one-dimensional analytical transport models have been CFITM (van

Genuchten, 1980b), CFITIM (van Genuchten, 1981), CXTFIT (Parker and van Genuchten, 1984), and

CXTFIT2 (Toride et al., 1995). While CFITM considers only one-dimensional equilibrium solute transport in both finite and infinite solute columns, CFITIM additionally considers physical and chemical
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 25 — #25

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-25

nonequilibrium transport (i.e., the two-region mobile–immobile model for physical nonequilibrium and

the two-site sorption model for chemical nonequilibrium). CXTFIT further expanded the capabilities of

CFITIM by considering more general initial and boundary conditions, as well as degradation processes.

CXTFIT2 (Toride et al., 1995), an updated version of CXTFIT, solves both direct and inverse problems for

three different one-dimensional transport models (1) the conventional advection–dispersion equation,

ADE; (2) the chemical and physical nonequilibrium ADEs; and (3) a stochastic stream tube model based

upon the local-scale equilibrium or nonequilibrium ADE. These three types of models all consider linear

adsorption, and include zero- and first-order decay/source terms. All of these models may be used also to

solve the inverse problem by fitting the analytical ADE solutions to experimental results.

22.4.2.2 Multi-Dimensional Models

Some of the more popular multi-dimensional analytical transport models have been AT123D (Yeh, 1981),

3DADE (Leij and Bradford, 1994), N3DADE (Leij and Toride, 1997), and MYGRT (Ungs et al., 1998). These

programs provide analytical solutions to transport problems in two- and three-dimensional domains.

3DADE also included parameter estimation capabilities.

A large number of analytical models for one-, two-, and three-dimensional solute transport problems

were recently incorporated into the public domain software package STANMOD (STudio of ANalytical

MODels) (Šim  ̊unek et al., 1999a) (http://www.hydrus2d.com). Figure 22.6 shows a STANMOD dialog

window where users can make a selection of the analytical program to be used. This Windows-based

computer software package includes not only programs for equilibrium advective-dispersive transport

(i.e., the CFITM of van Genuchten 1980b) for one-dimensional transport and 3DADE (Leij and Bradford,

1994) for three-dimensional problems (Figure 22.7), but also programs for more complex problems.

For example, STANMOD also incorporates the CFITIM (van Genuchten, 1981) and N3DADE (Leij

and Toride, 1997) programs for nonequilibrium transport (i.e., the two-region mobile–immobile model

for physical nonequilibrium and the two-site sorption model for chemical nonequilibrium) in one and

multiple dimensions, respectively.

FIGURE 22.6 Selection of an analytical model in the STANMOD software package (Šim  ̊unek et al., 1999a).
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 26 — #26

22-26 The Handbook of Groundwater Engineering

FIGURE 22.7 Selection of the geometry of the transport domain for the 3DADE code (Leij and Bradford, 1994)

within the STANMOD software package (Šim  ̊unek et al., 1999a).

Figure 22.7 shows a dialog window of STANMOD in which users can select the geometry of the

transport domain for use in the 3DADE program. STANMOD also includes CXTFIT2 (Toride et al.,

1995) as discussed earlier, and the CHAIN code of van Genuchten (1985) for analyzing the advective–

dispersive transport of up to four solutes involved in sequential first-order decay reactions. Examples of

the latter are the migration of radionuclides in which the chain members form first-order decay reactions,

and the simultaneous movement of various interacting nitrogen or organic chemicals. The latest version

of STANMOD also includes the screening model of Jury et al. (1983) for describing transport and

volatilization of soil-applied volatile organic chemicals.

An application of CFITIM within of the STANMOD program is demonstrated in Figure 22.8, which

shows measured and fitted breakthrough curves of a nonreactive (3H2O) solute for transport through

a saturated Glendale clay loam soil. In this example, a tritiated water pulse of 3.102 pore volumes was

applied to a 30 cm long column, with the breakthrough curve being determined from the effluent concentrations. An analytical solution for two-region (mobile–immobile) physical nonequilibrium transport

(van Genuchten, 1981) was used for the analysis. With the pore-water velocity known (v = 37.54 cm d−1)

and assuming that the retardation factor R for 3H2O was equal to one (no sorption or anion exclusion),

only three parameters were optimized against the breakthrough curve, that is, the dispersion coefficient

D (= 15.6 cm2 d−1), a dimensionless variable β (= 0.823) for partitioning the liquid phase in mobile and

immobile regions, and a dimensionless mass transfer coefficient ω (= 0.870).

An example application of the 3DADE program (Leij and Bradford, 1994) in the STANMOD package

is demonstrated in Figure 22.9, which shows calculated concentration contours for a transport problem

in which solute is applied at a concentration C0 = 1, assuming a third-type boundary condition, to a

rectangular area of the soil surface (or within groundwater) of 15 × 15 cm2. Other transport parameters

were as follows: v = 10 cm d−1, Dx = 100 and Dy = Dz = 10 cm2 d−1, and R = 1. Figure 22.9 shows

concentration contours in the xy-plane at t = 1 days, and z = 0.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 27 — #27

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-27

0

0246

0.2

0.4

0.6

0.8

1

1.2

8

Pore volumes [–]

Concentration [–]

Measured

Fitted

FIGURE 22.8 Breakthrough curve for nonreactive tritium transport as analyzed with the two-region physical

nonequilibrium program CXTFIT2 (Toride et al., 1995) within STANMOD (Šim  ̊unek et al., 1999a).

20

0

y

–20

0.0 0.2 0.4 0.6 0.8 1.0

0 x 30

FIGURE 22.9 Solute distribution in the xy-plane for continuous application at x = 0 and −7.5 < y < 7.5 cm after

1 days, with v = 10 cm d−1, Dx = 100 and Dy = Dz = 10 cm2 d−1 as calculated with 3DADE (Leij and Bradford,

1994) within STANMOD (Šim  ̊unek et al., 1999a).

22.5 Numerical Models

22.5.1 Numerical Approaches

While analytical and semi-analytical solutions are still popularly used for many relatively simple applications, the ever-increasing power of personal computers, and the development of more accurate and
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 28 — #28

22-28 The Handbook of Groundwater Engineering

stable numerical solution techniques, have led to the much wider use of numerical models in recent years.

Numerical methods are in general superior to analytical methods in terms of being able to solve much

more practical problems (Šim  ̊unek, 2005). They allow users to design complicated geometries that reflect

complex natural geologic and hydrologic conditions, control parameters in space and time, prescribe

more realistic initial and boundary conditions, and permit the implementation of nonlinear constitutive

relationships. Numerical methods usually subdivide the time and spatial coordinates into smaller pieces,

such as finite differences, finite elements, and finite volumes, and reformulate the continuous form of

governing partial differential equations in terms of a system of algebraic equations. To obtain solutions

at certain times, numerical methods generally require intermediate simulations (time-stepping) between

the initial condition and the points in time for which the solution is needed.

Reviews of the history of development of various numerical techniques used in vadose zone flow and

transport models are given by van Genuchten and Šim  ̊unek (1996) and Šim  ̊unek (2005). Here we will

summarize only the main principles behind some of the more popular techniques.

22.5.1.1 Finite Differences

The method of finite differences is generally very intuitive and relatively easy to implement. Time and

space are both divided into small increments t and z (or x and z) (Figure 22.10). Temporal and

spatial derivatives in the governing equation are then replaced with finite differences (formally using Taylor

series expansions). For example, the standard advection–dispersion equation for steady-state water flow:

∂c

∂t = −∂J∂z = D

∂2c

∂z2 − v

∂c

∂z (22.64)

can be approximated as follows using an explicit (forward-in-time) finite difference scheme:

c

j+1

i − c

j

i

t = −J

j

i+1/2 − J

j

i−1/2

z = D c

j

i+1 − 2c

j

i + c

j

i−1

(z)

2 − v

c

j

i+1 − c

j

i−1

2z (22.65)

where subscripts refer to spatial discretization and superscript to temporal discretization (j and j+1 are

for the previous and actual time levels, respectively; see Figure 22.10), t is the time step, and z is the

spatial step (assumed to be constant). Notice that this equation contains only one unknown variable (i.e.,

the concentration c

j+1

i at the new time level), which hence can be evaluated directly (explicitly) by solving

the equation.

zi+1 ci +1

ci ci

ci +1

ci – 1 ci – 1

j +1

j +1

j +1

j

j

j

t i+1 t

i

zi –1

zi

z

t

∆t

∆z

∆zi

∆x

∆xi–1 ∆xi

∆z

∆zj

x

zj

z

∆zi–1

zj – 1

zj+1

∆zj – 1

xi –1

ci –1, j +1

ci–1, j

ci–1, j –1

ci, j

ci +1, j +1

ci+1, j

ci, j–1 ci+1, j–1

ci, j +1

xi xi + 1 (a) (b)

FIGURE 22.10 Examples of the spatial and temporal finite difference discretization of a one-dimensional

problem (a), and the finite difference discretization of a two-dimensional domain (b).
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 29 — #29

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-29

By comparison, a fully implicit (backward-in-time) finite difference scheme can be written as follows:

c

j+1

i − c

j

i

t = −J

j+1

i+1/2 − J

j+1

i−1/2

z = D c

j+1

i+1 − 2c

j+1

i + c

j+1

i−1

(z)

2 − v

c

j+1

i+1 − c

j+1

i−1

2z (22.66)

and an implicit (weighted) finite difference scheme as:

c

j+1

i − c

j

i

t = D ε(c

j+1

i+1 − 2c

j+1

i + c

j+1

i−1 ) + (1 − ε)(c

j

i+1 − 2c

j

i + c

j

i−1)

(z)2

− v

ε(c

j+1

i+1 − c

j+1

i−1 ) + (1 − ε)(c

j

i+1 − c

j

i−1)

2z (22.67)

where ε is a temporal weighting coefficient. Different finite-difference schemes results depending upon

the value of ε, that is, an explicit scheme when ε = 0, a Crank–Nicholson time-centered scheme when

ε = 0.5, and a fully implicit scheme when ε = 1. Of these, the latter two cases lead to several unknown

concentrations in each equation at the new time level. Consequently, equations for all nodes of the

transport domain need to be assembled into an algebraic system of linear equations to produce a matrix

equation of the form

[P]

j+1{c}j+1 = {F} (22.68)

in which [P] is a tridiagonal matrix [P] given by:

[P] =

















































d1 e1 0 0

b2 d2 e2 0 0

0 b3 d3 e3 0 0

.

.

. .

.

. .

.

.

0 0 bN−2 dN−2 eN−2 0

0 0 bN−1 dN−1 eN−1

0 0 bN dN

















































(22.69)

Equation 22.69 is subsequently solved using a matrix equation solver to yield the concentrations at the

new time level.

The problem becomes more complex for problems involving transient flow. In that case one must

first obtain estimates of the water contents and fluxes at the new time level using a numerical solution

of the Richards equation. These water contents and fluxes are subsequently used for assembling the

matrix equations for solute transport. While finite differences are relatively straightforward to implement

for many transport problems, one major disadvantage is that this method cannot be easily applied to

relatively complicated (irregular) transport domains in two and three dimensions.

22.5.1.2 Finite Elements

Finite element methods can be implemented in very much the same way as finite differences for one-, two-,

and three-dimensional problems. A major advantage of the finite elements is that they are much easier

used to discretize complex two- and three-dimensional transport domains (Figure 22.11). As an example,

Figure 22.11 shows triangular unstructured finite element grids for a regular rectangular and an irregular

domain as generated with the automated MeshGen2D mesh generator of HYDRUS-2D (Šim  ̊unek et al.,

1999b). Note that even though the figure on the right (Figure 22.11) has an irregular soil surface, as well as

a tile drain within the transport domain, MeshGen2D could easily discretize/accommodate this transport

domain using an unstructured triangular finite element mesh.

The finite element method assumes that the dependent variable in the advection–dispersion equation

(e.g., the concentration function c(x, t) in (22.64)) can be approximated by a finite series c

(x, t) of the
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 30 — #30

22-30 The Handbook of Groundwater Engineering

Z

X

(a) (b)

FIGURE 22.11 Examples of the triangular finite element grids for regular (a) and irregular (b) two-dimensional

transport domains.

i-1 i i +1

1 2

∆xi

x

c

ci–1 ci+1

1

x

ci

i – 1

1

i i + 1 x

(a)

(b)

FIGURE 22.12 Example of one-dimensional linear basis functions (a) and their use to approximate solution c(x, t)

of the advection–dispersion equation (b).

form (Figure 22.12b):

c

(x, t) = 

N

m=1

φm(x)cm(t) (22.70)

where φm are the selected linear basis functions that fulfill the condition φm(xn) = δnm, δnm is Kronecker

delta ( δnm = 1 for m = n, and δnm = 0 for m = n), cm are the unknown time-dependent coefficients

which represent solutions of (22.64) at the finite element nodal points, and N is the total number of nodal

points. For example, for the one-dimensional finite element xi < x < xi+1, linear basis functions have
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 31 — #31

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-31

the following form (Figure 22.12a):

φ1 = 1 − x − xi

x

φ2 = x − xi

x







xi ≤ x ≤ xi+1 (22.71)

where x (= xi+1 − xi) is the size of a finite element [L], that is, the distance between two neighboring

nodal points. The approximate solution c

(x, t) converges to the correct solution c(x, t) as the number of

basis functions N increases.

Application of the Galerkin method, which postulates that the differential operator associated with the

transport equation is orthogonal to each of the N basis functions (e.g., Pinder and Gray, 1977), leads to

the following system of N time-dependent differential equations with N unknown values cn(t):



L

0



−∂c

∂t − D

∂2c

∂x2 − v

∂c

∂x



φn dx = 0 (22.72)

where φn are the selected linear weighting functions that can be the same as the basis functions φm.

Integrating by parts the terms containing the spatial derivatives to remove the second-order derivatives,

and substituting (22.70) for c(x, t) leads to the following equation:



e



Le

0



−∂cm

∂t φm



φn dx −

e



Le

0



Dcm

∂φm

∂x − vcmφm

∂φn

∂x dx − qsLφn(L) + qs0φn(0) = 0

(22.73)

where qs0 and qsL are solute fluxes across the left and right boundaries, respectively; e is the element index,

Le is the size of an element e, and L is the size of the transport domain (length of the soil profile). The

integrals in Equation 22.73 can be evaluated either numerically or analytically, depending upon the type

of invoked basis and weighting functions ( φm and φn, respectively). The equation can be rewritten in

matrix form similarly as (Equation 22.68) and then solved using a matrix equation solver. For a detailed

derivation of the matrix equation we refer to standard finite element textbooks (e.g., Pinder and Gray, 1977;

Zienkiewicz, 1977; Huyakorn and Pinder, 1983) or manuals of some of the models (e.g., the HYDRUS

models by Šim  ̊unek et al., 1998a,b).

22.5.1.3 Stability and Oscillations

Numerical solutions of the transport equation often exhibit oscillatory behavior and/or excessive numerical dispersion near relatively sharp concentration fronts (Huyakorn and Pinder, 1983; Šim  ̊unek et al.,

1998a). These problems can be especially serious for advection-dominated transport characterized by

small dispersivities. One way to partially circumvent numerical oscillations is to use upstream weighing

in which case the flux term of Equation 22.64 is not weighted using regular linear basis functions φn,

but instead with nonlinear quadratic functions φu

n . Undesired oscillations can often be prevented also by

selecting an appropriate combination of the space and time discretizations. Two dimensionless numbers

may be used to characterize the discretizations in space and time. One of these is the grid Peclet number, Pee, which defines the predominant type of solute transport (notably the ratio of the advective and

dispersive transport terms) in relation to coarseness of the finite element grid:

Pee = vx

D (22.74)
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 32 — #32

22-32 The Handbook of Groundwater Engineering

where x is the characteristic length of a finite element [L]. The Peclet number increases when the advective part of the transport equation dominates the dispersive part (i.e., when a relatively steep concentration

front is present). To achieve acceptable numerical results, the spatial discretization must be kept relatively

fine to maintain a low Peclet number. Numerical oscillations can be virtually eliminated when the local

Peclet numbers do not exceed about 5. However, acceptably small oscillations may be obtained with

local Peclet numbers as high as 10 (Huyakorn and Pinder, 1983). Undesired oscillations for higher Peclet

numbers can be effectively eliminated by using upstream weighing.

A second dimensionless number, the Courant number, Cre, may be used to characterize the relative

extent of numerical oscillations in the numerical solution. The Courant number is associated with the

time discretization, t [T], as follows:

Cre = vt

x (22.75)

Given a certain spatial discretization, the time step must be selected such that the Courant number remains

less than or equal to 1.

Perrochet and Berod (1993) developed a criterion for minimizing or eliminating numerical oscillations

that is based on the product of the Peclet and Courant numbers:

Pe · Cr ≤ ωs (=2) (22.76)

where ωs is the performance index [−]. This criterion indicates that advection-dominated transport

problems having large Pe numbers can be safely simulated provided Cr is reduced according to (22.76)

(Perrochet and Berod, 1993). When small oscillations in the solution can be tolerated, ωs can be increased

to about 5 or 10.

22.5.2 Existing Models

22.5.2.1 Single-Species Solute Transport Models

A large number of numerical models are now available for evaluating variably saturated water flow and

solute transport processes in the subsurface. Some of these models are in the public domain, such as

MACRO (Jarvis, 1994), SWAP (van Dam et al., 1997), UNSATH (Fayer, 2000), VS2DI (Healy, 1990), and

HYDRUS-1D (Šim  ̊unek et al., 1998a, 2005), while others are in the commercial domain, such as HYDRUS2D (Šim  ̊unek et al., 1999b) and MODFLOW-SURFACT (HydroGeoLogic, 1996). These models vary

widely in terms of their complexity, sophistication, and ease of use. Although some models are still being

run under the DOS operating system, with associated difficulties of preparing input files and interpreting

tabulated outputs, many others, especially those in the commercial domain, are supported by sophisticated

graphics-based interfaces that greatly simplify their use (Šim  ̊unek et al., 1998a; 1999b). Several studies

have recently reviewed and compared various numerical models for vadose zone applications (e.g., Wilson

et al., 1999; Scanlon et al., 2002; MDH Engineered Solutions Corp., 2003; Vanderborght et al., 2005). These

studies typically evaluated comparable precision, speed, and ease of use of the codes involved.

While many of the early numerical zone models were developed mostly for agricultural applications

(e.g., optimizing soil moisture conditions for crop production), this focus has increasingly shifted to

environmental research, with the primary concern now being the subsurface fate and transport of various agricultural and other contaminants, such as radionuclides, hydrocarbons, fumigants, pesticides,

nutrients, pathogens, pharmaceuticals, viruses, bacteria, colloids, and/or toxic trace elements. The earlier

models solved the governing flow and transport equations for relatively simplified system-independent

boundary conditions (i.e., specified pressure heads or fluxes, and free drainage). Models developed recently

can cope with much more complex system-dependent boundary conditions evaluating surface flow and

energy balances and accounting for the simultaneous movement of water, vapor, and heat. Examples are

DAISY (Hansen et al., 1990), TOUGH2 (Pruess, 1991), SHAW (Flerchinger et al., 1996), SWAP (van Dam

et al., 1997), HYDRUS-1D (Šim  ̊unek et al., 1998a, 2005), UNSATH (Fayer, 2000), and COUP (Jansson
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 33 — #33

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-33

and Karlberg, 2001). Several models now also account for the extremely nonlinear processes associated

with the freezing and thawing cycle (e.g., DAISY, SHAW, and COUP).

Models have recently also become increasingly sophisticated in terms of the type and complexity

of solute transport processes that can be simulated. Transport models are no longer being limited to

solutes undergoing relatively simple chemical reactions such as linear sorption and first-order decay,

but now consider also a variety of nonlinear sorption and exchange processes, physical and chemical

nonequilibrium transport, volatilization, gas diffusion, colloid attachment/ detachment, decay chain

reactions, and many other processes (e.g., the HYDRUS-1D and -2D codes of Šim  ̊unek et al., 1999b,

2005, or MODFLOW-SURFACT of HydroGeoLogic, Inc., 1996). For example, the general formulation of

the transport equations in the HYDRUS codes permit simulations of not only non-adsorbing or linearly

sorbing chemicals, but also of a variety of other contaminants, such a viruses (Schijven and Šim  ̊unek,

2001), colloids (Bradford et al., 2002), cadmium (Seuntjens et al., 2001), and hormones (Casey et al.,

2003, 2004), or chemicals involved in the sequential biodegradation of chlorinated aliphatic hydrocarbons

(Schaerlaekens et al., 1999; Casey and Šim  ̊unek, 2001).

Much effort has recently been directed also toward improving models for purposes of simulating

nonequilibrium and/or preferential flow. Examples of this are the TOUGH codes (Pruess, 1991, 2004),

MACRO (Jarvis, 1994), and HYDRUS-1D (Šim  ̊unek et al., 2003, 2005). These models typically assume the

presence of dual-porosity and dual-permeability regions, with different fluxes possible in the two regions.

As discussed in detail in Section 22.2b, dual-porosity, mobile–immobile water flow models result when the

Richards or kinematic wave equations are used for flow in the fractures, and immobile water is assumed

to exist in the matrix, whereas dual-permeability models assume that water is mobile in both the matrix

and fracture domains, while invoking terms accounting for the exchange of water and solutes between

the two regions. Example applications of these dual-porosity and dual-permeability models are given by

Šim  ̊unek et al. (2001), Haws et al. (2005), Köhne et al. (2004b, 2005), and Pot et al. (2005), among many

others.

As an example of available vadose zone flow and transport models, we discuss here in somewhat greater

detail the HYDRUS software packages of Šim  ̊unek et al. (1999b, 2005) and the MODFLOW-SURFACT

model developed by HydroGeoLogic (1996).

22.5.2.1.1 The HYDRUS Software Packages

HYDRUS-1D (Šim  ̊unek et al., 2005) and HYDRUS-2D (Šim  ̊unek et al., 1999b) are software packages

(www.hydrus2d.com) that simulate the one- and two-dimensional movement of water, heat, and multiple

solutes in variably saturated media, respectively. Both programs use finite elements to numerically solve the

Richards equation for saturatedunsaturated water flow and Fickian-based advectiondispersion equations

for both heat and solute transport. The flow equation includes a sink term to account for water uptake

by plant roots as a function of both water and salinity stress. The flow equation may also consider

dual-porosity-type flow with a fraction of water content being mobile, and fraction immobile. The heat

transport equation considers conduction as well as advection with flowing water. The solute transport

equations assume advective-dispersive transport in the liquid phase, and diffusion in the gaseous phase.

The transport equations also include provisions for nonlinear and/or nonequilibrium reactions between

the solid and liquid phases, linear equilibrium reactions between the liquid and gaseous phases, zeroorder

production, and two firstorder degradation reactions: one which is independent of other solutes, and

one which provides the coupling between solutes involved in the sequential first-order decay reactions.

In addition, physical nonequilibrium solute transport can be accounted for by assuming a two-region,

dual-porosity type formulation which partitions the liquid phase into mobile and immobile regions.

Alternatively, the transport equations include provisions for kinetic attachment/detachment of solute to

the solid phase and thus can be used to simulate the transport of viruses, colloids, or bacteria (Schijven

and Šim  ̊unek, 2002; Bradford et al., 2003).

Both programs may be used to analyze water and solute movement in unsaturated, partially saturated,

or fully saturated media. The flow region itself may consist of nonuniform (layered) soils. HYDRUS-2D

additionally can handle flow regions delineated by irregular boundaries. The flow region itself may be
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 34 — #34

22-34 The Handbook of Groundwater Engineering

composed of nonuniform soils having an arbitrary degree of local anisotropy. Flow and transport can occur

in the vertical plane, the horizontal plane, or in a three-dimensional region exhibiting radial symmetry

about the vertical axis. The unsaturated soil hydraulic properties (the constitutive relationships) can be

described using van Genuchten (1980a), Brooks and Corey (1964), Kosugi (1996), and Durner (1994)

type analytical functions, or modified van Genuchten type functions that produce a better description of

the hydraulic properties near saturation. HYDRUS-1D incorporates hysteresis by assuming that drying

scanning curves are scaled from the main drying curve, and wetting scanning curves from the main wetting

curve. Root growth is simulated by means of a logistic growth function, while root water uptake can be

simulated as a function of both water and salinity stress.

The HYDRUS-1D software package additionally includes modules for simulating carbon dioxide and

major ion solute movement (Šim  ̊unek et al., 1996). Diffusion in both the liquid and gas phases and advection in the liquid phase are considered as the main CO2 transport mechanisms (Šim  ̊unek and Suarez,

1993). The major variables of the chemical system are Ca, Mg, Na, K, SO4, Cl, NO3, H4SiO4, alkalinity, and CO2(Šim  ̊unek and Suarez, 1994; Šim  ̊unek et al., 1996). The model accounts for equilibrium

chemical reactions between these components such as complexation, cation exchange and precipitationdissolution. For the precipitation-dissolution of calcite and dissolution of dolomite, either equilibrium

or multicomponent kinetic expressions are used, involving both forward and back reactions. Other

dissolution-precipitation reactions considered include gypsum, hydromagnesite, nesquehonite, and sepiolite. Since the ionic strength of soil solutions can vary considerably with time and space and often reach

high values, both modified Debye–Hückel and Pitzer expressions were incorporated into the model as

options to calculate single ion activities.

In addition, both HYDRUS programs implement a Marquardt-Levenberg type parameter estimation

technique (Marquardt, 1963; Šim  ̊unek and Hopmans, 2002) for inverse estimation of soil hydraulic

(Šim  ̊unek et al., 1998b,c; Hopmans et al., 2002) and/or solute transport and reaction (Šim  ̊unek et al.,

2002) parameters from measured transient or steady-state flow and/or transport data.

The HYDRUS packages use a Microsoft Windows based Graphics User Interface (GUI) to manage the

input data required to run the program, as well as for nodal discretization and editing, parameter allocation,

problem execution, and visualization of results. All spatially distributed parameters, such as those for

various soil horizons, the root water uptake distribution, and the initial conditions for water, heat and

solute movement, are specified in a graphical environment. HYDRUS-2D includes the MeshGen2D mesh

generator, specifically designed for variably-saturated subsurface flow transport problems, for defining

more general domain geometries, and for discretizing the transport domain into an unstructured finite

element mesh. The program offers graphs of the distributions of the pressure head, water content, water

and solute fluxes, root water uptake, temperature, and solute concentrations in the soil profile at preselected times. Also included is a small catalog of unsaturated soil hydraulic properties (Carsel and Parish,

1988), as well as pedotransfer functions based on neural networks (Schaap et al., 2001). Many different

applications of both software packages were referenced throughout this book chapter.

22.5.2.1.2 MODFLOW-SURFACT

MODFLOW-SURFACT (HydroGeoLogic, 1996) is a fully integrated flow and transport code, based on

the MODFLOW groundwater modeling software package of the U.S. Geological Survey (McDonald and

Harbaugh, 1988). While MODFLOW (see detailed discussion in Chapter 20) deals only with fully saturated

groundwater flow, MODFLOW-SURFACT expands the applicability of the code to unsaturated domains.

MODFLOW-SURFACT includes new flow packages that enhance schemes for performing unconfined

simulations to rigorously model desaturation and resaturation of aquifers. It provides an option for

discretizing the domain using an axi-symmetric geometry for efficient simulation of pumping tests,

recovery tests, etc. In addition to normally allowed external stresses in MODFLOW (i.e., constant head,

constant flux, areal recharge, evapotranspiration, drains, and streams), MODFLOW-SURFACT provides a

rigorous well withdrawal package, unconfined recharge boundary conditions, and seepage face boundary

conditions. Additionally, MODFLOW-SURFACT also includes options for adaptive time-stepping and
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 35 — #35

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-35

output control procedures, and a new solution package based on the preconditioned conjugate gradient

method.

MODFLOW-SURFACT includes new Analysis of Contaminant Transport (ACT) modules that allow

simulations of single-species and multi-component contaminant transport. The program can handle up to

5 chemical species, which may undergo linear or nonlinear retardation, first-order decay and biochemical

degradation, as well as generate transformation products. MODFLOW-SURFACT also permits mass

partitioning of a single or multicomponent contaminant with diffusive mass movement in the inactive

phase, and mass partitioning of a single or multicomponent contaminant from a depleting immobile

NAPL phase with advective and dispersive transport in the active phase and diffusive transport in the

inactive phase.

22.5.2.2 Geochemical Transport Models

Significant efforts have recently been recently carried out also in coupling physical flow and transport

models with geochemical models to simulate ever more complex reactions, such as surface complexation,

precipitation/dissolution, cation exchange, and biological reactions. Recent reviews of the development

of hydrogeochemical transport models involving reactive multiple components are given by Mangold and

Tsang (1991), Lichtner (1996) and Steefel and MacQuarrie (1996), Šim  ̊unek and Valocchi (2002), and Bell

and Binning (2004). Most modeling efforts involving multicomponent transport have thus far focused on

the saturated zone, where changes in the water flow velocity, temperature and pH are often much more

gradual and hence less important than in the unsaturated zone. Consequently, most multicomponent

transport models assumed one- or two-dimensional steady-state saturated water flow with a fixed value of

the flow velocity, temperature and pH. Only recently have several multicomponent transport models been

published which also consider variably saturated flow; these include DYNAMIX (Liu and Narasimhan,

1989), HYDROGEOCHEM (Yeh and Tripathi, 1990), TOUGH-REACT (Pruess, 1991), UNSATCHEM

(Šim  ̊unek and Suarez, 1994; Šim  ̊unek et al., 1996, 1997), FEHM (Zyvoloski et al., 1997), MULTIFLO

(Lichtner and Seth, 1996), OS3D/GIMRT (Steefel and Yabusaki, 1996), HYDROBIOGEOCHEM (Yeh

et al., 1998), FLOTRAN (Lichtner, 2000), MIN3P (Mayer et al., 2002), HP1 (Jacques et al., 2002; Jacques

and Šim  ̊unek, 2005), and HYDRUS-1D (Šim  ̊unek et al., 2005).

Geochemical models can be divided into two major groups: those with specific chemistry and general

models (Šim  ̊unek and Valocchi, 2002). Models with specific chemistry are limited in the number of species

they can handle, while their application is restricted to problems having a prescribed chemical system.

They are, however, much easier to use and computationally can be much more efficient than general

models. Typical examples of models with specified chemistry are those simulating the transport of major

ions, such as LEACHM (Wagenet and Hutson, 1987), UNSATCHEM (Šim  ̊unek and Suarez, 1994; Šim  ̊unek

et al., 1996), and HYDRUS-1D (Šim  ̊unek et al., 2005). Models with generalized chemistry (DYNAMIX,

HYDROGEOCHEM, MULTIFLO, FLOTRAN, OS3D/GIMRT, and HP1, all referenced above) provide

users with much more freedom in designing a particular chemical system; possible applications of these

models are also much wider.

22.5.2.2.1 HP1

HYDRUS-1D was recently coupled with the PHREEQC geochemical code (Parkhurst and Appelo, 1999)

to create a new comprehensive simulation tool, HP1 (acronym for HYDRUS1D-PHREEQC) (Jacques

et al., 2003; Jacques and Šim  ̊unek, 2005). The combined code contains modules simulating (1) transient

water flow in variably saturated media, (2) the transport of multiple components, (3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heat transport. HP1 is a significant expansion of the

individual HYDRUS-1D and PHREEQC programs by preserving most of their original features and

capabilities. The code still uses the Richards equation for simulating variably saturated water flow and

advection–dispersion type equations for heat and solute transport. However, the program can now simulate also a broad range of low-temperature biogeochemical reactions in water, the vadose zone and in

ground water systems, including interactions with minerals, gases, exchangers, and sorption surfaces,

based on thermodynamic equilibrium, kinetics, or mixed equilibrium-kinetic reactions.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 36 — #36

22-36 The Handbook of Groundwater Engineering

0.0008

0.0006

0.0004

0.0002

0

0.012

0.01

0.008

0.006

0.004

0.002

0

0 3 6 9 12 15 0 3 6 9 12 15

Time [days]

(a) (b)

Time [days]

Na

Ca

Cd

Pb

Zn

Cl

Al

Concentrations [mol/L]

Concentrations [mol/L]

FIGURE 22.13 Major ions (a) and heavy metals (Zn, Pb, and Cd) (b) concentrations in the effluent of an 8-cm soil

long column.

Jacques et al. (2003) and Jacques and Šim  ̊unek (2005) demonstrated the versatility of the HP1 model

on several examples such as (a) the transport of heavy metals (Zn2+, Pb2+, and Cd2+) subject to multiple cation exchange reactions, (b) transport with mineral dissolution of amorphous SiO2 and gibbsite

(Al(OH)3), (c) heavy metal transport in a medium with a pH-dependent cation exchange complex,

(d) infiltration of a hyperalkaline solution in a clay sample (this example considers kinetic precipitation–

dissolution of kaolinite, illite, quartz, calcite, dolomite, gypsum, hydrotalcite, and sepiolite), (e) long-term

transient flow and transport of major cations (Na+, K+, Ca2+, and Mg2+) and heavy metals (Cd2+, Zn2+,

and Pb2+) in a soil profile, (f) cadmium leaching in acid sandy soils, (g) radionuclide transport (U and

its aqueous complexes), and (h) the fate and subsurface transport of explosives (TNT and its daughter

products 2ADNT, 4ADNT, and TAT) (Šim  ̊unek et al., 2006).

22.5.2.2.2 Leaching of Heavy Metals from a Soil Column

As an example application of HP1, Figure 22.13 shows calculated effluent concentrations of heavy metals

(Zn2+, Pb2+, and Cd2+) leached from an 8-cm long soil column having an initial solution as defined in

Table 22.2, and with its ion-exchange complex in equilibrium with this solution. Heavy metals initially

present in the soil column are leached from the exchange complex by major ions (Ca2+, Mg2+, and Al3+).

Water was applied to the top of the column at a steady flow rate of 2 cm day−1 and having a chemical

composition as given in Table 22.2. Details are given by Jacques and Šim  ̊unek (2005).

22.6 Concluding Remarks

This chapter demonstrates the abundance of models and modeling approaches that are currently available

for simulating variably saturated water flow and contaminant transport at various levels of approximation

and for different applications. Models range from relatively simple analytical approaches for analyzing

solute transport problems during one-dimensional steady-state flow, to sophisticated numerical models

for addressing multi-dimensional variably saturated flow and contaminant transport problems at the field

scale.

One may expect that unsaturated zone flow and transport models will be used increasingly as tools for

developing cost-effective, yet technically sound strategies for resource management and pollution remediation and prevention. Improved understanding of the key underlying processes, continued advances in

numerical methods, and the introduction of more and more powerful computers make such simulations increasingly practical for many field-scale problems. For example, models can be helpful tools for

designing, testing and implementing soil, water and crop management practices that minimize soil and

water pollution. Models are equally needed for designing or remediating industrial waste disposal sites
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 37 — #37

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-37

TABLE 22.2 Chemical Components (mmol/L) and Species Considered, and Elemental

Compositions of Initial and Boundary Solutions Used in the Column Simulation. X Refers

to Ion Exchanger

Solutions

Components Species Boundary Initial

Al Al3+, Al(OH)2+, Al(OH)+

2 , Al(OH)0

3, Al(OH)−

4 0.1 0.5

Br Br− 3.7 11.9

Cl Cl− (and Cd, Pb, and Zn-species, see below) 10 0

Ca Ca2+, Ca(OH)+ 5 0

K K+, KOH0 0 2

Na Na+, NaOH0 0 6

Mg Mg2+, Mg(OH)+ 1 0.75

Cd Cd2+, Cd(OH)+, Cd(OH)0

2, Cd(OH)−

3 , Cd(OH)2−4 , CdCl+,

CdCl2, CdCl−

3 0 0.9

Pb Pb2+, Pb(OH)+, Pb(OH)0

2, Pb(OH)−

3 , Pb(OH)2−4 , PbCl+, PbCl0

2,

PbCl−

3 , PbCl2−4 0 0.1

Zn Zn2+, Zn(OH)+, Zn(OH)0

2, Zn(OH)−

3 , Zn(OH)2−4 , ZnCl+, ZnCl0

2,

ZnCl−

3 , ZnCl2−4 0 0.25

X AlX3, AlOHX2, CaX2, CdX2, KX, NaX, MgX2, PbX2, ZnX2 (mmol) NA 11

and landfills, for predicting contaminant transport from mining wastes, or for long-term stewardship of

nuclear waste repositories. The main challenge is to make the models as realistic as possible for the various

applications.

Continued progress in subsurface flow and transport modeling requires equal advances in both numerical techniques as well as the underlying science. Addressing preferential flow phenomena, and the related

problems of subsurface heterogeneity, including the stochastic nature of boundary conditions (precipitation and/or evapotranspiration), will continute to pose formidable challenges. The same is true for

improving multicomponent geochemical transport modeling for the vadose zone. For example, numerical

algorithms and databases for multicomponent transport models must be extended to higher temperatures

and ionic strengths, complex contaminant mixtures (including especially mixed organic and inorganic

wastes), multiphase flow, redox disequilibria for low-temperature systems, and coupled physico-chemical

systems to account for possible changes in the soil water retention and hydraulic conductivity functions.

Better integration is also needed between variably satuarted subsurface and existing larger-scale surface

numerical models, which in turn requires further research on such issues as spatial and temporal scaling

of hydrological, chemical and biological processes and properties, linking constitutive (soil hydraulic)

relationships to measurements scales, preferential flow, and issues of parameter and model uncertainty.

Many science questions related to colloid and colloid-facilitated transport are also still largely unresolved. This is an area of research where our understanding lags far behind current numerical capabilities.

Much work is needed to better understand the processes of filtration, straining, size exclusion, colloid–

colloid interactions, mobilization of colloids and microorganisms; accumulation at air–water interfaces,

interactions between microorganisms and contaminants (including biodegradation), the effects of both

physical factors (water content, flow velocity, textural interfaces) and chemical processes (ionic strength,

solution composition, pH) on colloid retention and mobilization, and modeling colloid-facilitated

transport during conditions of transient flow.

Also, to the best of our knowledge, no models are currentely available that consider all the various

processes simultaneously and in their full complexity, including their mutual interactions. That is, no

models exist that consider transient preferential flow and transport in structured soils or fractured rocks,

while simultaneously considering complex biogeochemical reactions between contaminants, organic and

inorganic colloids and/or organic complexes, and solid and air phases of a soil, including widely varying rates of these various reactions. Further integration of the different types of numerical models is
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 38 — #38

22-38 The Handbook of Groundwater Engineering

needed to address practical problems of contaminant transport (trace elements, radionuclides, organic

contaminants) in complex vadose zone environments.

References

Abriola, L.M., Pennell, K.D., Weber, Jr., W.J. Lang, J.R. and Wilkins, M.D., Persistence and interphase

mass transfer of organic contaminants in the unsaturated zone: Experimental observations and

mathematical modeling. In J.Y. Parlange and J.W. Hopmans (Eds.) Vadose Zone Hydrology: Cutting

Across Disciplines. Oxford University Press, New York, pp. 210–234, 1999.

Adamczyk, Z., Siwek, B., Zembala, M. and Belouschek, P., Kinetics of localized adsorption of colloid

particles. Adv. Coll. Interf. Sci., 48, 151–280, 1994.

Ahuja, L.R. and Hebson, C., Root Zone Water Quality Model. GPSR Tech. Rep. No. 2, USDA, ARS, Fort

Collins, CO, 1992.

Anderson, M.P., Movement of contaminants in groundwater: Groundwater transport–advection and

dispersion. In Groundwater Contamination, Studies in Geophysics, National Academy Press,

Washington, D.C. pp. 37–45, 1984.

Appelo, C.A.J. and Postma, D., Geochemistry, Groundwater and Pollution. A. A. Balkema. Rotterdam, 1993.

Bache, B.W. and Williams, E.G., A phosphate sorption index for soils. J. Soil Sci. 27, 289–296, 1971.

Barry, D.A., Modelling contaminant transport in the subsurface: theory and computer programs. In

H. Ghadiri and C.W. Rose (Eds.), Modelling Chemical Transport in Soil: Natural and Applied

Contaminants, Lewis Publishers, Chelsea, MI, pp. 105–144, 1992.

Bateman, H., The solution of a system of differential equations occurring in the theory of radioactive

transformation. Proc. Cambridge Philos. Soc., 15 (pt V) 423–427, 1910.

Bear, J., Dynamics of Fluids in Porous Media, Elsevier Science, New York, 1972.

Bell, L.S.J. and Binning, P.J., A split operator approach to reactive transport with the forward particle

tracking Eulerian–Lagrangian localized adjoint method. Adv. Water Resour. 27, 323–334, 2004.

Bethke, C. M., Geochemical Reaction Modeling — Concepts and Applications. Oxford University Press,

New York, p. 397, 1996.

Biggar, J.W. and Nielsen, D.R., Miscible displacement and leaching phenomena. Agronomy 11, 254–274,

1967.

Bradford, S.A., Yates, S.R., Bettahar, M. and Šim  ̊unek, J., Physical factors affecting the transport and fate of

colloids in saturated porous media. Water Resour. Res. 38, 1327, doi:10.1029/2002WR001340, 2002.

Bradford, S. A., Šim  ̊unek, J., Bettahar, M., van Genuchten, M. Th. and Yates, S. R., Modeling colloid

attachment, straining, and exclusion in saturated porous media. Environ. Sci. Technol., 37, 2242–

2250, 2003.

Bradford, S. A., Bettahar, M., Šim  ̊unek, J. and van Genuchten, M. Th., Straining and attachment of colloids

in physically heterogeneous porous media. Vadose Zone J. 3, 384–394, 2004.

Bresler, E. and Dagan G., Solute dispersion in unsaturated heterogeneous soil at field scale: I Applications.

Soil Sci. Soc. Am. J. 43, 467–472, 1979.

Brooks, R.H. and Corey, A. T., Hydraulic properties of porous media. Hydrol. Paper No. 3, Colorado State

Univ., Fort Collins, CO, 1964.

Brusseau, M.L., Nonideal transport of reactive solutes in porous media; Cutting across history and

disciplines. In Parlange, M.B. and Hopmans, J.W. (Eds.), Vadose Zone Hydrology: Cutting Across

Disciplines, Oxford University Press, New York, pp. 130–154, 1999.

Carsel, R.F. and Parrish, R. S., Developing joint probability distributions of soil water retention

characteristics. Water Resour. Res. 24, 755–769, 1988.

Casey, F.X.M. and Šim  ̊unek, J., Inverse analyses of the transport of chlorinated hydrocarbons subject to

sequential transformation reactions. J. Environ. Qual., 30, 1354–1360, 2001.

Casey, F.X.M., Larsen, G.L., Hakk, H. and Šim  ̊unek, J., Fate and transport of 17β-Estradiol in soil-water

systems. Environ. Sci. Technol. 37, 2400–2409, 2003.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 39 — #39

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-39

Casey, F.X.M., Larsen, G.L., Hakk, H. and Šim  ̊unek, J., Fate and transport of testosterone in agriculturally

significant soils. Environ. Sci. Technol. 38, 790–798, 2004.

Corapcioglu, M.Y. and Jiang, S., Colloid-facilitated ground-water contaminant transport. Water Resour.

Res. 29, 2215–2226, 1993.

Corapcioglu, M.Y. and Kim, S., Modeling facilitated contaminant transport by mobile bacteria. Water

Resour. Res. 31, 2639–2644, 1995.

AQ: Please update

Cotta et al. Cotta, R.M., Goldstein, R., Heilbron, P.F.L. and Ungs, M.J. (Eds.), Hybrid Methods in Environmental Transport Phenomena, Center for Analysis and Simulation in Environmental Engineering,

POLI/COPPE, Federal University of Rio de Janeiro, Brazil, 2005 (in press).

Dagan G. and E. Bresler, Solute dispersion in unsaturated heterogeneous soil at field scale: I Theory. Soil

Sci. Soc. Am. J. 43, 461–467, 1979.

Dagan G., Flow and Transport in Porous Formations, Springer-Verlag, Berlin, 1989.

Durner, W., Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour.

Res. 32, 211–223, 1994.

Fayer, M.J., UNSAT-H Version 3.0: Unsaturated Soil Water and Heat Flow Model. Theory, User Manual,

and Examples. Pacific Northwest National Laboratory 13249, 2000.

Fitter, A.H. and Sutton, S.D., The use of the Freundlich isotherm for soil phosphate sorption data. J. Soil

Sci. 26, 241–246, 1975.

Flerchinger, G.N., Hanson, C.L. and Wight, J.R., Modeling evapotranspiration and surface energy budgets

across a watershed. Water Resour. Res. 32, 2539–2548, 1996.

Flury, M., Fluhler, H., Jury, W. and Leuenberger, J., Susceptibility of soils to preferential flow of water:

a field study. Water Resour. Res. 30, 1945–1954, 1994.

Freundlich, H., Kapillarchemie; eine Darstellung der Chemie der Kolloide und verwandter Gebiete,

Akademische Verlagsgellschaft, Leipzig, Germany, 1909.

Gelhar, L.W., Mantaglou, A., Welty, C. and Rehfeldt, K.R., A review of Field Scale Physical Solute Transport

Processes in Unsaturated and Saturated Porous Media, EPRI Topicl Report EA-4190, Electric Power

Research Institute, Palo Alto, CA, 1985.

Gerke, H.H. and van Genuchten, M.Th., A dual-porosity model for simulating the preferential movement

of water and solutes in structured porous media. Water Resour. Res. 29, 305–319, 1993a.

Gerke, H.H. and van Genuchten, M.Th., Evaluation of a first-order water transfer term for variably

saturated dual-porosity flow models. Water Resour. Res. 29, 1225–1238, 1993b.

Gerke, H.H. and van Genuchten, M.Th., Macroscopic representation of structural geometry for simulating

water and solute movement in dual-porosity media. Adv. Water Resour. 19, 343–357, 1996.

Germann, P.F., Kinematic wave approach to infiltration and drainage into and from soil macropores.

Trans. ASAE 28, 745–749, 1985.

Germann, P.F. and Beven, K., Kinematic wave approximation to infiltration into soils with sorbing

macropores. Water Resour. Res. 21, 990–996, 1985.

Glotfelty, D.E. and Schomburg, C.J., Volatilization of pesticides. In B.L. Sawhney and K. Brown (Eds.),

Reactions and Movement of Organic Chemicals in Soils, Special Publ. No. 22, Soil Science Society of

America, Madison, WI, pp. 181–208, 1989.

Grolimund, D., Borkovec, M., Barmettler, K. and Sticher, H., Colloid-facilitated transport of strongly

sorbing contaminants in natural porous media; a laboratory column study. Environ. Sci. Technol.,

30, 3118–3123, 1996.

Gunary, D., A new adsorption isotherm for phosphate in soil. J. Soil Sci. 21, 72–77, 1970.

Gwo, J.P., Jardine, P.M., Wilson, G.V. and Yeh, G.T., A multiple-pore-region concept to modeling mass

transfer in subsurface media. J. Hydrol. 164, 217–237, 1995.

Hanselman, T.A., Graetz, D.A. and Wilkie, A.C., Manure-borne estrogens as potential environmental

contaminants: A review. Environ. Sci. Technol. 37, 5471–5478, 2003.

Hansen, S., Jensen, H.E., Nielsen, N.E. and Svendsen, H., DAISY: Soil Plant Atmosphere System

Model. NPO Report No. A 10. The National Agency for Environmental Protection,Copenhagen,

272 pp, 1990.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 40 — #40

22-40 The Handbook of Groundwater Engineering

Haws, N.W., Rao, P.S.C. and Šim  ̊unek, J., Single-porosity and dual-porosity modeling of water flow and

solute transport in subsurface-drained fields using effective field-scale parameters. J. Hydrol. 2005

(in press).

AQ: Please update

Haws et al.

Healy, R.W., Simulation of solute transport in variably saturated porous media with supplemental information on modifications to the U.S. Geological Survey’s computer program VS2D. U.S. Geological

Survey, Water-Resources Investigation Report 90-4025: 125, 1990.

Hendrickx, J.M.H., Dekker, L.W. and Boersma, O.H., Unstable wetting fronts in water repellent field soils.

J. Env. Qual. 22, 109–118, 1993.

Hendrickx, J.M.H. and Flury, M., Uniform and preferential flow, Mechanisms in the vadose zone. In

Conceptual Models of Flow and Transport in the Fractured Vadose zone. National Research Council,

National Academy Press, Washington, D.C., pp. 149–187, 2001.

Hopmans, J.W., Šim  ̊unek, J., Romano, N. and Durner, W., Inverse modeling of transient water flow. In

J.H. Dane and G.C. Topp (Eds.), Methods of Soil Analysis, Part 1, Physical Methods, Chapter 3.6.2,

Third edition, SSSA, Madison, WI, pp. 963–1008, 2002.

Hutson, J.L. and Wagenet, R.J., A multiregion model describing water flow and solute transport in

heterogeneous soils. Soil Sci. Soc. Am. J. 59, 743–751, 1995.

Huyakorn, P.S. and Pinder, G.F., Computational Methods in Subsurface Flow, Academic Press, London,

1983.

HydroGeoLogic, Inc., MODFLOW-SURFACT, A Comprehensive MODFLOW-Based Flow and Transport

Simulator, Version 2.1, HydroGeoLogic, Inc., 1996.

Jacques, D., Šim  ̊unek, J., Mallants, D. and van Genuchten, M. Th., Multicomponent transport model

for variably-saturated porous media: Application to the transport of heavy metals in soils. In

S.M., Hassanizadeh, R.J. Schotting, W.G. Gray, and G.F. Pinder (Eds.) Computational Methods in

Water Resources, XIVth International Conference on Computational Methods in Water Resources,

June 23–28, 2002, Delft, The Netherlands. Elsevier, Developments in Water Science, 47, 555–562,

2002.

Jacques, D., Šim  ̊unek, J., Mallants, D. and van Genuchten, M. Th., The HYDRUS-PHREEQC multicomponent transport model for variably-saturated porous media: Code verification and application,

MODFLOW and More 2003: Understanding through Modeling, Conference Proceedings, E. Poeter,

Ch. Zheng, M. Hill, and J. Doherty, (Eds.), Int. Ground Water Modeling Center, Colorado School

of Mines, pp. 23–27, 2003.

Jacques, D. and Šim  ̊unek, J., Multicomponent — Variable Saturated Transport Model, Description,

Manual, Verification and Examples, Waste and Disposal, SCKCEN, BLG-998, Mol, Belgium, 79 pp,

2005.

Jansson, P.-E. and Karlberg, L., Coupled heat and mass transfer model for soil-plant-atmosphere systems.

Royal Institute of Technology, Dept of Civil and Environmental Engineering, Stockholm 325 pp,

2001.

Jarvis, N.J., The MACRO model (Version 3.1), Technical description and sample simulations. Reports and

Dissertations 19. Dept. Soil Sci., Swedish Univ. Agric. Sci., Uppsala, Sweden, 51 pp, 1994.

Jarvis, N. J., Modeling the impact of preferential flow on nonpoint source pollution. In H.M. Selim and

L. Ma (Eds.), Physical Nonequilibrium in Soils: Modeling and Application. Ann Arbor Press, Chelsea,

MI, pp. 195–221, 1998.

Jiang, S. and Corapcioglu, M.Y., A hybrid equilibrium model of solute transport in porous media in the

presence of colloids. Coll. Surf. A: Physicochem. Eng. Aspects 73, 275–286, 1993.

Jury, W.A., Simulation of solute transport using a transfer function model. Water Resour. Res. 18, 363–368,

1982.

Jury, W.A., Spencer, W.F. and Farmer, W.J., Behavior assessment model for trace organics in soil: I.

Description of model. J. Env. Qual. 12, 558–564, 1983.

Jury, W.A., Spencer W.F. and Farmer, W.J., Behavior assessment model for trace organics in soil: I. Chemical

classification and parameter sensitivity. J. Env. Qual. 13, 567–572, 1984.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 41 — #41

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-41

Jury W.A. and Sposito, G., Field calibration and validation of solute transport models for the unsaturated

zone. Soil Sci. Soc. Am. J. 49, 1331–1341, 1985.

Jury, W.A. and Horton, R., Soil Physics, Sixth ed, John Wiley & Sons, Inc., New York, 2004.

Kan, A.T. and Tomson, M.B., Ground water transport of hydrophobic organic compounds in the presence

of dissolved organic matter. Env. Toxicol. Chem. 9, 253–263, 1990.

Kodešová, R., Kozák, J., Šim  ̊unek, J. and Vacek, O., Field and numerical study of chlorotoluron transport

in the soil profile: Comparison of single and dual-permeability model. Plant, Soil Environ. 51,

310–315, 2005.

Köhne, J.M., Köhne, S., Mohanty, B.P. and Šim  ̊unek, J., Inverse mobile–immobile modeling of transport

during transient flow: Effect of between-domain transfer and initial water content. Vadose Zone J.

3, 1309–1321, 2004a.

Köhne, J.M., B. Mohanty, J. Šim  ̊unek and H.H. Gerke, Numerical evaluation of a second-order

water transfer term for variably saturated dual-permeability models,Water Resour. Res. 40,

doi:10.1029/2004WR00385, 2004b.

Köhne, J.M. and Mohanty, B.P., Water flow processes in a soil column with a cylindrical macropore: Experiment and hierarchical modeling. Water Resour. Res., 41, W03010, doi:10,1029/2004WR003303,

2005.

Köhne, S., Lennartz, B., Köhne, J.M. and Šim  ̊unek, J., Bromide transport at a tile-drained field site:

Experiment, one- and two-dimensional equilibrium and non-equilibrium numerical modeling.

J. Hydrol. 2005 (in press).

AQ: Please update

Köhne et al. Köhne, J.M., Köhne, S. and Šim  ̊unek, J., Multi-process herbicide transport in structured soil columns:

Experiment and model analysis. J. Contamin. Hydrol. (submitted).

Kosugi, K., Lognormal distribution model for unsaturated soil hydraulic properties, Water Resour. Res.,

32, 2697–2703, 1996.

Kung, K.-J.S., Preferential flow in a sandy vadose zone. 2. Mechanism and implications. Geoderma 46,

59–71, 1990.

Lai, Sung-Ho and Jurinak, J.J., Numerical approximation of cation exchange in miscible displacement

through soil columns. Soil Sci. Soc. Am. Proc. 35, 894–899, 1971.

Langmuir, I., The adsorption of gases on plane surfaces of glass, mica, and platinum. J. Am. Chem. Soc.

40, 1361–1403, 1918.

Lapidus, L. and Amundson, N.R., Mathematics of adsorption in beds: VI. The effect of longitudinal

diffusion in ion exchange and chromatographic columns. J. Phys. Chem. 56, 984–988, 1952.

Leij, F.J. and Bradford, S.A., 3DADE: A computer program for evaluating three-dimensional equilibrium

solute transport in porous media. Research Report No. 134, U.S. Salinity Laboratory, USDA, ARS,

Riverside, CA, 1994.

Leij, F.J. and Toride, N., N3DADE: A computer program for evaluating nonequilibrium three-dimensional

equilibrium solute transport in porous media. Research Report No. 143, U.S. Salinity Laboratory,

USDA, ARS, Riverside, CA, 1997.

Leij, F.J., Priesack, E. and Schaap, M.G., Solute transport modeled with Green’s functions with application

to persistent solute sources. J. Contam. Hydrol., 41, 155–173, 2000.

Lenhard, R.J., Oostrom, M. and Dane, J.H., Multi-fluid flow. In J.H. Dane and G.C. Topp (Eds.), Methods

of Soil Analysis, Part 4, Physical Methods, SSSA Book Series 5, Soil Sci. Soc, Am., Madison, WI,

pp. 1537–1619, 2002.

Lichtner, P.C., Continuum formulation of multicomponent-multiphase reactive transport. In P.C. Lichtner, C.I. Steefel and E.H. Oelkers (Eds.), Reactive Transport in Porous Media, Reviews in Mineralogy.

Chapter 1, Vol. 34. Mineralogical Society of America, pp. 1–81, 1996.

Lichtner, P.C. and Seth, M.S., User’s manual for MULTIFLO: Part II MULTIFLOW 1.0 and GEM 1.0.

Multicomponent-multiphase reactive transport model. CNWRA 96-010, Center for Nuclear Waste

Regulatory Analyses, San Antonio, TX, 1996.

Lichtner, P.C., FLOWTRAN User’s manual. Los Alamos National Laboratory Document. NM, 2000.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 42 — #42

22-42 The Handbook of Groundwater Engineering

Lindstrom, F.T., Haque, R., Freed, V.H. and Boersma, L., Theory on the movement of some herbicides in

soils. Environ. Sci. Technol. 1, 561–565, 1967.

Lindstrom, F.T., Boersma, L. and Stockard, D., A theory on the mass transport of previously distributed

chemicals in a water saturated sorbing porous medium. Isothermal cases, Soil Sci. 112, 291–300,

1971.

Lindqvist, R. and Enfield, C.G., Biosorption of dichlorodiphenyltri-chloroethane and hexachlorobenzene in groundwater and its implications for facilitated transport. Appl. Environ. Microbiol. 58,

2211–2218, 1992.

Liu, C.W. and Narasimhan, T.N., Redox-controlled multiple-species reactive chemical transport. 1. Model

development. Water Resour. Res. 25:869–882, 1989.

Liu, C., Szecsody, J.E., Zachara, J.M. and Ball, W.P., Use of the generalized integral transform method

for solving equations of solute transport in porous media. Adv. Water. Resour. 23, 483–492,

2000.

Logan, B.E., Jewett, D.G., Arnold, R.G., Bouwer, E.J. and O’Melia, C.R., Clarification of clean-bed filtration

models. J. Environ. Eng. 121, 869–873, 1995.

Magee, B.R., Lion, L.W. and Lemley, A.T., Transport of dissolved organic macromolecules and their effect

on the transport of phenanthrene in porous media. Environ. Sci. Technol. 25, 323–331, 1991.

Mangold, D.C. and Tsang, A summary of subsurface hydrological and hydrochemical models. Rev.

Geophys. 29, 51–79, 1991.

Mansfeldt, T., Leyer, H., Barmettler, K. and Kretzschmar, R., Cyanide leaching from soil developed from

coking plant purifier waste as influenced by citrate. Vadose Zone J. 3, 471–479, 2004.

Marquardt, D.W., An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math.

11, 431–441, 1963.

Mayer, U., Frind, E.O. and Blowes, D.W., Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water

Resour. Res. 38, 1174, 2002.

McDonald, M.G. and Harbaugh, A.W., A modular three-dimensional finite-difference ground-water flow

model, U.S. Geological Survey Techniques of Water-Resources Investigations Book 6, Chapter A1,

586 p, 1988.

MDH Engineered Solutions Corp., Evaluation of Computer Models for Predicting the Fate and Transport

of Salt in Soil and Groundwater, Phase I Report. Pub. No: T/403, ISBN No. 0-7785-2493-0 (Printed

Edition), ISBN No. 0-7785-2494-9 (On-line Edition), 2003.

Millington, R.J. and J.M. Quirk, Permeability of porous solids. Trans. Faraday Soc. 57, 1200–1207, 1961.

Mills W.M., Liu S. and Fong, F.K., Literature review and model (COMET) for colloid/metals transport in

porous media. Groundwater 29, 199–208, 1991.

Moldrup, O., Poulsen, T.G., Schjonning, P., Olesen, T. and Yamaguchi, T., Gas permeability in undisturbed

soils: measurements and predictive models. Soil Sci. 163, 180–189, 1998.

Narasimhan, T.N., Buckingham, 1907; An appreciation. Vadose Zone J. 4, 434–441, 2005,

Nielsen, D.R., van Genuchten, M. Th. and Biggar, J.W., Water flow and solute transport processes in the

unsaturated zone. Water Resour. Res. 22, 89S–108S, 1986.

Noell, A.L., Thompson, J.L., Corapcioglu, M.Y. and Triay, I.R., The role of silica colloids on facilitated

cesium transport through glass bead columns and modeling. J. Contam. Hydrol. 31, 23–56, 1998.

Parker, J.C. and van Genuchten, M.Th., Determining transport parameters from laboratory and field

tracer experiments. Bull. 84-3, Va Agric. Exp. St., Blacksburg, VA, 1984.

Parkhurst, D.L. and Appelo, C.A.J., User’s guide to PHREEQC (version 2) — A computer program for

speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Water

Resources Investigation, Report 99-4259, Denver, Co, USA, 312 pp, 1999.

Perrochet, P. and Berod, D., Stability of the standard Crank–Nicolson–Galerkin scheme applied to the

diffusion-convection equation: some new insights. Water Resour. Res. 29, 3291–3297, 1993.

Pinder, G.F. and Gray, W.G., Finite Element Simulation in Surface and Subsurface Hydrology, Academic

Press, New York, 1977.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 43 — #43

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-43

Pitzer, K.S., Activity coefficients in electrolyte solutions. In R.M. Pitkowicz, Activity Coefficients in

Electrolyte Solutions. CRC Press, Inc., Boca Raton, FL, pp. 157–208, 1979.

AQ: Please update

Pot, V., et al. Pot, V., Šim  ̊unek, J., Benoit, P., Coquet, Y., Yra, A. and Martínez-Cordón, M.-J., Impact of rainfall intensity

on the transport of two herbicides in undisturbed grassed filter strip soil cores. J. Contam. Hydrol.

2005 (in press).

Pruess, K., Tough2 — A general-purpose numerical simulator for multiphase fluid and heat flow. Rep.

LBL-29400, Lawrence Berkeley Laboratory, Berkeley, CA, 1991.

Pruess, K., The TOUGH codes — A family of simulation tools for multiphase flow and transport processes

in permeable media. Vadose Zone J. 3, 738–746, 2004.

Pruess, K. and Wang, J.S.Y., Numerical modeling of isothermal and non-isothermal flow in unsaturated

fractured rock — a review. In D.D. Evans and T.J. Nicholson (Eds.), Flow and Transport Through

Unsaturated Fractured Rock, Geophysics Monograph 42, American Geophysical Union, Washington,

D.C., pp. 11–22, 1987.

Rajagopalan, R. and Tien, C., Trajectory analysis of deep-bed filtration with the sphere-in-cell porous

media model. AIChE. J. 22, 523–533, 1976.

Rathfelder, K.M., Lang, J.R. and Abriola, L.M., A numerical model (MISER) for the simulation of coupled

physical, chemical and biological processes in soil vapor extraction and bioventing systems. J.

Contam. Hydrol. 43, 239–270, 2000.

Richards, L.A., Capillary conduction of fluid through porous mediums. Physics 1, 318–333, 1931.

Saiers, J.E., Hornberger, G.M. and Harvey, C., Colloidal silicakaolinite in the transport of cesium through

laboratory sand columns. Water Resour. Res. 32, 33–41, 1996.

Saiers, J.E. and Hornberger, G.M., The role of colloidal kaolinite in the transport of cesium through

laboratory sand columns. Water Resour. Res. 32, 33–41, 1996.

Scanlon, B.R., Christman, M., Reedy, R., Porro, I., Šim  ̊unek, J. and Flerchinger, G., Intercode comparisons

for simulating water balance of surficial sediments in semiarid regions. Water Resour. Res. 38, 1323,

doi:10.1029/2001WR001233, 59.1–59.16, 2002.

Schaerlaekens, J., Mallants, D., Šim  ̊unek, J. van Genuchten, M.Th. and Feyen, J. Numerical simulation of

transport and sequential biodegradation of chlorinated aliphatic hydrocarbons using CHAIN_2D.

J. Hydrologi. Proc. 13, 2847–2859, 1999.

Schaap, M.G., Leij, F.J. and van Genuchten, M.Th., Rosetta: A computer program for estimating soil

hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251, 163–176, 2001.

Schijven, J.F. and Hassanizadeh, S.M., Removal of viruses by soil passage: Overview of modeling, processes

and parameters. Crit. Rev. Environ. Sci.Technol. 30, 49–127, 2000.

Schijven, J. and Šim  ̊unek, J., Kinetic modeling of virus transport at field scale. J. Contam. Hydrol. 55,

113–135, 2002.

Selim, H.M., Davidson, J.M. and Mansell, R.S., Evaluation of a two-site adsorption-desorption model for

describing solute transport in soils, paper presented at Proceedings, Summer Computer Simulation

Conference, Nat. Sci. Found., Washington, DC, 1976.

Seuntjens, P., Tirez, K., Šim  ̊unek, J., van Genuchten, M.Th., Cornelis, C. and Geuzens, P., Aging effects

on cadmium transport in undisturbed contaminated sandy soil columns. J. Environ. Qual. 30,

1040–1050, 2001.

Shapiro, R.E. and Fried M., Relative release and retentiveness of soil phosphates. Soil Sci. Soc. Am. Proc.

23, 195–198, 1959.

Sibbesen, E., Some new equations to describe phosphate sorption by soils. J. Soil Sci. 32, 67–74, 1981.

Šim  ̊unek, J. and Suarez, D.L., Modeling of carbon dioxide transport and production in soil: 1. Model

development. Water Resour. Res. 29, 487–497, 1993.

Šim  ̊unek, J. and Suarez, D.L., Two-dimensional transport model for variably saturated porous media with

major ion chemistry. Water Resour. Res. 30, 1115–1133, 1994.

Šim  ̊unek, J. and van Genuchten, M.Th., The CHAIN_2D code for simulating two-dimensional movement

of water flow, heat, and multiple solutes in variably-saturated porous media, Version 1.1, Research

Report No 136, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 205 pp, 1994.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 44 — #44

22-44 The Handbook of Groundwater Engineering

Šim  ̊unek, J., Suarez, D.L. and Šejna, M., The UNSATCHEM software package for simulating onedimensional variably saturated water flow, heat transport, carbon dioxide production and transport,

and multicomponent solute transport with major ion equilibrium and kinetic chemistry. Version

2.0. Research Report No. 141, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 186 pp,

1996.

Šim  ̊unek, J. and Suarez, D.L., Sodic soil reclamation using multicomponent transport modeling. ASCE J.

Irrig. Drain. Eng. 123, 367–376, 1997.

Šim  ̊unek, J., Šejna, M. and van Genuchten, M.Th. The HYDRUS-1D software package for simulating

the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media.

Version 2.0, IGWMC-TPS-70, International Ground Water Modeling Center, Colorado School of

Mines, Golden, Colorado, 202pp, 1998a.

Šim  ̊unek, J., Wendroth, O. and van Genuchten, M.Th., A parameter estimation analysis of the evaporation

method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62, 894–905, 1998b.

Šim  ̊unek, J., van Genuchten, M.Th., Gribb, M.M. and Hopmans, J.W., Parameter estimation of unsaturated

soil hydraulic properties from transient flow processes. Soil Tillage Res. 47, Special issue "State of

the art in soil physics and in soil technology of antrophic soils", 27–36, 1998c.

Šim  ̊unek, J., van Genuchten, M.Th., Šejna, M., Toride, N. and Leij, F.J., The STANMOD computer software

for evaluating solute transport in porous media using analytical solutions of convection-dispersion

equation. Versions 1.0 and 2.0, IGWMC-TPS-71, International Ground Water Modeling Center,

Colorado School of Mines, Golden, Colorado, 32 pp, 1999a.

Šim  ̊unek, J., Šejna, M. and van Genuchten, M.Th., The HYDRUS-2D software package for simulating twodimensional movement of water, heat, and multiple solutes in variably saturated media. Version

2.0, IGWMC-TPS-53, International Ground Water Modeling Center, Colorado School of Mines,

Golden, Colorado, 251pp, 1999b.

Šim  ̊unek, J., Wendroth, O., Wypler, N. and van Genuchten, M.Th., Nonequilibrium water flow

characterized from an upward infiltration experiment. Eur. J. Soil Sci. 52, 13–24, 2001.

Šim  ̊unek, J. and Hopmans, J.W., Parameter optimization and nonlinear fitting. In: J.H. Dane and

G.C. Topp, (Eds.) Methods of Soil Analysis, Part 1, Physical Methods, Chapter 1.7, 3rd ed., SSSA,

Madison, WI, pp. 139–157, 2002.

Šim  ̊unek, J., Jacques, D., Hopmans, J.W., Inoue, M., Flury, M. and van Genuchten, M.Th., Solute transport

during variably-saturated flow — inverse methods. In J.H. Dane and G.C. Topp, (Eds.), Methods

of Soil Analysis, Part 1, Physical Methods, Chapter 6.6, 3rd ed., SSSA, Madison, WI, pp. 1435–1449,

2002.

Šim  ̊unek, J. and Valocchi, A.J., Geochemical transport, In J. H. Dane and G.C. Topp, (Eds.), Methods of

Soil Analysis, Part 1, Physical Methods, Chapter 6.9, 3rd ed., SSSA, Madison, WI, pp. 1511–1536,

2002.

Šim  ̊unek, J., Jarvis, N.J., van Genuchten, M.Th. and Gärdenäs, A., Nonequilibrium and preferential flow

and transport in the vadose zone: review and case study. J. Hydrol. 272, 14–35, 2003.

Šim  ̊unek, J., van Genuchten, M.Th. and Šejna, M., The HYDRUS-1D software package for simulating

the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media.

Version 3.0, HYDRUS Software Series 1, Department of Environmental Sciences, University of

California Riverside, Riverside, CA, 270pp, 2005.

AQ: Please update

Šim  ̊unek

Šim  ̊unek, J., Models of water flow and solute transport in the unsaturated zone. In: M.G. Anderson (Ed.),

The Encyclopedia of Hydrological Sciences, John Wiley & Sons, Ltd., Chichester, England, 25 pp,

2005 (in press).

AQ: Please update

Šim  ̊unek et al.

Šim  ̊unek, J., Jacques, D., van Genuchten, M.Th. and Mallants, D., Multicomponent geochemical transport

modeling using the HYDRUS computer software packages. J. Am. Water Resour. Assoc. 2006 (in

press).

Sips, R., On the structure of the catalyst surface, II. J. Chem. Phys. 18, 1024–1026, 1950.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 45 — #45

Contaminant Transport in the Unsaturated Zone Theory and Modeling 22-45

Skaggs, T.H. and Leij, F.J., Solute transport: Theoretical background. In J.H. Dane and G.C. Topp,

(Eds.), Methods of Soil Analysis, Part 1, Physical Methods, Chapter 6.3, 3rd ed., SSSA, Madison, WI,

pp. 1353–1380, 2002.

Sposito, G. and Barry, D.A., On the Dagan model of solute transport I groundwater: Foundational aspects.

Water Resour. Res. 23, 1867–1875, 1987.

Steefel, C.I. and Yabusaki, S.B., OS3D/GIMRT, software for multicomponent–multidimensional reactive transport. User Manual and Programmer’s Guide, PNL-11166, Pacific Northwest National

Laboratory, Richland, WA, 1996.

Steefel, C.I. and MacQuarrie, K.T.B., Approaches to modeling of reactive transport in porous media. In

P.C. Lichtner, C.I. Steefel and E.H. Oelkers (Eds.), Reactive Transport in Porous Media. Reviews in

Mineralogy. Chapter 2, Vol. 34, Mineralogical Society of America, pp. 83–129, 1996.

Thiele-Bruhn, S., Pharmaceutical antibiotic compounds in soils — a review. J. Plant Nutr. Soil Sci. 166,

145–167, 2003.

Tolls, J., Sorption of veterinary pharmaceuticals in soils: A review. Environ. Sci. Technol. 35, 3397–3406.

2001.

Toride, N., Leij, F.J. and van Genuchten, M.Th., A comprehensive set of analytical solutions for nonequilibrium solute transport with first-order decay and zero-order production. Water Resour. Res. 29,

2167–2182, 1993.

Toride, N., Leij, F.J. and van Genuchten, M.Th., The CXTFIT code for estimating transport parameters

from laboratory or field tracer experiments. Version 2.0, Research Report No. 137, U. S. Salinity

Laboratory, USDA, ARS, Riverside, CA, 1995.

Truesdell, A.H. and Jones, B.F., Wateq, a computer program for calculating chemical equilibria of natural

waters. J. Res. U.S. Geol. Surv. 2, 233–248, 1974.

Ungs, M., Wilkinson, K. and Summers, K., Users guide for MYGRT, Version 3.0: Software for simulating

migration of organic and inorganic chemicals in grioundwater, EPRI Report TR-111748, Electric

Power Research Institute, Palo Alto, CA, 1998.

van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabat, P. van Valsum, P.E.V. and Groenendijk,

P. and van Diepen, C.A., Theory of SWAP, version 2.0. Simulation of water flow, solute transport

and plant growth in the Soil–Water–Atmosphere–Plant environment. Dept. Water Resources, WAU,

Report 71, DLO Winand Staring Centre, Wageningen, Technical Document 45, 1997.

Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiery, D.,Vanclooster, M., Mouvet, C. andVereecken,

H., A set of analytical benchmarks to test numerical models of flow and transport in soils. Vadose

Zone J. 4, 206–221, 2005.

van de Weerd, H., Leijnse, A. and van Riemsdijk, W.H., Transport of reactive colloids and contaminants

in groundwater: effect of nonlinear kinetic reactions. J. Contam. Hydrol. 32, 313–331, 1998.

van Genuchten, M.Th. and Wierenga, P.J., Mass transfer studies in sorbing porous media, I. Analytical

solutions. Soil Sci. Soc. Am. J. 40, 473–481, 1976.

van Genuchten, M.Th., A closed-form equation for predicting the hydraulic conductivity of unsaturated

soils. Soil Sci. Soc. Am. J. 44, 892–898, 1980a.

van Genuchten, M.Th. and J. Šim  ̊unek, Evaluation of pollutant transport in the unsaturated zone, In:

P.E. Rijtema and V. Eliáš (Eds.), Regional Approaches to Water Pollution in the Environment, NATO

ASI Series: 2. Environment, Kluwer, Dordrecht, The Netherlands, pp. 139–172, 1996.

van Genuchten, M.Th., Davidson, J.M. and Wierenga, P.J., An evaluation of kinetic and equilibrium

equations for the prediction of pesticide movement through porous media. Soil Sci. Soc. Am. Proc.

38, 29–35, 1974.

van Genuchten, M.Th., Determining transport parameters from solute displacement experiments.

Research Report No. 118, U. S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1980b.

van Genuchten, M.Th., Non-equilibrium transport parameters from miscible displacement experiments.

Research Report No. 119, U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1981.
JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 46 — #46

22-46 The Handbook of Groundwater Engineering

van Genuchten, M.Th. and Wagenet, R.J., Two-site/two-region models for pesticide transport and degradation: Theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53, 1303–1310,

1989.

van Genuchten, M.Th. and Parker, J.C., Boundary conditions for displacement experiments through short

laboratory soil columns. Soil Sci. Soc. Am. J. 48, 703–708, 1984.

van Genuchten, M.Th., Convective–dispersive transport of solutes involved in sequential first-order decay

reactions. Comp. & Geosci. 11, 129–147, 1985.

van Genuchten, M.Th. and Šim  ̊unek, J., Integrated modeling of vadose zone flow and transport processes.

Proc. Unsaturated Zone Modelling: Progress, Challenges and Applications, Eds. R. A. Feddes, G. H. de

Rooij, and J. C. van Dam, Wageningen, The Netherlands, October 3–5, 2004, pp. 37–69, 2004.

Vinten, A.J.A., Yaron, B. and Nye, P.H., Vertical transport of pesticides into soil when adsorbed on

suspended particles. J. Agric. Food Chem. 31, 662–664, 1983.

Von Gunten, H.R., Weber, V.E. and Krahenbuhl, U., The reactor accident at Chernobyl: A possibility to test

colloid-controlled transport of radionuclides in a shallow aquifer. J. Contam. Hydrol. 2, 237–247,

1998.

Wagenet, R.J. and Hutson, J.L., LEACHM: Leaching Estimation And CHemistry Model, A process-based

model of water and solute movement, transformations, plant uptake and chemical reactions in the

unsaturated zone. Continuum 2. Water Resour. Inst., Cornell University, Ithaca, New York, 1987.

Westall, J.C., Zachary, J.L. and Morel, F.M.M. MINEQL: A computer program for the calculation of

chemical equilibrium composition of aqueous systems. Tech. Note 18, 91pp, Dep. of Civ. Eng.,

Mass. Inst. of Technol., Cambridge, MA, 1976.

White, N. and Zelazny, L.W., Charge properties in soil colloids. In D.L. Sparks (Ed.) Soil Physical Chemistry.

CRC Press, Boca Raten, Florida, pp. 39–81, 1986.

Wheatcraft, S.W. and Cushman, J.H., Hierarchical approaches to transport in heterogeneous porous

media. Rev. Geophys. 29 (Suppl.), 261–267, 1991.

Wilson, G.V., Albright, W.H., Gee, G.W., Fayer, M.J. and Ogan, B.D., Alternative cover assessment project,

Phase 1 Report. Final Technical Report prepared by Desert Research Institute, NV and Pacific

Northwest National Laboratory, WA for U.S. Environmental Protection Agency under Contract No.

44-0000-2038, 202p, 1999.

Yeh, G.T., Analytical transient one-, two-, and three-dimensional simulation of waste transport in the

aquifer system, ORNL-5602, Oak Ridge National Laboratory, Oak Ridge, TN, 1981.

Yeh, G.T. and Tripathi, V.S., HYDROGEOCHEM: A coupled model of hydrological and geochemical

equilibrium of multicomponent systems. Rep. ORNL-6371, Oak Ridge Natl. Lab., Oak Ridge, TN,

1990.

Yeh, G.T., Salvage, K.M., Gwo, J.P., Zachara, J.M. and Szecsody, J.E., HYDROBIOGEOCHEM: A coupled

model of hydrological transport and mixed biochemical kinetic/equilibrium reactions in saturated–

unsaturated media. Rep. ORNL/TM-13668, Oak Ridge Natl. Lab., Oak Ridge, TN, 1998.

Zhang, P., Šim  ̊unek, J. and Bowman, R.S., Nonideal transport of solute and colloidal tracers through

reactive zeolite/iron pellets. Water Resour. Res. 40, doi:10.1029/2003WR002445, 2004.

Zienkiewicz, O.C., The Finite Element Method, 3rd ed., McGraw-Hill, London, 1977.

Zyvoloski, G.A., Robinson, B.A., Dash, Z.V. and Trease, L.L., Summary of the models and methods

for the FEHM application — a finite element heat- and mass-transfer code. Los Alamos National

Laboratory Rept. LA-13307-MS, NM, 1997.