Chapter 18

Leaching of Contaminants to Groundwater
Dirk Mallants, Martinus Th. Van Genuchten, Jiří Šimůnek, Diederik Jacques,
and Suresh Seetharam

Abstract In this chapter the water flow and contaminant transport processes in theunsaturated or vadose zone are described. These processes include water retention and hydraulic conductivity, evapotranspiration, preferential flow, root water
uptake (water flow) and diffusion, dispersion, advection and volatilization (contaminant transport). The equation governing transport of dissolved contaminants in thevadose zone is obtained by combining the contaminant mass balance with equations
defining the total concentration of the contaminant and the contaminant flux density. Further attention is this chapter is given to nonequilibrium transport, stochastic
models, multicomponent reactive solute transport, multiphase flow and transport.
Mathematical models should be critical components of any effort to understandand predict site-specific subsurface water flow and contaminant transport processes.
Generally, models range from relatively simple analytical approaches for analyzing contaminant transport problems during one-dimensional steady-state flow, tosophisticated numerical models for addressing multi-dimensional variably-saturated
flow and contaminant transport problems at the field scale. An overview is given of
several existing analytical and numerical models. Moreover, several applications to
unsaturated flow and geochemical transport modeling are presented in this chapter.

Contents
18.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.2 Variably Saturated Water Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.2.1 Water Retention and Hydraulic Conductivity . . . . . . . . . . . . . . . . . .
18.2.2 Mass Balance Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.2.3 Preferential Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

788
790
791
797
799

D. Mallants (B)Belgian Nuclear Research Centre SCK•CEN, Mol, Belgium
e-mail: dmallant@sckcen.be

F.A. Swartjes (ed.), Dealing with Contaminated Sites,
C Springer Science+Business Media B.V. 2011
DOI 10.1007/978-90-481-9757-6_18, 

787

788

D. Mallants et al.

The Evapotranspiration Process . . . . . . . . . . . . . . . . . . . . . . . . . .
Penman-Monteith Equation for Evapotranspiration . . . . . . . . . . . . . .
FAO56 Reference Evapotranspiration . . . . . . . . . . . . . . . . . . . . . .
Root Water Uptake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Application: Numerical Simulations of Variably Saturated Flowin a Soil Profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.3 Contaminant Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.3.1 Transport Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.3.2 Advection-Dispersion Equations . . . . . . . . . . . . . . . . . . . . . . . . .
18.3.3 Nonequilibrium Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.3.4 Stochastic Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.3.5 Multicomponent Reactive Solute Transport . . . . . . . . . . . . . . . . . . .
18.3.6 Multiphase Flow and Transport . . . . . . . . . . . . . . . . . . . . . . . . . .
18.4 Analytical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.4.1 Analytical Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.4.2 Existing Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.5 Numerical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.5.1 Numerical Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.5.2 Existing Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.6 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18.2.418.2.518.2.618.2.718.2.8

799
801
802
803
804
808
809
813
817
821
826
827
828
828
828
829
829
831
839
841

18.1 Introduction
Many contaminants are finding their way into the subsurface environment, either
intentionally applied, for example in agricultural operations, or unintentionallyreleased from leaking industrial and municipal waste disposal sites, or from othersources (Šimůnek and Van Genuchten 2006). 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 contamination.
The same is true for salts and toxic trace elements, which are often an unintended
consequence of irrigation in arid and semi-arid regions. Agriculture also increasingly uses a variety of pharmaceuticals and hormones in animal production manyof which, along with pathogenic microorganisms, are being released to the environment through animal waste. Similar problems arise with point-source contaminationfrom industrial and municipal waste disposal sites, leaking underground storage
tanks, chemicals spills, nuclear waste repositories, and mine tailings, among othersources.
Mathematical models should be critical components of any effort to understand and predict site-specific subsurface water flow and contaminant transport
processes. For example, models can be helpful tools for designing, testing andimplementing soil, water and crop management practices that minimize soil andwater contamination. Models are equally needed for designing or remediating industrial waste disposal sites and landfills, or for long-term stewardship of nuclear waste

18

Leaching of Contaminants to Groundwater789

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 canbe eliminated before it contaminates groundwater, surface and subsurface water
resources. Sources of contamination often can be more easily remediated in thevadose zone, before contaminants reach the underlying groundwater. The focus of
this chapter thus will be on conceptual and mathematical descriptions of water flow
and especially transport processes in the predominantly unsaturated or variablysaturated vadose zone. The vadose zone is defined here as the zone between the land
surface and the permanent (seasonal) groundwater table. The vadose zone is usually only partially saturated, although saturated regions may exist, such as perchedwater above a low-permeable fine-textured (clay) layer or a saturated zone behind
the infiltration front during or after a high-intensity rainfall event.
Since the transport of contaminants is closely linked with the water flux in soils
and rocks making up the vadose zone, any quantitative analysis of contaminanttransport must first evaluate water fluxes into and through the vadose zone. Water
typically enters the vadose zone in the form of rainfall or irrigation (Fig. 18.1), orby 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 bysurface runoff, or will accumulate at the soil surface until it evaporates back to theatmosphere or infiltrates into the soil. Some of the water that infiltrates into the soil
profile may be taken up by plant roots and eventually returned to the atmosphere byplant transpiration. The processes of evaporation and transpiration are often combined into the single process of evapotranspiration. Only water that is not returnedto the atmosphere by evapotranspiration may percolate to the deeper vadose zone
and eventually reach the groundwater table. If the water table is close enough tothe soil surface, the process of capillary rise may move water from the groundwatertable through the capillary fringe towards the root zone and the soil surface.
Because of the close linkage between water flow and contaminant transport, we
will first briefly focus on the physics and mathematical description of water flowin the vadose zone (Section 18.2). An overview is given of the governing equations for water flow, while a comprehensive example is used to illustrate watercontent and pressure head distributions in single- and multi-layered soil profiles
following steady-state infiltration (Section 18.2.8). This is followed by a discussionof the equations governing contaminant transport (Section 18.3) where we review
the standard equilibrium formulations for contaminant transport (Section 18.3.2) aswell as alternative non-equilibrium models (Section 18.3.3). We also briefly discusspossible formulations for colloid-facilitated transport (Section 18.3.3.3), stochastic
approaches for contaminant transport (Section 18.3.4), and multicomponent geochemical transport (Section 18.3.5). This is followed by a discussion of analytical(Section 18.4) and numerical (Section 18.5) approaches for solving the governing flow and/or transport equations, and an overview of computer models currently

790

D. Mallants et al.

Fig. 18.1 Schematic of water fluxes and various hydrologic components in the vadose zone(Šimůnek and Van Genuchten 2006)

available for simulating vadose zone flow and transport processes (Sections 18.4.2
and 18.5.2). Several sections of this text are adapted from a recent book chapter byŠimůnek and Van Genuchten (2006).

18.2 Variably Saturated Water Flow
In this Section we briefly present the equations governing unsaturated water flowin the subsurface. 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 assuming the presence of spatially variable soil hydraulic properties (e.g., as dictated by the presence of different soil horizons, but possibly also
varying laterally). Unfortunately, the vadose zone can be extremely heterogeneousat 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

18

Leaching of Contaminants to Groundwater791

and dissolved contaminants through macropores (e.g., between soil aggregates, orcreated by earthworms or decayed root channels, see e.g. Mallants et al. (1996a)) orrock fractures, with much of the water bypassing (short-circuiting) the soil or rockmatrix. However, many other causes of preferential flow exist, such as flow instabilities caused by soil textural changes or water repellency (Hendrickx and Flury 2001;
Ritsema and Dekker 2000; Šimůnek et al. 2003), 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.

18.2.1 Water Retention and Hydraulic Conductivity
18.2.1.1 Water Retention
Above the groundwater table, a zone of a few to several tens of meters occurs wherepart of the pore space is occupied by the air phase. This is the unsaturated or vadosezone, of which the upper part typically contains a soil profile. When the soil becomesdrier owing to internal drainage and/or evapotranspiration, air replaces water first in
the coarse parts of the pore space and at successively lower (negative) values of thewater potential (see further) also in the finer pores. In the unsaturated zone, water isheld in the soil pores by so-called surface-tension forces. In other words, capillary
forces (and to a lesser extent also adsorption) bind water to solids. This leads to theexistence of a negative pressure, also referred to as the suction, tension, or matrichead (by definition pressures less than atmospheric are considered negative).
Capillary forces are the result of a complex set of interactions between the solid(particles) and liquid (pore water) phases involving the surface tension of the liquid phase, the contact angle between the solid phase and the liquid phase, and thediameter of pores. As a consequence of these forces, water will rise to a height H[L] when a capillary tube of radius R is placed into a water reservoir open to theatmosphere (Fig. 18.2). This capillary rise is given by the Laplace equation:H=

2σ cos γ
ρw gR(18.1)

where σ is the surface tension [MT−2 ] (7.27×10−2 kg/s2 at 20◦ C), γ is the contact
angle, ρ w is the density of the liquid phase [ML−3 ] (998 kg/m3 at 20◦ C), and g is
the gravitational acceleration [LT−2 ] (9.81 m/s2 ). For water at 20◦ C with γ = 0, Eq.
(18.1) can be simplified to H ∼= 1.5×10−5 /R (with both H and R in m).
Since a soil can be viewed as a complex system containing pores of various diameters, water in those pores will rise to different heights (Fig. 18.3), and hence will
be held with different potential energies. Because each soil has a different distribution of pore sizes, the distribution of water above the water reservoir will also be
different. This simple conceptual model assumes that soil pores can be represented

792

D. Mallants et al.

Fig. 18.2 Capillary rise H ofwater in a cylindrical tube ofradius R with a contact
angle γFig. 18.3 Capillary rise ofwater in a set of cylindricaltubes of various radiiby an equivalent bundle of capillaries with identical water retention properties as
the real soil. Such an approach with capillary tubes and a water reservoir can beused to evaluate the water content distribution in a soil above the groundwater tableat equilibrium. The experimental curve that describes this relationship between thewater content versus the height above the water table is called here the water retention curve. Many other names may be found in the literature, including pore water
characteristic curve, capillary pressure-saturation relationship, and pF curve. The
retention curve historically was often given in terms of pF, which is defined as thenegative logarithm of the absolute value of the pressure head in cm. The water retention curve provides information on how tightly water is held in pores and how muchwork would need to be exerted to extract it from the different pores. The water

18

Leaching of Contaminants to Groundwater793

retention curve thus characterizes the energy status of water in the soil, and is one of
the two soil hydraulic functions necessary to describe the status and movement ofwater in the vadose zone. Since capillary rise depends on the radius of a particular
capillary, the retention curve may also be interpreted as a curve that characterizes
the distribution of pores of different radii in the soil.
The dependency of the fraction of water-filled porosity (i.e., the water content
θ ) on the soil matric potential ψ m is formalized in the soil water retention curve,θ (ψ m ). This relation is of fundamental importance for the hydraulic characterizationof a soil since it relates an energy density (potential) to a capacity quantity (watercontent). One may view the water retention curve as the curve that results whenwater is slowly removed from an initially water saturated soil until the soil becomes
air saturated. This is the desorption curve. The adsorption (or absorption) curve
describes the reverse process. Note that air can only enter the porous medium after
the matric potential has fallen below a certain value ψ m0 , the so-called air-entry
value. This value is determined by the largest pore of the porous medium open to air.
Once air has entered the porous medium, the water content decreases monotonicallywith increasingly negative matric potentials ψ m .
Rather than using the matric potential in the parameterization of the water
characteristic, the pressure head h is often used instead:
h=−

m
ρw g

(18.2)

where h has SI units of meter (m) and dimensions of L; h represents the energy stateof pore water and is expressed as energy per unit weight.
Figure 18.4 shows water retention curves for three different textural classes, i.e.,
for a clay, a loam, and a sand, as given by Carsel and Parish (1988) using the soil
hydraulic parameters given in Table 18.1 (to be discussed later). Notice that the sandloses its water relatively quickly (at small negative pressure heads) and abruptly
above the water table (the pressure head h = 0 cm at the water table), while themore fine-textured loam and especially the clay soil lose their water much more
gradually. This reflects the pore-size distribution of a particular soil textural class.
While the majority of pores in coarse-textured soils (such as sand and gravel) haverelatively large diameters and thus drain at relatively small negative pressures, themajority of pores in fine-textured soils (such as clays, silty clays and clay loams) do
not drain until very large tensions (negative pressures) are applied.
Commonly used mathematical expressions for the retention curve, θ (h), arethe Van Genuchten (Van Genuchten 1980a) and Brooks and Corey (1964) equations since they permit a relatively good description of θ (h) for many soils using
only a limited number of parameters. The Van Genuchten soil moisture retention
characteristic is defined as:θ (h) = θr +θs − θr(1 + |ah|n )m(18.3)

794

D. Mallants et al.
4
Sand
Loam
Clay
Log10(matric head [m])Fig. 18.4 Example ofretention curves for sand,loam and clay textural classes(based on Carsel and Parish(1988) soil hydraulicparameters, Table 18.1)

2

0

–2
0

0.1

0.2

0.3

0.4

0.5

Water content [–]Table 18.1 Soil hydraulic parameters for the analytical functions of Van Genuchten (1980a) for
the twelve textural classes of the USDA soil textural triangle (Carsel and Parish 1988)
Textural classθ r [L3 L−3 ]θ s [L3 L−3 ]α [cm−1 ]n [−]Ks [cm d−1 ]SandLoamy Sand
Sandy Loam
Loam
Silt
Silty Loam
Sandy Clay Loam
Clay Loam
Silty Clay Loam
Sandy Clay
Silty Clay
Clay

0.045
0.057
0.065
0.078
0.034
0.067
0.100
0.095
0.089
0.100
0.070
0.068

0.430
0.410
0.410
0.430
0.460
0.450
0.390
0.410
0.430
0.380
0.360
0.380

0.145
0.124
0.075
0.036
0.016
0.020
0.059
0.019
0.010
0.027
0.005
0.008

2.68
2.28
1.89
1.56
1.37
1.41
1.48
1.31
1.23
1.23
1.09
1.09

712.8
350.2
106.1
24.96
6.00
10.80
31.44
6.24
1.68
2.88
0.48
4.80

where θ r is the residual water content [L3 L−3 ], θ s is the saturated water content[L3 L−3 ], and α [L−1 ], n [−] and m (= 1−1/n) [–] are shape parameters. The
dependency of the Van Genuchten model on the parameters α and n is shown in
Fig. 18.5. The curves in this figure are plotted in terms of effective saturation given
by Se =(θ –θ r )/(θ s −θ r ).
As a first approximation and on intuitive grounds, θ s = η (total porosity) andθ r = 0. In reality, however, the saturated water content θ s of soils is often smaller
than the total porosity η because of entrapped and dissolved air. The residual watercontent θ r is generally larger than zero, because of the presence of adsorbed water.

18

Leaching of Contaminants to Groundwater

795
10000

Matric head h (m)

10000
1000
100
10

α = 0.1

1000

1

100

n = 1.25

10

10

1

1

0.1

0.1

2
5

0.01

0.01
0

0.2

0.4

0.6

0.8

Saturation degree Se

1

0

0.2

0.4

0.6

0.8

1

Saturation degree Se

Fig. 18.5 Dependency of soil water characteristic curves on Van Genuchten parameters α (keeping
n fixed at 1.25) (left) and n (keeping α fixed at 1) (right). Saturation degree Se is used rather thanwater content

While some definitions of θ r have been proposed (e.g., the water content where
the hydraulic conductivity becomes first zero during drying), θ r in practice is besttreated as a fitting parameter without any physical significance.
Another parametric description of the moisture retention characteristic is theequation of Brooks and Corey (1964):θ (h) =

8 8n
8 8
θr + (θs − θr ) 8 hhe 8 h < heθsh ≥ he(18.4)

where he is the air-entry value [L] and n is the pore-size distribution index thatcharacterizes the width of the pore-size distribution (or the steepness of the retentionfunction). Like n in Eq. (18.3), the value of λ is relatively large for soils with arelatively uniform pore-size distribution, but small for soils having a wide range ofpore sizes.
A large number of experimental techniques can be used to measure the water
retention curves. These techniques include methods using a hanging water column,
pressure cells, pressure plate extractors, suction tables, soil freezing, and many otherapproaches (Dane and Hopmans 2002). Alternatively, soil hydraulic parameters may
be obtained with pedotransfer functions using information on textural class or textural fractions (and some additional information). Pedotransfer functions are statistical
relationships that use easily available soil data to estimate the soil hydraulic parameters (Bouma 1989; Leij et al. 2002; Pachepsky and Rawls 2004; Pachepsky et al.
2006).
18.2.1.2 Hydraulic Conductivity
Another important soil hydraulic property is the unsaturated hydraulic conductivity
function. The hydraulic conductivity characterizes the ability of a soil to transmit
water, and as such is inversely related to the resistance to water flow. The hydraulic

796

D. Mallants et al.

conductivity depends on many factors, including the pore-size distribution of theporous medium, and the tortuosity, shape, roughness, and degree of interconnectedness of pores. The hydraulic conductivity decreases considerably as soil becomesunsaturated and less pore space is filled with water. The unsaturated hydraulic conductivity function gives the dependency of the hydraulic conductivity on the watercontent, K(θ ), or pressure head, K(h). The conceptual model that views the soil
as a bundle of capillaries of different radii, as used earlier to explain water retention properties (Fig. 18.3), may be used also to evaluate the hydraulic conductivity
function. The ability of each capillary tube to conduct water can be calculated using
Poiseuille’s flow equation. By adding the conductivity of all capillaries that are filledwith water at a particular water content or pressure head, one obtains the hydraulicconductivity of the complete set of capillaries, and consequently of the soil itself.
Figure 18.6 presents examples of hydraulic conductivity functions for the sand,loam, and clay textural classes, again using the Carsel and Parish (1988) parametervalues as listed in Table 18.1. The hydraulic conductivity curves in Fig. 18.6 are presented as functions of both the pressure head (left) and water content (right). Noticethat the hydraulic conductivity at saturation is significantly larger for coarse-textured
soils (sands and gravels) than for loams and clays. This difference is often several orders of magnitude. Also notice that the hydraulic conductivity may decrease
several orders of magnitude as the soil becomes unsaturated. This decrease, when
expressed as a function of the pressure head (Fig. 18.6, left), is much more significant for sands than for loams, and even more so than for clays. The decrease forcoarse-textured soils (sands and gravel) is so large that at a certain pressure head thehydraulic conductivity becomes smaller than for heavy-textured soils (clays). These
properties of the hydraulic conductivity function are often used in the design of

4
Sand
Loam
Clay

2

Log10(Hydraulic conductivity (cm/d))

Log10(Hydraulic conductivity (cm/d))

4

0
–2
–4
–6
–8

–10

2
0
–2
–4
–6

Sand
Loam
Clay

–8

–10
–2

0

2

Log10 (Matric head [m])

4

0

0.1

0.2

0.3

0.4

0.5

Water content [-]Fig. 18.6 Example of hydraulic conductivity curves for the sand, loam and clay soil texturalclasses (based on Carsel and Parish (1988) soil hydraulic parameters, Table 18.1)

18Leaching of Contaminants to Groundwater

797

engineered structures, such as capillary barriers (finer-textured soils above coarsertextured soils) in landfill covers to divert water from flowing through the underlying
waste (e.g., Mallants et al. 1999), or having narrow sand or gravel layers at the soil
surface to prevent or limit evaporation.
Similarly as for the water retention curve, analytical models are often used also
for the hydraulic conductivity function. Many of these functions were derived using
the pore-size distribution models of Burdine (1953) or Mualem (1976) in combination with one of the retention functions given earlier. The Brooks and Corey (1964)
retention (Eq. (18.4)) is commonly associated with Burdine’s pore-size distribution
model, leading to the hydraulic conductivity function:
K = Ks Se2/λ+l+1(18.5)

in which Ks is the saturated hydraulic conductivity [LT−1 ], Se =(θ –θ r )/(θ s –θ r ) is
effective saturation, and l a pore-connectivity parameter assumed to be 2.0 in theoriginal study of Brooks and Corey (1964).
The Van Genuchten (1980a) retention function (Eq. (18.3)) is similarly coupledmostly with the model of Mualem (1976) to giveK(h) = Ks Sel [1 − (1 − Se1/m )m ]2(18.6)

m = 1 − 1/n, n > 1(18.7)

where

The pore-connectivity parameter l in (Eq. (18.7)) was estimated by Mualem(1976) to be about 0.5 as an average for many soils. However, more recently, Schaapand Leij (2000) recommended using l equal to –1 as an appropriate value of mostsoil textures.
The saturated hydraulic conductivity Ks can be measured in the laboratory using
either the constant or falling head method, and in the field using single or double
ring infiltrometers, a constant head well permeameter, the auger-hole method, or the
piezometer method (Dane and Topp 2002). Many laboratory and field methods arealso available for measuring unsaturated hydraulic conductivities (Dane and Topp
2002).

18.2.2 Mass Balance Equation
Water flow in variably saturated rigid porous media (such as soils) is usuallyformulated in terms of a mass balance equation of the form
∂q
∂θ
=−
−S
∂t
∂z

(18.8)

798

D. Mallants et al.

where θ is the volumetric water content [L3 L−3 ], t is time [T], z is the spatial coordinate [L], q is the volumetric flux [LT−1 ], and S is a general sink/source term
[L3 L−3 T−1 ], for example to account for root water uptake (transpiration).
Equation (18.8) is often referred to also 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 thewater flux (i.e., fluxes in and out of some small 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 (q) to produce the governing
equation for variably saturated flow.
For a soil that can be saturated or unsaturated, the flux is given by the equation:q = −k(h)∂h+ K(h)∂z(18.9)

where K is the unsaturated hydraulic conductivity [LT−1 ], and h the pressure head.
Eq. (18.9) is often referred to as the Darcy-Buckingham equation. The hydraulicconductivity in this equation is a function of the pressure head, h. In saturated systems, the conductivity becomes independent of the pressure head, in which case theequation reduces to Darcy law as:
q = −Ks

∂h
+ Ks
∂z

(18.10)where Ks is the saturated hydraulic conductivity [LT−1 ]. The Darcy-Buckingham
equation hence 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 (Eq. (18.8)) with Darcy-Buckingham’s law (Eq.
(18.9)) leads to()∂∂h∂θ (h)=K(h)+ K(h) − S(h)
∂t
∂z
∂h

(18.11)

which was first formulated by Richards (1931) and is now generally referred to as
the Richards equation. This partial differential equation is the equation governing
water flow in the variably-saturated vadose zone. Because of its strongly nonlinear
makeup, only a relatively few simplified analytical solutions can be derived. Most
practical applications of Eq. (18.11) require a numerical solution, which can beobtained using a variety of methods such as finite differences or finite elements.
Equation (18.11) is generally referred to as the mixed form of the Richards equation
since it contains two dependent variables, i.e., the water content and the pressure
head. Various other formulations of the Richards equation are possible.

18

Leaching of Contaminants to Groundwater799

18.2.3 Preferential Flow
Increasing evidence exists that water in many field soils does not move according
to the uniform flow pattern typically predicted with the Richards equations (Fluryet 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 andcontaminants may move preferentially, while bypassing a large part of the matrix
pore-space, called preferential flow. Preferential flow and transport processes aremost frustrating in terms of hampering accurate predictions of contaminant transport
in soils and fractured rocks (Šimůnek and Van Genuchten 2006). Contrary to uniform flow, preferential flow results in irregular wetting of the soil profile as a directconsequence 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 contaminants move along certain pathways, while bypassing afraction of the porous matrix. Water and contaminants for these reasons can movequickly to far greater depths, and much faster, than would be predicted with theRichards 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), andfunneling 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) areusually caused by textural differences and other factors at scales significantly larger
than the pore scale, macropore flow and transport are usually generated at the poreor slightly larger scales, including scales where soil structure first manifests itself,
mostly the pedon scale (Šimůnek 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/or multipermeability models (Gerke and Van Genuchten 1993; Gwo et al. 1995; Jarvis
et al. 1998; Mallants et al. 1996b, 1996c; Pruess and Wang 1987; Richards 1931;Šimůnek et al. 2003, 2005). While single-porosity models assume that a single poresystem exists that is fully accessible to both water and contaminant, dual-porosityand dual-permeability models 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 soil or rock matrix is stagnant, dual-permeability models allow
also for water flow within the matrix.

18.2.4 The Evapotranspiration Process
Evapotranspiration (ET) is the process by which water is transported from theearth surface (i.e., the plant-soil system) to the atmosphere by evaporation (E) from

800

D. Mallants et al.

surfaces (soils and wet vegetation) and by transpiration (T) from plants throughstomata in the plant leaves (Fig. 18.1). Evaporation is the physical process by whichwater is transformed into water vapour and is removed from the evaporating surface by mass transfer. Transpiration is the process by which water in plant tissues
is transformed into water vapour and removed towards the atmosphere. Important
weather parameters affecting ET are radiation, air temperature, humidity and windspeed (see further). Another important factor is the amount of water available at theevaporating surface (e.g., the soil surface), or in the soil for uptake by plant roots.
The available water is determined by such factors as soil type (texture), depth toground water, irrigation and management practices. Root water uptake is negatively
influenced by water logging (i.e., when the soil is completely or almost completely
saturated), and pore water salinity. Crop characteristics such as crop type, variety,development stage, crop height, crop roughness, reflection, ground cover and crop
rooting characteristics all influence the resistance to transpiration, and hence ET.
To remove the effects due to soil type, management and crop factors on calculations of the evaporative demand of the atmosphere, the evapotranspiration rate is
generally calculated for a reference surface not short of water. This ET is called
the reference evapotranspiration, ETo , and usually calculated following guidelinesof the FAO56 paper by Allen et al. (1998) using the Penman-Monteith equation(Monteith 1965) (see Section 18.2.6). The reference surface is a hypothetical reference crop (grass) with an assumed height of 0.12 m, a fixed surface resistance of
70 s m−1 and an albedo of 0.23. ETo thus only depends on climatic parameters andcan be considered as a climatic parameter expressing the evaporating power of theatmosphere at a specific location and time of the year (Allen et al. 1998).
Since actual ground cover, canopy properties and aerodynamic resistance of thecrop are different from those used for calculating ETo , the evapotranspiration rate
under standard conditions (i.e., of a large field under excellent agronomic and porewater conditions) for a specific crop, ETc , is required for specific applications.
ETc can be obtained either by using specific crop parameters (e.g., albedo, aerodynamic and canopy surface resistances) in the Penman-Monteith equation, or bymultiplying ETo with a crop coefficient Kc . The crop coefficient incorporates four
primary characteristics that distinguish a specific crop from the reference crop (cropheight, albedo, canopy resistance, and evaporation from the soil) (Allen et al. 1998)
and is determined by crop type, climate, soil evaporation and crop growth stages.
Consequently, Kc coefficients change during a growing season.
Actual evapotranspiration under field conditions, ETa , takes into account nonideal (non-standard) conditions such as water or salinity stress. Under dry soil
conditions, water flow in a soil can be too slow to satisfy the evaporative demand.
Similarly as for transpiration, very dry or very wet soil conditions, or high salt concentrations, impose water and salinity stresses and reduce root water uptake. Allenet al. (1998) related ETa to ETo by means of a water stress coefficient and/or by
adjusting Kc for all kinds of stresses. In this chapter, a mechanistic approach todetermine ETa is described using water flow in the soil and water uptake by plantroots, the latter incorporating empirical water stress reduction functions (see Section
18.2.5).

18

Leaching of Contaminants to Groundwater801

18.2.5 Penman-Monteith Equation for Evapotranspiration
For evapotranspiration to occur, three conditions are needed in the soil-plantatmosphere system (Jensen et al. 1990):
1. a supply of water must be available;
2. energy must be available to convert liquid water into vapour water;
3. a vapour pressure gradient must exist to create a flux from the evaporating surface
to the atmosphere.
Penman (1948) proposed a combination method by introducing an energy balance (condition 2) and a mass transfer term in an aerodynamic formula (condition
3) into a single equation to calculate ET. Penman’s method was developed to calculate E as open water evaporation. Written as the weighted sum of the ratesof evaporation due to net radiation, Er (MJ m−2 d−1 ), and turbulent mass transfer, Ea (MJ m−2 d−1 ), Penman’s equation for the evaporative latent heat flux, λE(MJ m−2 d−1 ), is:
λE =

+γ

Er +

γ
Ea
+γ

(18.12)where λ is the latent heat of vaporization (MJ kg−1 ), is the slope of the vapour
pressure curve (kPa ◦ C−1 ) and γ is the psychrometric constant (kPa ◦ C−1 ). Er isgiven by:
Er = Rn − G(18.13)where Rn is the net radiation flux (MJ m−2 d−1 ) and G is the sensible heat flux into
the soil (MJ m−2 d−1 ), and Ea by:
Ea = Wf (es − ea )(18.14)

where es and ea are the saturation and actual vapour pressures, respectively (kPa),(es –ea ) is the saturation vapour pressure deficit, and Wf is a wind function (MJ d−1
kPa−1 ). A linear wind function was found to be adequate, defined as (Allen 2001):
Wf = Kw (aw + bw u2 )(18.15)where Kw is a units conversion factor [6.43 for ET0 in mm d−1 ], and aw and bw are
empirical wind function coefficients often obtained by regional or local calibration.
The Penman method to estimate the evaporation from open water is then:
E=

1
λ

+γ

(Rn − G) +(γ
Kw (aw + bw u2 )(es − ea )+ γ)

(18.16)

802

D. Mallants et al.

Penman (1948) derived this equation for open water evaporation. Evaporationfrom bare soil, wet soil and grasses is obtained as a fraction of E. Bulk surface
resistances from the soil and crop is not explicitly accounted for, but are incorporated in the wind function. Resistance factors are incorporated in Penman-based
equations to include the resistance of vapour flow through stomata openings, total
leaf area and the soil surface (the surface resistance, rs ), and the resistance fromthe vegetation upwards involving friction from air flowing over vegetative surfaces(aerodynamic resistance, ra ) (Allen et al. 1998). The Penman-Monteith equation(Monteith 1965) for evaporation from bare soil, wet soil and grasses, ET, is given in
the ASCE standard form (Allen et al. 1998) as:⎞⎛(es −ea )1 ⎝ (Rn − G) + Ktime ρa cp ra ⎠


ET =
λ
+ γ 1 + rs

(18.17)rawhere ρ a is the mean air density at constant pressure (kg m−3 ), cp is the specific heatof air (MJ kg−1 ◦ C−1 ), ra and rs are the aerodynamic and (bulk) surface resistances,respectively (s m−1 ), and Ktime is a units conversion factor (86400 s d−1 when ET isexpressed in mm d−1 ). For a more detailed definition and discussion of Eq. (18.17)
and its parameters the reader is referred to Allen et al. (1998).

18.2.6 FAO56 Reference Evapotranspiration
Allen et al. (1998) calculated the reference evapotranspiration, ET0 , using the ASCE
Penman-Monteith equation (Eq. (18.17)) for a hypothetical reference surface or reference crop defined as “a cropped soil with an assumed crop height of 0.12 m, a fixedsurface resistance of 70 s m−1 and an albedo of 0.23”, with climatological parameters measured at a reference level of 2 m above the soil surface. The popularly used
FAO56 Penman-Monteith equation is defined as (Allen et al. 1998):

ET0 =

900
u2 (es − ea )T + 273+ γ (1 + 0.34u2 )0.408 (Rn − G) + γ(18.18)To calculate ET0 using the FAO56 Penman-Monteith equation on a daily basis,such as implemented for instance in the ET-REF program (Allen 2001), site-specificdata (altitude above sea level and latitude) and climatological data (temperature,humidity, radiation, and wind speed) are required. The altitude above sea level determines the local average value of atmospheric pressure. The latitude is needed tocompute extraterrestrial radiation.
The potential evapotranspiration of a particular crop or vegetation, ETc , is thenobtained by multiplying ET0 with a crop coefficient, Kc . ETc is divided between

18

Leaching of Contaminants to Groundwater803

evaporation of the intercepted water (Ep,I ), potential soil evaporation (Ep,s ) andpotential transpiration (Tp ) using procedures outlined below.
Daily values of the actual soil evaporation rate, Ea,t [L], and actual transpiration
rate, Ta,t [L], may be calculated using HYDRUS-1D (Šimůnek et al. 2005). The
input variables for HYDRUS-1D are time series of daily values for the throughfall (e.g., precipitation or irrigation reaching the soil surface), T, potential soil
evaporation, Ep,s , and potential transpiration, Tp .

18.2.7 Root Water Uptake
The sink term S in Eq. (18.8) is defined as the volume of water extracted from a unit
volume of soil per unit time by the roots. The potential root water uptake rate Sp (z)is often obtained by multiplying a normalized water uptake distribution b(z) [L−1 ]with the potential transpiration rate Tp [LT−1 ] as follows:
Sp (z) = b(z)Tp(18.19)The function b(z) may be obtained from the root distribution with depth:b (z)b(z) = 9b (z)(18.20)LR

where b (z) is the root distribution function and LR is the soil root zone. Note thatb’(z) can be of any form. Tp depends on climate conditions and vegetation (leaf areaindex, crop coefficients − see further). The actual root water uptake rate S(z) may
be obtained by multiplying Sp (z) with a root water stress response function (e.g.,
Feddes et al. 1978; Van Genuchten 1987) to account for a possible reduction in root
water uptake due to water stress conditions in the soil profile:
S(h, z) = α(h)Sp = α(h)b(z)Tp(18.21)

where α(h) is the water stress response function as a function of the pressure head.
To obtain the actual transpiration rate Ta of the vegetation, the actual root water
uptake S(h,z) (Eq. (18.21)) is integrated over the rooting depth::
Ta = Tpα(h, z)b(z)dz
LR(18.22)

804

D. Mallants et al.

18.2.8 Application: Numerical Simulations of Variably Saturated
Flow in a Soil Profile
Analysis of water flow in soils can be based on experimental exploration and/or
numerical simulation. Experimental studies yield valuable insight into the physical processes, but their implementation is often complicated, time consuming, andcostly. As an alternative, or a useful complement, numerical simulations of the water
regime can provide necessary information to analyse flow phenomena at a very highspatio-temporal resolution and at reasonable cost. Also, simulations may be used
to plan and interpret experiments and to predict future flow behaviour for complex
conditions.
In the following examples, numerically simulated water flow will be discussed
for different conditions of increasing complexity. We start the discussion with onedimensional flow in a hypothetical single-layer (homogeneous) soil profile, and thenmove to a two-layered soil. Note that we will not calculate water uptake by roots
explicitly, but rather apply a net rainfall rate, which is equal to the total rainfall
rate minus the actual evapotranspiration (runoff is assumed zero here, but could be
accounted for if needed).
The simulations for the hypothetical soil made use of the Van GenuchtenMualem hydraulic functions without hysteresis. Figure 18.7 shows the soil water
retention and hydraulic conductivity functions. All simulations were carried out
with the HYDRUS-1D finite element code (Šimůnek et al. 2005). As top boundary condition for the hypothetical soil profile we implemented a constant rainfall
rate of 2×10−6 m/s (or 7.2×10−3 m/h). This value corresponds to one tenth of
the saturated hydraulic conductivity of the sand material. As bottom boundary weassumed a groundwater table 4 m below surface. This condition is mathematically
implemented by fixing the pressure head equal to zero at the water table. As initialcondition for the pressure head we used h = – z – 4 m, where z is the depth below
surface (the vertical axis is positive upward, with the reference depth z = 0 at thesoil surface, see Fig. 18.8). At the soil surface the pressure head initially is hence
0 − 4 = −4 m, with h decreasing linearly to h = −(−4) − 4 = 0 m at the groundwater table. The soil profile was discretized into 0.05-m thick elements, which resulted
in 81 nodes for the one-dimensional simulation.
18.2.8.1 Single-Layer Soil
The numerically calculated infiltration process for a homogeneous sand is shown in
Fig. 18.9 in terms of pressure head (h − z) and water content (θ − z) distributions
versus depth. At t < 0, the initial pressure head is equal to the height above thewater table, i.e. h = 0 at the water table and h = −4 at the soil surface (also see
Fig. 18.8). The corresponding water content profile has the same shape as the soil
water retention curve. At t > 0, water starts to infiltrate and the infiltration front
moves downward at a nearly constant rate. After 70 h following the start of theinfiltration, the water content in the profile reaches a steady-state value, while aconstant water flux is established in the entire soil profile.

Hydraulic conductivity, K(h) (m/s)

18

Leaching of Contaminants to Groundwater

805

1E–4

1E–4

1E–5

1E–5

1E–6

1E–6

1E–7

1E–7

1E–8

1E–8

1E–9

1E–9
0

1

2

3

4

0

0.2

5

0.6

0.8

1

5

Sand
Silt

4

log10(-h, cm)

0.4

Saturation degree, Se

log10(h, cm)

4

3

3

2

2

1

1

0

0
0

0.1

0.2

0.3

0.4

0.5

0

0.2

Water content, θ

0.4

0.6

0.8

1

Saturation degree, Se

Fig. 18.7 Hydraulic functions for a sand and a silt soilzh, z, H

Fig. 18.8 Variation of thecomponents h and z of thetotal head H along a verticalprofile of an unsaturated soil.
The total head H = −4 meverywhere–4 m

Soil surface

H

z
h

H

Reference level z = 0hzGroundwater tableD. Mallants et al.
0

0

–1

–1

–2

–2

–3

–3

–4

soil depth, z (m)soil depth, z (m)

806

–4
–4

–3

–2

–1

pressure head, h (m)

0

0

0.1

0.2

0.3

0.4

0.5

water content, θFig. 18.9 Numerically calculated pressure heads (left) and water contents (right) versus depth for
infiltration in a homogeneous sand assuming a water flux q of 2×10−6 m/s. The solid line is for
t = 0, while the distributions at other times (1, 5, 10, 20, 30, 40, 50, 60, 70 h) are indicated bydashes of decreasing length

Having a constant water content between 0 and 3 m would greatly facilitate any
sampling effort when the water content has to be measured in a deep homogeneous field soil profile. As water infiltrated, a less negative pressure head developed,
reaching a constant value of −0.46 m from the soil surface to close to the watertable. In other words, with a pressure head gradient dh/dz ≈ 0, the hydraulic gradient d(h+z)/dz becomes one, which corresponds to a unit gradient. The water fluxis then determined by gravity only, often referred to as gravity flow. In that case,the unsaturated hydraulic conductivity equals the applied water flux and the watercontent in the profile is equal to the water content of hydraulic conductivity atthat flux.
Results for the homogeneous silt soil using the same water flux of 2×10−6 m/s
are shown in Fig. 18.10. The flux is now approximately five times lower than thesaturated hydraulic conductivity. Compared to the sand, the degree of saturation
degree will now be much higher. In contrast to the sand, the initial water content
profile is now much smoother, consistent with the reduced nonlinearity of the soil
water retention curve. The higher initial (at t=0) degree of saturation for the siltmeans that a smaller volume of pore space has to be filled with water during theinfiltration process as compared to the sand. As a result, the infiltration front moves
faster downward in the silt than in the sand.

Leaching of Contaminants to Groundwater

807

0

0

–1

–1

–2

–2

–3

–3

–4

soil depth, z (m)soil depth, z (m)

18

–4
–4

–3

–2

–1

pressure head, h (m)

0

0

0.1

0.2

0.3

0.4

0.5

water content, θFig. 18.10 Numerically calculated pressure head (left) and water content (right) versus depth fora homogeneous silt, assuming a water flux q of 2×10−6 m/s. Solid line is for t = 0, while the
distributions at other times (1, 5, 10, 20, 30, 40, 50, 60, 70 h) are indicated by dashes of decreasing
length

18.2.8.2 Two-Layer Soil
Heterogeneous soil profiles are the rule rather than the exception. In general, two
or more layers (soil horizons) are present parallel to the surface. Although waterflow is not necessarily one-dimensional in such cases, one can often approximate
the infiltration process by a one-dimensional model.
The first example considers a sand over silt layer. Hydraulic properties, initialand boundary conditions are those used for the single-layer soil described above. For
the water content, the initial condition displays now a discontinuity at the interface
between both layers (Fig. 18.11). The pressure head, on the other hand, is continuous
across the interface. This behaviour follows from the Buckingham-Darcy equation(Eq. 18.9) which has a finite q only when dh/dz is finite, or when h is continuous
everywhere. Up to t = 30 h, the infiltration into the sand behaves identically as in thehomogeneous sand. The pressure head behind the infiltration front is −0.46 m andthe corresponding conductivity in the sand is 2×10−6 m/s. Owing to the continuityof the pressure head across the interface, the silt soil has a conductivity of 1×10−6m/s, which is two times smaller than in the sand. As a result, water flow is retardedin the silt and the water content above the interface increases. This in turn leads
to a higher K in the silt layer. At steady-state, a unit gradient condition is present

D. Mallants et al.
0

0

–1

–1

–2

–2

–3

–3

–4

soil depth, z (m)soil depth, z (m)

808

–4
–4

–3

–2

–1

0

0

0.1

0.2

0.3

0.4

0.5

water content,θ

pressure head, h (m)Fig. 18.11 Numerically calculated pressure head (left) and water content (right) versus depth for
infiltration in a two-layer (sand over silt) soil profile, assuming a flux q of 2×10−6 m/s. The solidline is for t = 0, while distributions at other times (1, 5, 10, 20, 30, 40, 50, 60, 70 h) are indicatedby dashes of decreasing length

in both layers, except in a fairly thin zone near the interface. The consequence of
having several soil layers on sampling density is obvious from Fig. 18.11: very
few samples are required within each homogeneous layer, but the zone above theinterface requires a more dense sampling scenario.

18.3 Contaminant Transport
Similarly as Eq. (18.8) for water flow, mathematical formulations for contaminanttransport are based on a mass balance equation of the form:
∂CT
∂JT
=−
−φ
∂t
∂z

(18.23)where CT is the total concentration of contaminant in all forms [ML−3 ], JT is
the total contaminant mass flux density (mass flux per unit area per unit time)[ML−2 T−1 ], and φ is the rate of change of mass per unit volume by reactions orother sources (negative) or sinks (positive) such as plant uptake [ML−3 T−1 ].
In its most general interpretation, (Eq. (18.23)) allows the contaminant to exist
in all three phases of the soil (i.e., soil gas, liquid phase, and solid phase). The

18

Leaching of Contaminants to Groundwater809

equation further allows for a broad range of transport processes (including advective
transport, diffusion, and hydrodynamic dispersion in both liquid phase and the soil
gaseous phase), while accommodating any type of chemical reaction that leads tolosses or gains in the total concentration.
Whereas the majority of contaminants in soils are usually present only in the liquid phase and solid phase, and as such are transported in the vadose zone mostly
only in and by water, some contaminants such as many organic contaminants (e.g.,
VOCs, BETX), ammonium, fumigants, carbon dioxide, methane, and naturally
occurring radioactive gases including radon can have a significant portion of theirmass in the gaseous phase, and are thus also subject to transport in the gaseous
phase. The total contaminant concentration CT can be described as (Jury et al.
1991):CT = ρb s + θ c + ag

(18.24)

where ρ b is the bulk density [ML−3 ], θ is the volumetric water content [L3 L−3 ], ais the volumetric air content [L3 L−3 ], and s [MM−1 ], c [ML−3 ], and g [ML−3 ] areconcentrations in the solid phase, liquid phase, and soil gas, respectively. The solid
phase concentration represents contaminants sorbed onto sorption sites of the solid
phase, but can also include contaminants sorbed onto colloids attached to the solid
phase or strained by the porous system, and contaminants precipitated onto or into
the solid phase.
Various chemical or biological reactions that lead to a loss or gain of contaminant
in the soil system may be represented by the reaction term φ of Eq. (18.23). This
includes radionuclide decay, biological degradation, and/or mineral precipitationdissolution. In simulation models these reactions are most commonly expressedusing zero- and/or first-order reaction rates:
φ = ρb sμs + θ cμw + agμg − ρb γs − θ γw − aγg(18.25)where μs , μw , and μg are first-order degradation constants for the solid phase, liquid phase, and soil gas [T−1 ], respectively, and γ s [T−1 ], γ w [ML−3 T−1 ], and γ g[ML−3 T−1 ] are zero-order production constants for the solid phase, liquid phase,and soil, respectively.

18.3.1 Transport Processes
For contaminants that are present in both liquid phase and soil gas, the total contaminant mass flux, JT , is the result of various transport processes in both of thesephases:JT = Jl + Jg

(18.26)

810

D. Mallants et al.

where Jl and Jg represent contaminant fluxes in the liquid phase and soil gas[ML−2 T−1 ], respectively. The three main transport processes that can be active inboth liquid phase and soil gas are molecular diffusion, hydrodynamic dispersion,and advection (often also called convection). The contaminant fluxes in the waterand gas phases are then the sum of fluxes due to these different processes:
Jl = Jlc + Jld + Jlh
Jg = Jgc + Jgd + Jgh(18.27)

where the subscripts c, d, and h denote convection (or advection), molecular
diffusion, and hydrodynamic dispersion, respectively.
18.3.1.1 Diffusion
Diffusion is a result of the random motion of molecules. This process causes a contaminant to move from a location with a higher concentration to a location with alower concentration. Diffusive transport can be described using Fick’s first law:
Jd = −θ D∂c∂c= −θ ξ (θ )Do∂z∂z(18.28)where D is liquid phase diffusion coefficient [L2 T−1 ], Do is the diffusion coefficientof the contaminant in free water [L2 T−1 ], and ξ is a tortuosity factor to accountfor the increased path lengths and decreased cross-sectional areas of the diffusing
contaminant in both phases (Jury and Horton 2004).
The parameter D is related to the effective diffusion coefficient Deff (diffusion
through the pore space of the porous medium) following Deff = D×θ . Values for
Do vary between 2×10−9 m2 /s for anions like Cl− and Br− and 0.78×10−9 m2 /s forcations like Ca2+ . Since contaminant diffusion in the liquid phase is severely hampered by both air and solid particles, the tortuosity factor as defined here increases
strongly with water content. Many empirical models have been suggested in theliterature 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 ofMillington and Quirk (1961) given by:ξ (θ ) =

θ 7/3
θs2

(18.29)where θ s is the saturated water content [L3 L−3 ].
Jin and Jury (1996) suggested that the following version of the Millington-Quirk
model may be in better agreement with experimental data:
D=

θ
2/3
θs

Do

(18.30)

18

Leaching of Contaminants to Groundwater811

As an example, taking θ s = 0.37 and Do = 10−9 m2 /s, and assuming saturated
or nearly saturated conditions (i.e. θ = θ s = 0.37), then D = 0.72×10−9 m2 /s.

18.3.1.2 Dispersion
Dispersive transport of contaminants results from the heterogeneous distributionof water flow velocities within and between different soil pores (Fig. 18.12, left).
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 (Fig. 18.12, left).
For this reason contaminants in the middle of a pore will travel faster than contaminants that are farther from the center. Since the distribution of contaminant ions
within a pore depends on their charge, as well as on the charge of pore walls,
some contaminants may move considerably faster than others. In some circumstances (e.g., in fine-textured soils), anion exclusion may occur. Since the anions
in such soils are located predominantly in the faster moving liquid in the centerof a pore away from the negatively charge solid surface, anionic contaminants
may travel faster than the average velocity of water (e.g., Nielsen et al. 1986).
Based on 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. Since soils consist of pores of many different
diameters, contaminant fluxes in pores of different diameters will be significantly
different, with some contaminants again traveling faster than others (Fig. 18.12,
right). Furthermore, contaminants may travel according to pathways of different
length. All these factors result in a bell-shaped distribution of velocities and thus ofarrival times, typical of a breakthrough curve.

slow velocity
fast velocity

shortpathwaylong
pathway

Fig. 18.12 Distribution of velocities in a single pore (left) and distribution of velocities in a poresystem (right) resulting in different arrival times of contaminants (modified from Šimůnek and Van
Genuchten 2006)

812

D. Mallants et al.

The above pore-scale dispersion processes lead to an overall (macroscopic)hydrodynamic dispersion process that can be described mathematically in approximately the same way as molecular diffusion using Fick’s first law. Adding thedispersion and diffusion processes leads then to:
Jh = −θ Dh∂c∂c= −θ (Dm + D)∂c∂z

(18.31)

where Dh is the hydrodynamic dispersion coefficient [L2 T−1 ] that accounts for both
molecular diffusion and mechanical dispersion (Fetter 1999), Dm is the mechanical dispersion coefficient [L2 T−1 ], and D is the liquid phase diffusion coefficient[L2 T−1 ]. The mechanical dispersion coefficient in one-dimensional systems has
been found to be approximately proportional to the average pore-water velocity v(= q/θ ) [LT−1 ], with the proportionality constant generally referred to as the (longitudinal) dispersivity λ (Biggar and Nielsen 1967). The discussion above holds for
one-dimensional transport; multi-dimensional applications require the use of a more
complicated dispersion tensor involving longitudinal and transverse dispersivities(e.g., Bear 1972).
Dispersivity is a transport parameter that is often obtained experimentally byfitting measured breakthrough curves with analytical solutions of the advectiondispersion equation (discussed further below). The dispersivity often changes withthe distance over which contaminants travel. Values of the longitudinal dispersivity
usually range from about 1 cm for relatively short, packed laboratory columns, toabout 5 or 10 cm for field soils. Longitudinal dispersivities can be significantly larger(up to hundreds of meters) for regional groundwater transport problems (Gelhar
et al. 1985). If no other information is available, a good first approximation is touse a value of one-tenth of the transport distance for the longitudinal dispersivity(e.g, Anderson 1984), and a value of one-hundred of the transport distance for thetransverse dispersivity when multi-dimensional applications are considered.
18.3.1.3 Advection
Advective transport refers to contaminants being transported with the moving fluid,
either in the liquid phase (Jlc ) or the soil gas (Jgc ), i.e.:
Jc = qc(18.32)where we now use a more general notation by omitting the subscripts l and g.
Advective transport generally is dominant in the liquid phase. Advective transport
may also occur in the gaseous phase, but is often neglected since its contribution is
generally negligible compared to gaseous diffusion.
The total contaminant flux density in both the liquid phase and soil gas is
obtained by incorporating contributions from the various transport processes into
Eq. (18.27) to obtain

18

Leaching of Contaminants to Groundwater813

J = qc − θ Dh∂c∂z

(18.33)where Dh is the hydrodynamic dispersion coefficient [L2 T−1 ] that accounts for both
molecular diffusion and mechanical dispersion.
Mechanical dispersion in most subsurface transport problems dominates molecular diffusion in the liquid phase, except when the fluid velocity becomes relatively
small or is negligible. Diffusion dominated transport occurs in low permeability
media, such as clays, rock matrices, and man-made structures such as concrete.

18.3.2 Advection-Dispersion Equations
18.3.2.1 Transport Equations
The equation governing transport of dissolved contaminants in the vadose zone
is obtained by combining the contaminant mass balance (Eq. (18.23)) with equations defining the total concentration of the contaminant (Eq. (18.24)) and thecontaminant flux density (Eq. (18.33)) to give
∂
∂
∂(qc)
∂c
∂g
∂(ρb s + θ c + ag)= (θ Dh ) + (aDsg ) −
−φ
∂t
∂z
∂z
∂z
∂z
∂z

(18.34)

where Dh and Ds g are the hydrodynamic dispersion coefficient in the liquid andgaseous phases [L2 T−1 ], respectively.
Several alternative formulations of Eq. (18.34) can be found in the literature. For
example, for one-dimensional transport of non-volatile contaminants, Eq. (18.34)
simplifies to∂(θ Rc)∂∂(ρb s + θ c)
∂c
=
=
θ Dh
∂t
∂t
∂z
∂z

−

∂(qc)
−φ
∂z

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

R=1+

ρb ds(c)θ dc(18.36)

For transport of inert, non-adsorbing contaminants during steady-state water flow
a further simplification is possible:
∂c
∂ 2c
∂c
= Dh 2 − v
∂t
∂z
∂z

(18.37)

The above equations are usually referred to as advection-dispersion equations(ADEs), or convection-dispersion equations (CDEs).

814

D. Mallants et al.

18.3.2.2 Linear and Non-Linear Sorption
The advection-dispersion equation given by Eq. (18.34) contains three unknown
concentrations (those for the soil phase, liquid phase, and soil gas), while Eq. (18.35)
contains two unknowns. To be able to solve these equations, additional informationis needed that relates solid phase and liquid phase concentrations to each other.
A common approach 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 a linear isotherm, sometimes very appropriate at low
contaminant concentrations, given by:
s = Kd c(18.38)where Kd is the distribution or partitioning coefficient [L3 M−1 ]. Equation (18.38)
assumes reversible sorption (adsorption equals desorption). Substitution of thisequation into Eq. (18.36) leads to a constant value for the retardation factor, i.e.,

R=1+

ρKd
θ

(18.39)Whereas the use of a linear isotherm greatly simplifies the mathematical description of contaminant transport, sorption and exchange are generally non-linear andmost often depend also on the presence of competing species in the liquid phase.
Unlike linear adsorption, the contaminant retardation factor for non-linear adsorption is not constant, but changes as a function of concentration. Many models havebeen used in the past to describe non-linear sorption. Two widely used non-linear
sorption models are those by Freundlich (1909) and Langmuir (1918) given by:
s = Kf cβ(18.40)Kd c
1 + ηc(18.41)s=respectively, where Kf [L3 β M−β ] and β [−] are coefficients in the Freundlich
isotherm, and η [L3 M−1 ] is a coefficient in the Langmuir isotherm. Examples oflinear, Freundlich and Langmuir adsorption isotherms are given in Fig. 18.13.
18.3.2.3 Volatilization
Volatilization is increasingly recognized as an important process affecting the fateof many organic contaminants, including pesticides, fumigants, and explosives insoils (Jury et al. 1983a, 1984; Glotfelty and Schomburg 1989). While many organic
contaminants 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Leaching of Contaminants to Groundwater

Sorbed concentration [M M –1]

1.6

1.6

Freundlich model
β-parameter
0.5
0.75
1
1.5
2
3
5

1.2

815

Sorbed concentration [M M–1]

18

0.8

0.4

Langmuir model
η-parameter
0.5
0.75
1
1.5
2
3
5

1.2

0.8

0.4

0

0
0

0.4

0.8

1.2

1.6

0

–3

0.4

0.8

1.2

1.6

Dissolved concentration [M L–3]

Dissolved concentration [M L ]Fig. 18.13 Plots of the Freundlich adsorption isotherm given by Eq. (18.40), with Kf = 1 andβ defined in the caption (left), and the Langmuir adsorption isotherm given by Eq. (18.41), with
Kd = 1 and η defined in the caption (right)

physicochemical properties of the contaminant itself as well as such environmental variables as temperature and solar energy. Even though only a small fraction ofa pesticide may exist in the gas phase, air-phase diffusion rates can sometimes becomparable to diffusion in the liquid phase since gaseous diffusion coefficients areabout 4 orders of magnitude greater than liquid phase diffusion coefficients.
The general transport equation given by Eq. 18.34 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 (Eq. (18.38)) and KH , the latter appearinging = KH c(18.42)where KH is the dimensionless Henry constant [−].
Assuming linear partitioning, Eq. (18.34) for one-dimensional transport then has
the form∂∂c∂(ρb Kd + θ + aKH )c
=
θ Dh
∂t
∂z
∂z

+

∂
∂c
aDsg KH
∂z
∂z

−

∂(qc)
− φ (18.43)
∂x

or
∂
∂c
∂θ Rc
=
θ DE
∂t
∂z
∂z

−

∂(qc)
−φ
∂x

(18.44)

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

816

D. Mallants et al.

ρb Kd + aKH
θaDsg KH
DE = Dh +θ

R=1+(18.45)Jury et al. (1983b, 1984) provided for many organic contaminants their distribution coefficients Kd and Henry’s constants KH , from which calculated percent masspresent in each phase may be obtained.
18.3.2.4 Initial and Boundary Conditions
Initial Conditions
The governing equations for one-dimensional contaminant transport can be solved
analytically or numerically if the initial and boundary conditions for the problem being considered are specified. Initial conditions need to be provided for one
equilibrium phase concentration, i.e.,c(z, t) = ci (z, 0)(18.46)

where ci is the initial concentration [ML−3 ], as well as for all non-equilibriumphases if nonequilibrium transport is to be considered (see further).
Boundary Conditions
Complex interactions between the transport domain and its environment often must
be considered for the water flow part of the problem at hand since these interactions
determine the magnitude of water fluxes across the domain boundaries. By comparison, the contaminant transport part of most analytical and numerical models usuallyconsiders only three types of boundary conditions. When the concentration at theboundary is known, one can use a first-type (or Dirichlet type) boundary conditionof the form:c(z, t) = c0 (z, t)(18.47)

where c0 is a prescribed concentration [ML−3 ] at the top or bottom of the soil
profile. 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:− θD∂c+ qc = qc0∂x(18.48)in which q represents the boundary fluid flux [LT−1 ], and c0 is the concentrationof the incoming fluid [ML−3 ]. In some cases, for example when a boundary isimpermeable (q=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:

18

Leaching of Contaminants to Groundwater

θD

∂c
=0
∂x

817

(18.49)

Most applications require a Cauchy boundary condition rather than Dirichlet(or concentration) boundary condition. Since Cauchy boundary conditions define
the contaminant flux across a boundary, the contaminant flux entering the transport domain will be known exactly (as specified). This specified contaminant fluxis then in the transport domain divided into advective and dispersive components.
On the other hand, a Dirichlet boundary condition controls only the concentrationon the boundary, but not the contaminant flux into the domain, 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 early times, especially for relative
short transport domains (Van Genuchten and Parker 1984).

18.3.3 Nonequilibrium Transport
Because equilibrium contaminant transport models have frequently been unable
to describe experimental data accurately, a large number of diffusion-controlled
physical nonequilibrium and chemical-kinetic models have been developed andused to describe the transport of both non-adsorbing and adsorbing contaminants(Šimůnek and Van Genuchten 2008). Efforts to model nonequilibrium transport usually involve relatively simple first-order rate equations. Nonequilibrium models haveused the assumptions of two-region (dual-porosity) type transport involving contaminant exchange between mobile and immobile liquid phase transport regions, and/or
one-, two- or multi-site sorption formulations (e.g., Brusseau 1999; Nielsen et al.
1986). Models simulating the transport of particle-type contaminants, such as colloids, viruses, and bacteria, often also use first-order rate equations to describe suchprocesses as attachment, detachment, and straining. In many cases nonequilibrium
models have resulted in better descriptions of measured laboratory and field contaminant transport data, in part by providing additional degrees of freedom for fitting
measured concentration distributions.

18.3.3.1 Physical Nonequilibrium
Dual-Porosity and Mobile-Immobile Water Models
Two-region transport models (Fig. 18.14b and c) assume that the liquid phasecan be partitioned into distinct mobile (flowing) and immobile (stagnant) liquidpore regions, and that contaminant exchange between the two liquid regions canbe modeled as a first-order exchange process. Using the same notation as before,
the two-region contaminant transport model is given by (Toride et al. 1993; Van
Genuchten and Wagenet 1989):

818

D. Mallants et al.

a. Uniform Flowb. Mobile-immobile water

Water

Water

Solutec. Dual-porosity

Water
Imob.

SoluteImob.

θd. Dual-permeability

MobileSoluteMobileImob.

θ = θim + θmo

Mobile

Water
Slow

Fast

SoluteSlowθ = θim + θmFast

θ = θm + θf

Fig. 18.14 Conceptual models of water flow and contaminant transport (θ is the water content, θ moand θ im in (b) and (c) are water contents of the mobile and immobile flow regions, respectively,
with θ mo /(θ mo + θ im ) = φ m , and θ m and θ f in (d) are water contents of the matrix and macropore(fracture) regions, respectively) (after Šimůnek and Van Genuchten (2006))∂f ρsmo∂∂θmo cmo∂cmo+=θmo Dmo∂t∂t∂z∂z∂(1 − f )ρsim∂θim cim+= −φim + s
∂t
∂t

−

∂qcmo
− φmo − s∂z(18.50)for the mobile (macropores, subscript mo) and immobile (matrix, subscript im)
domains, respectively, where f is the dimensionless fraction of sorption sites incontact with the mobile water [−], φ mo and φ im are reactions in the mobile andimmobile domains [ML3 T−1 ], respectively, and Γ s is the contaminant transfer rate
between the two regions [ML3 T−1 ]. The same Eq. (18.50) can be used to describe
contaminant transport considering both the mobile-immobile and dual-porosity
models shown in Fig. 18.14b and c, respectively.
Dual-Permeability Model
One approach for implementing a dual-permeability formulation for contaminanttransport is to assume advection-dispersion type equations for transport in both thefracture and matrix regions as follows (Gerke and Van Genuchten 1993):∂ρsf∂θf cf
∂cf
∂
+
=
θf Df
∂t
∂t
∂z
∂z

−

∂qf cfs− φf −∂zw(18.51)∂θm cm
∂ρsm
∂
∂cm
+
=
θm Dm
∂t
∂t
∂z
∂z

−

∂qm cms− φm −∂z1−w(18.52)where the subscript f and m refer to the macroporous (fracture) and matrix poresystems, respectively; φ f and φ m represent sources or sinks in the macroporous andmatrix domains [ML3 T−1 ], respectively; and w is the ratio of the volume of themacropore (inter-aggregate) domain and that of the total soil system [-]. Equations(18.51) and (18.52) assume advective-dispersive type transport descriptions for both

18

Leaching of Contaminants to Groundwater819

the fractures and the matrix. Several authors simplified transport in the macropore
domain, for example by ignoring diffusion and dispersion in the macropores andconsidering only piston displacement there (Ahuja and Hebson 1992; Jarvis et al.
1994).
Mass Transfer
The transfer rate, Γ s , in Eq. (18.50) for contaminants between the mobile and immobile domains in the dual-porosity models can be given as the sum of diffusive andadvective fluxes as follows:s = αs (cmo − cim ) + w c ∗(18.53)where c∗ is equal to cmo for Γ w > 0 and cim for Γ w < 0, and α s is the first-ordercontaminant mass transfer coefficient [T−1 ]. Notice that the advection term of Eq.
18.53 is equal to zero for the mobile-immobile model (Fig. 18.14b) since the immobile water content in this model is assumed to be constant. However, Γ w may havea nonzero value in the dual-porosity model depicted in Fig. 18.14c.
The transfer rate, Γ s , in Eqs. (18.51) and (18.52) for contaminants between
the fracture and matrix regions is also usually given as the sum of diffusive andadvective fluxes as follows (e.g., Gerke and Van Genuchten 1996):
s = αs (1 − wm )(cf − cm ) + w c ∗(18.54)in which the mass transfer coefficient, Γ s [T−1 ], is of the form:
αs =

βg
Da
d2

(18.55)where β g is a dimensionless geometry-dependent coefficient, d is the characteristic
length of the matrix structure (L) (e.g., the radius of a spherical or solid cylindricalaggregate, or half the fracture spacing in the case of parallel rectangular voids), Da
is an effective diffusion coefficient [L2 T−1 ] representing the diffusion properties ofthe fracture-matrix interface.
18.3.3.2 Chemical Nonequilibrium
Kinetic Sorption Models
A substitute to expressing sorption as an instantaneous process using algebraic equations (e.g., Eqs. (18.38), (18.40) or (18.41)) is to quantify the reaction kinetics basedon ordinary differential equations. A popular and simple formulation of a chemicallycontrolled kinetic reaction arises when first-order linear kinetics is assumed:∂s= αk (Kd c − s)∂t

(18.56)

820

D. Mallants et al.

where α κ is a first-order kinetic rate coefficient [T−1 ]. Several other nonequilibrium
adsorption expressions have been used in the past (see Table 2 in Van Genuchtenand Šimůnek 1996), and are often referred to as one-site sorption models.
Because contaminant transport models assuming chemically controlled nonequilibrium (one-site sorption) often did not lead to significant improvements inpredictive capability when used to describe column transport experiments, theone-site first-order kinetic model was extended into a two-site sorption conceptthat divides the available sorption sites into two fractions (Selim et al. 1976; Van
Genuchten and Wagenet 1989). Conceptually, 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 linear sorption, the two-site transportmodel is given by (Van Genuchten and Wagenet 1989):
∂
∂(qc)
∂c
∂(f ρb Kd + θ )c
=
θ Dh
−
− φe
∂t
∂z
∂z
∂z
∂sk
= αk [(1 − f )Kd c − sk ] − φk∂t(18.57)

where f is the fraction of exchange sites assumed to be at equilibrium [−], φ e[ML3 T−1 ] and φ k [MM−1 T−1 ] are reactions in the equilibrium and nonequilibriumphases, respectively, and the subscript k refers to kinetic (type-2) sorption sites.
When 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). However, if f = 1, thetwo-site sorption model reduces to the equilibrium sorption model for which only
type-1 equilibrium sites are present.
18.3.3.3 Colloid-Facilitated Solute Transport
There is considerable evidence that many contaminants, including radionuclides(Noell et al. 1998; Von Gunten et al. 1988), pesticides (Kan and Tomson
1990; Lindqvist and Enfield 1992; Vinten et al. 1983), heavy metals (Grolimund
et al. 1996), viruses, pharmaceuticals (Thiele-Bruhn 2003; Tolls 2001), hormones(Hanselman et al. 2003), and other contaminants (Magee et al. 1991; Mansfeldt et al.
2004) in the subsurface are transported not only with moving water, but also sorbed
to mobile colloids. Because many colloids and microbes are negatively charged andthus electrostatically repelled by negatively-charged solid surfaces, the process ofanion exclusion may occur. As a result, contaminant transport could be slightly
enhanced relative to water 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). Sorption of contaminants onto mobile colloids can thus
significantly accelerate their transport relative to more standard advection-transport
circumstances.
Colloid-facilitated transport is a relatively complicated process that requires
knowledge of water flow, colloid transport, dissolved contaminant transport, andcolloid-contaminant interaction. This requires formulation of transport and/or mass

18

Leaching of Contaminants to Groundwater821

balance equations for both water flow and colloid transport, and in addition for
the total contaminant, for contaminants sorbed kinetically or instantaneously to thesolid phase, and for contaminants sorbed to mobile colloids, to colloids attachedto 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. The interested
reader will find additional documentation from de Jonge et al. (2004), Flury and
Qiu (2008), Massoudieh and Ginn (2009), Šimůnek et al. (2006a), and Totsche andKögel-Knabner (2004).

18.3.4 Stochastic Models
Several studies have demonstrated that solutions of classical contaminant transport models often fail to accurately describe transport processes in most field soils.
A major reason for the inability to capture field-scale contaminant transport is thefact that the subsurface environment is very heterogeneous. Heterogeneity occurs
at different levels of spatial and time scales (Wheatcraft and Cushman 1991),
ranging from microscopic scales involving time-dependent chemical sorption andprecipitation/dissolution reactions, to intermediate scales involving the preferential
movement of water and contaminants through macropores or fractures, and further
to large scales involving the spatial variability of soils across the landscape. Soil 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 tackled using stochastic approaches which incorporate certain assumptions about thetransport process in the heterogeneous system (e.g., Dagan 1989; Sposito and Barry
1987). In this Section we briefly review flow and transport parameter heterogeneityand discuss several stochastic transport approaches, notably those using stream tube
models and the transfer function approach.
18.3.4.1 Flow and Transport Parameter Heterogeneity
Spatial Variation of Physical and Chemical Properties
Hydrological and geological processes are known to vary in space (Delhomme
1979; Nielsen et al. 1973). This is due to the natural evolution of the subsurface
over geological time scales and anthropogenic activities in recent times. Knowledgeof spatial variability of physical and chemical properties is of utmost importancenot only from the point of view of characterisation of the heterogeneous subsurface, but also in the context of long term prediction of water flow and contaminanttransport in the subsurface. In this context, the last four decades has seen significant
number of studies carried out in order to comprehend the spatial variability of physical and chemical properties of soil with the help of both conventional statistics andgeostatistics.
From the vast body of literature, a wealth of data have been meticulously compiled by Jury (1985) and Thibault et al. (1990) on the spatial variability of physicaland chemical properties respectively; salient data from these compilations are presented in Tables 18.2 and 18.3. Note that the properties which are treated statistically

Sand
Clay

Static properties
Porosity
Porosity
Bulk density
Bulk density
Water content (0.1 bar, low suction)
Water content (0.1 bar, low suction)
Water content (15 bar, high suction)
Water content (15 bar, high suction)
Water content
Water content
Van Genuchten parameter, θ r
Van Genuchten parameter, θ s
Van Genuchten parameter, α
Van Genuchten parameter, nSoil pH
Soil pH
Kd
CEC
Langmuir parameter, Kf
Langmuir parameter, nSand, Loam, Clay and OC

Parameter
9.6
11.4
2.3
26.0
4.3
20.0
14
45
7.0
23.0
54.9–156.4
5.1–7.6
45.0–53.0
11.0–22.0
2.0
15.0
31.0
35.5
101.198.06
82.6, 37.2, 32.4
and 78.1
3.0
55.0

CV (%)

Loamy sand
Clay

Loamy sandClay loam
Sand
Sandy loam
Clay loam
Sandy loam
Clay loam
Sandy loam
Sandy loam
Sandy loam
Sandy loam
Sandy loam
Sandy loam
Sandy loam
Clay loam
Clay loam
Loamy sandLoamy clay
Clay loam
Clay loam
Clay loam

Soil texture

5
480

12
120
5
30
100
64
900
172
180
180
180
180
180
180
100
1040
36
93
93
93
93

Number of observations

Babalola (1978)
Nielsen et al. (1973)

Cassel (1983)
Russo and Bresler (1981)
Babalola (1978)
Courtin et al. (1983)
Gajem et al. (1981)
Gumaa (1978)
Gajem et al. (1981)
Cassel and Bauer (1975)
Mallants et al. (1996c)
Mallants et al. (1996c)
Mallants et al. (1996c)
Mallants et al. (1996c)
Mallants et al. (1996c)
Mallants et al. (1996c)
Gajem et al. (1981)
Cameron et al. (1971)
El Abd (1984)
Jacques et al. (1999)
Jacques et al. (1999)
Jacques et al. (1999)
Jacques et al. (1999)

References

Table 18.2 Replicate measurements of field transport parameters and their variability expressed by the coefficient of variation CV

822
D. Mallants et al.

48320322 (C1 horizon)
and 897 (C2horizon)284 (90%saturation),
412 (60%)
61
127
36
194

CV (%)

Loamy sandLoamy sand
Sandy
Loamy clay

Clay loam

Silty clay loam
Very coarse
Sandy loam

Soil texture

Partly extracted from Jury (1985); extreme values have been presented

Solute concentration
Solute concentration
Solute travel time
Solute travel time

Water diffusivity

Dynamic properties
Ks
Ks
Ks

Parameter

Table 18.2 (continued)

32
36
14 sites
20 plots

20

33
36
180

Number of observations

Wild and Babiker (1976)
Jury et al. (1983b)
Jury et al. (1982)
Biggar and Nielsen (1976)

Nielsen et al. (1973)

Willardson and Hurst (1965)
Willardson and Hurst (1965)
Mallants et al. (1996c)

References

18
Leaching of Contaminants to Groundwater
823

824

D. Mallants et al.

Table 18.3 Kd variability based on statistical analysis of literature data (excerpt from Thibault
et al. 1990)
Parameterσ∗

Soil texture

Range of Kd values – L Kg−1

Number of observations

Cd

1.6
1.5
0.9
2.9
2.1
–
1.4
2.3
–
2.4
2.6
1.4

Silt
Sand
Clay
Silt
Sand
Clay
Silt
Sand
Clay
Silt
Sand
Clay

7–962
2.7−625112−2,450
2.2−1,000
1.7−1,729
1,500 (exp(μ))
3,500–59,000
19–1,405
550 (exp(μ))
3.6–11,000
0.1–8,000
200–100,000

8
14
10
4
15
−
3
3
−
12
22
23

Cr

Pb

Zn

∗ Standard deviation of the logarithms of the observed values

in this section can be loosely categorized as static properties and dynamic properties. Static properties include porosity, bulk density, compositional properties suchas particle size fractions, water retention properties, soil pH, cation exchange capacity (CEC) and such distribution coefficients as the linear Kd , the Freundlich Kf andthe Langmuir Kd . Dynamic properties include the soil water diffusivity, hydraulicconductivity, contaminant breakthrough curves, dispersivity and contaminant concentrations (Jury 1985). The following sections summarize salient aspects of someof the static and dynamic properties from these compilations, supplemented withrecent data on spatial variability.
Static Properties
Jury’s (1985) compilations suggest that the range of the coefficient of variation(CV=standard deviation/mean) for static properties is surprisingly small given thevariety of soils and field sizes where the experiments were conducted. For instance,
the CV for porosity ranged between 9.6–11.4% and soil bulk density varied between
2.3–26%. In some studies, replicate measurements of the volumetric or gravimetricwater content in equilibrium with an imposed suction of 0.1 bar showed CV’s ranging from 4.3 to 20%; however, when the imposed suction was raised to 15 bars, theCV was higher and had a range of 14–45%. This is similar to observations reportedby Mallants et al. (1996c) who obtained a CV of 23% at 2.3 pF compared to 7% at
0 pF, although the trend reversed at higher pFs.
From Jury’s (1985) compilations, soil pH was generally found to be quite uniform with a maximum CV of 15%. For the Kd values they reported only one study
based on batch and column measurements, which showed a low CV of 31% on soil
samples taken across a sandy field of 0.64 ha.
Thibault et al. (1990) presented a comprehensive compilation of Kd values forvarious soil textures and elements. Table 18.3 only presents the values for some

18

Leaching of Contaminants to Groundwater825

heavy metals for three soil textures (the total number of elements reported byThibault et al. (1990) was 49). It is evident that the range of values encounteredfor various soil types varies significantly and over several orders of magnitude. In
general, Kd values for sand are lower than for silt or clay for the elements shown.
Although not shown in Table 18.3, Thibault et al. (1990) noted that elements suchas silver, nickel, americium and thallium showed higher mean Kd values for organicsoils than for clay.
In yet another study, Jacques et al. (1999) reported the spatial variability of theFreundlich Kf and n parameters for atrazine and their correlation with soil texturalvariables, cation exchange capacity and organic carbon content of a stagnic podzoluvisol based on 93 samples. The CV values (of non-transformed variables) for sand,loam, clay and organic carbon content were 82.6, 37.2, 32.4, and 78.1%, respectively. The CV values for CEC, Kf and n were reported to be 35.5, 101.2, and 8.06%
respectively. These properties therefore have CV values that vary from moderate tohigh, with Kf showing the highest variation.
Dynamic Properties
Jury (1985) noted that, in general, the CV values of dynamic properties are higher
than those of the static properties. For instance, based on thirteen studies reporting replicated measurements of the saturated Ks , the CV was found to be between
48 and 320%. Interestingly, both the lowest (48%) and the second highest (190%)
CV’s were measured on finer-textured soils. Thus, Jury (1982) commented that there
appears to be no apparent relationship between the variability of Ks and soil type,at least within the soil textural classes considered in their compilation. His findings
have been contradicted by others (e.g. Cosby et al. 1984).
As regards the contaminant transport properties, observations of contaminantconcentration variability have also been compiled by Jury (1985). The results basedon four field scale studies in which enough replicates were taken to allow a meaningful determination of sample concentration variance revealed a range of 61–127%(Table 18.2). He also reported the log variance of contaminant velocity as anothertransport property of interest. Based on six studies, a relatively low CV of 36% anda high CV of 194% were reported.
18.3.4.2 Stream Tube Models
The downward movement of contaminants 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” (Fig. 18.15),
while contaminant mixing between the stream tubes is assumed to be negligible. Transport in each tube may be described deterministically with the standard
advection-dispersion equation (ADE), or modifications thereof to include additional
geochemical and microbiological processes. Transport at the field scale is thenimplemented by considering the column parameters as realizations of a stochastic
process, having a random distribution (Toride et al. 1995).

826

D. Mallants et al.

Fig. 18.15 Schematic
illustration of the stream tube
model (left), (Toride et al.
1995)

The stream tube model was implemented into the CXTFIT 2.0 code (Toride
et al. 1995) for a variety of transport scenarios in which the pore water velocity
in combination with either the hydrodynamic dispersion coefficient, Dh , the distribution coefficient for linear adsorption, Kd , or the first-order rate coefficient fornonequilibrium adsorption, α k , are stochastic variables (Toride et al. 1995).

18.3.5 Multicomponent Reactive Solute Transport
The various mathematical descriptions of contaminant transport presented thus far
all considered contaminants that would move independently of other contaminantsin the subsurface. In reality, the transport of reactive contaminants is more often
than not affected by interactive physico-chemical and even biochemical processes.
Simulating these processes requires a more comprehensive approach that couples
the physical processes of water flow and advective-dispersive transport with a rangeof biogeochemical processes. The liquid phase is always a mixture of many ionswhich may be involved in mutually dependent chemical processes, such as complexation reactions (Lichtner 1996; Yeh and Tripathi 1990), cation exchange (White andZelazny 1986), precipitation-dissolution (Šimůnek and Valocchi 2002), sorptiondesorption, volatilization, redox reactions, and degradation, among other reactions(Šimůnek and Valocchi 2002). Transport and transformation of many 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 theelectron acceptor. By catalyzing such reactions, bacteria gain energy and organiccarbon to produce new biomass. These and related processes can be simulatedusing integrated reactive transport codes that couple the physical processes of waterflow and advective-dispersive contaminant transport with a range of biogeochemical
processes (Jacques et al. 2003; Šimůnek et al. 2006b).
Once the various chemical reactions are defined, the final system of governing
equations usually consists of several partial differential equations for contaminant18

Leaching of Contaminants to Groundwater827

transport (i.e., advection-dispersion equations for each component) plus a set of nonlinear algebraic and ordinary differential equations describing the equilibrium andkinetic reactions, respectively. Each contaminant and/or biological reaction must be
represented by corresponding algebraic or ordinary differential equations depending upon the rate of the reaction. Since the reaction of one species depends upon theconcentration of many other species, the final sets of equations are tightly coupled.
For complex geochemical systems, consisting of many components and multidimensional transport, numerical solution of these coupled equations is challenging(Šimůnek and Valocchi 2002). As an alternative, more general models have recently
been developed that more loosely couple transport and chemistry using a variety
of sequential iterative or non-iterative operator-splitting approaches (e.g., Bell andBinning 2004; Jacques and Šimůnek 2005; Jacques et al. 2006). Models based onthese various approaches are further discussed in Section 18.5.2.2.

18.3.6 Multiphase Flow and Transport
While the transport of contaminants in variably saturated media generally involves
two phases (i.e., the liquid phase and soil gas, with advection in the gaseous phaseoften 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 manyindustrial 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
dense nonaqueous phase liquids, or DNAPLs) or lighter than water (light nonaqueous phase liquids, or LNAPLs). Their fate and dynamics in the subsurface is affectedby 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 alsousing separate equations for air and water. This shows that the standard Richards
equation (Eq. (18.11)) is a simplification of a more complete multiphase (airwater) approach in that the air phase is assumed to have a negligible effect onvariably-saturated flow, and that the air pressure varies only little in space andtime. This assumption appears adequate for most variably-saturated flow problems.
Similar assumptions, however, are generally not possible when NAPLs are present.
Mathematical descriptions of multiphase flow and transport in general hence require
separate flow equations for each of the three fluid phases, mass transport equations
for all organic components (including those associated with the solid phase), andappropriate equations to account for interphase mass transfer processes. We referreaders to reviews by Abriola et al. (1999) and Rathfelder et al. (2000) for discussions of the complexities involved in modeling systems subject to multiphase flow,
multicomponent transport and interphase mass transfer. A useful overview of a variety of experimental approaches for measuring the physical and hydraulic propertiesof multi-fluid systems is given by Lenhard et al. (2002).

828

D. Mallants et al.

18.4 Analytical Models
18.4.1 Analytical Approaches
Many analytical solutions have been derived in the past of the contaminant transport
equations, and they are now widely used for analyzing contaminant transport during
steady-state flow (Šimůnek 2005). Although a large number of analytical solutions
also exist for the unsaturated flow equation, they generally can be applied only torelatively simple flow problems. The majority of applications for water flow in thevadose zone require a numerical solution of the Richards equation.
Analytical methods are representative of the classical mathematical approach forsolving differential equations to produce an exact solution for a particular problem. Analytical models usually lead to an explicit equation for the concentration (orthe pressure head, water content, or temperature) at a particular time and location.
One hence can evaluate the concentration directly without time stepping typicalof numerical methods. While exceptions exist (e.g., Liu et al. 2000), analyticalsolutions usually can be derived only for simplified transport systems involving
linearized governing equations, homogeneous soils, simplified geometries of thetransport 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 contaminant transport with nonlinear reactions,
are generally not available and/or cannot be derived, in which case numerical models
must be adopted (Šimůnek 2005).
Analytical solutions are usually obtained by applying various transformations(e.g., Laplace, Fourier or other transforms) to the governing equations, invoking aseparation of variables, and/or using the Green’s function approach (e.g., Leij et al.
2000).

18.4.2 Existing Models
18.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 transport in both finite and semi-infinite
domains, CFITIM additionally considers physical and chemical nonequilibriumtransport (i.e., the two-region mobile-immobile model for physical nonequilibrium
and the two-site sorption model for chemical nonequilibrium). CXTFIT 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:18

Leaching of Contaminants to Groundwater829• the conventional advection-dispersion equation, ADE;
• the chemical and physical nonequilibrium ADEs;• 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- andfirst-order decay/source terms.
18.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 includes
parameter estimation capabilities.
A large number of analytical models for one-, two-, and three-dimensional
contaminant transport problems were incorporated into the public domain software package STANMOD (STudio of ANalytical MODels) (Šimůnek et al.
1999a) (http://www.hydrus2d.com). This Windows-based computer software package includes not only programs for equilibrium advective-dispersive transport suchas the CFITM code of Van Genuchten (1980b) for one-dimensional transport and
3DADE (Leij and Bradford 1994) for three-dimensional problems, but also programs for more complex problems. For example, STANMOD also incorporates the
CFITIM (Van Genuchten 1981) and N3DADE (Leij and Toride 1997) programs fornonequilibrium transport (i.e., the two-region mobile-immobile model for physicalnonequilibrium and the two-site sorption model for chemical nonequilibrium) in
one and multiple dimensions, respectively. A more recent version of STANMOD
includes additionally the screening model of Jury et al. (1983a) for transport andvolatilization of soil-applied organic contaminants.

18.5 Numerical Models
18.5.1 Numerical Approaches
Although analytical and semi-analytical solutions are still popularly used for solving
many relatively simple problems, the ever-increasing power of personal computers
and the development of more accurate and stable numerical solution techniques have
led to the much wider use of numerical models over the past ten years. Numericalmethods in general are superior to analytical methods in terms of their ability tosolve much more realistic problems (Šimůnek 2005). They allow users to design
complicated geometries that reflect complex natural pedological and hydrologicalconditions, control parameters in space and time, prescribe more realistic initialand boundary conditions, and permit the implementation of nonlinear constitutive

830

D. Mallants et al.

relationships (Šimůnek and Van Genuchten 2006). Numerical methods subdivide
the time and spatial coordinates into smaller pieces, such as finite differences, finite
elements, and/or finite volumes, and reformulate the continuous form of governing
partial differential equations in terms of a system of algebraic equations. In orderto obtain solutions at prescribed times, numerical methods generally require intermediate simulations (time-stepping) between the initial condition and the points intime for which the solution is needed.
Reviews of the history of development of various numerical techniques used invadose zone flow and contaminant transport models are given by Van Genuchtenand Šimůnek (1996) and Šimůnek (2005).
18.5.1.1 Finite Differences
Finite difference methods are generally very intuitive and relatively easy to implement. Time and space are divided into small increments t and z (or x and z)(Fig. 18.16). Temporal and spatial derivatives in the governing equations are thenreplaced with finite differences (formally using Taylor series expansions). For example, the standard advection-dispersion equation for steady-state water flow (Wangand Anderson 1982) given by
∂J
∂ 2c
∂c
∂c
=−
=D 2 −v
∂t
∂z
∂z
∂z

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

ci

j

ti

zi+1

j

j

j

j

j

j

j
Ji+1/2 − Ji−1/2
c − 2ci + ci−1
c − ci−1− ci=−= D i+1− v i+1
2
t
z
2 z
( z)

ti+1
j
ci+1

zj+1

j+1
ci+1

xi-1

xi

ci–1, j+1

ci, j+1(18.59)xi+1
ci+1, j+1

Δzi

zi

cij

cij+1

j
i−1

j+1
i−1

z

Δzj

Δz

zj

c

c

ci+1, j

ci, j

Δz
Δzj–1

Δzi–1

z
zi–1

ci–1, jzj–1

ci–1, j–1

ci, j–1

ci+1, j−1

Δx

Δt

Δxi–1

t

x

Δxi

Fig. 18.16 Examples of the spatial and temporal difference discretization of a one-dimensional
problem (left), and the finite difference discretization of a two-dimensional domain (right) (afterŠimůnek and Van Genuchten (2006))

18

Leaching of Contaminants to Groundwater831

where subscripts refer to spatial discretization and superscript to temporal
discretization (e.g., j and j+1 are for the previous and actual time levels, respectively; see Fig. 18.16), t is the time step, and z is the spatial step (assumed tobe constant). Notice that this equation contains only one unknown variable (i.e., theconcentration ci j+1 at the new time level), which hence can be evaluated directly(explicitly) by solving the equation.
By comparison, a fully implicit (backward-in-time) finite difference scheme canbe written as follows
j+1

j+1

j+1

j+1

j+1

j+1

j+1

j
Ji+1/2 − Ji−1/2− 2ci + ci−1
− ci−1cc− ci− v i+1=−= D i+1
2
t
z
2 z
( z)
(18.60)and an implicit (weighted) finite difference scheme as:
j+1

ci

j+1
j
ci − ci

t

=D





j+1
j+1
j+1
j
j
j
ε ci+1 − 2ci + ci−1 + (1 − ε) ci+1 − 2ci + ci−1

−v

ε



j+1
j+1
ci+1 − ci−1



( z)2


j
j
+ (1 − ε) ci+1 − ci−1

(18.61)2 z

where ε is a temporal weighting coefficient. Different finite difference schemes
result depending upon the value of ε, i.e., an explicit scheme when ε = 0, a CrankNicholson time-centered scheme when ε = 0.5, and a fully implicit scheme when
ε=1.
18.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 ofthe finite elements is that they are much easier used to discretize complex twoand three-dimensional transport domains (Fig. 18.17). As an example, Fig. 18.17
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 (Šimůnek et al. 1999b). Notice that even though the figure on theright (Fig. 18.17) has an irregular soil surface, as well as a tile drain within thetransport domain, MeshGen2D could easily discretize/accommodate this transportdomain using an unstructured triangular finite element mesh.

18.5.2 Existing Models
18.5.2.1 Single-Species Solute Transport Models
A large number of numerical models are now available for evaluating variablysaturated water flow and contaminant transport processes in the subsurface. Someof these models are in the public domain, such as MACRO (Jarvis et al. 1994),
SWAP (van Dam et al. 1997), UNSATH (Fayer 2000), VS2DI (Healy 1990), and832

D. Mallants et al.

Fig. 18.17 Examples of triangular finite element grids for regular (left) and irregular (right) twodimensional transport domains. The problem on the left displays a 1-m wide and 2-m high soil
profile with a 0.01-cm wide macropore in the middle and a non-uniform steady-state water content
distribution. The project on the right shows a drain in between two furrows (one dry, one filled withwater) and the resulting steady-state pressure head distribution (groundwater table corresponds to
h = 0)HYDRUS-1D (Šimůnek and Van Genuchten 2008; Šimůnek et al. 1998a, 2005),
while others are in the commercial domain, such as HYDRUS-2D (Šimůnek et al.
1999b), HYDRUS (2D/3D) (Šimůnek et al. 2006c), and MODFLOW-SURFACT(HydroGeoLogic 1996). The 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 andinterpreting tabulated outputs, many others, especially those in the commercial
domain, are supported by sophisticated graphics-based interfaces that greatly simplify their use (Šimůnek et al. 1998a, 1999b). Several studies have recently reviewed
and compared various numerical models for vadose zone applications (e.g., MDH
Engineered Solutions Corp. 2003; Scanlon et al. 2002; Vanderborght et al. 2005;
Wilson et al. 1999). These studies typically compared the precision, speed, and easeof use of the codes involved.

18

Leaching of Contaminants to Groundwater

833

While 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), the more recent models can cope with much
more complex system-dependent boundary conditions evaluating surface flow andenergy 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(Šimůnek et al. 1998a, 2005), UNSATH (Fayer 2000), and COUP (Jansson andKarlberg 2001). Several models now also account for the extremely nonlinear processes associated with the freezing and thawing cycle (e.g., DAISY, SHAW, andCOUP).
Contaminant transport models have also become more sophisticated in terms ofthe type and complexity of processes that can be simulated. Transport models are
no longer being limited to contaminants 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 nonequilibriumtransport, volatilization, gas diffusion, colloid attachment/ detachment, decay chain
reactions, and many other processes (e.g., the HYDRUS-1D, –2D, and (2D/3D)
codes of Šimůnek et al. (1999b, 2005, 2006c), or MODFLOW-SURFACT of
HydroGeoLogic, Inc. (1996)). For example, the general formulation of the transport
equations in the HYDRUS codes permit simulations of non-adsorbing or linearlysorbing contaminants, in addition to a variety of other contaminants, such a viruses(Schijven and Šimůnek 2002), colloids (Bradford et al. 2002), cadmium (Seuntjens
et al. 2001), and hormones (Casey et al. 2003, 2004), or contaminants involvedin the sequential biodegradation of chlorinated aliphatic hydrocarbons (Casey andŠimůnek 2001; Schaerlaekens et al. 1999).
Much effort has been directed also toward improving models for purposes of simulating nonequilibrium and/or preferential flow. Examples are the TOUGH codes(Pruess 1991, 2004), MACRO (Jarvis et al. 1994), and HYDRUS-1D (Šimůnek
and Van Genuchten 2008). These models typically assume the presence of dualporosity and dual-permeability regions, with different fluxes possible in the tworegions. Example applications of these dual-porosity and dual-permeability models
are given by Mallants et al. (1997), Šimůnek et al. (2001), Haws et al. (2005), Köhne
et al. (2004, 2006, 2009a, b), and Pot et al. (2005), among many others.
As an example of available vadose zone flow and transport models, we briefly
discuss here the HYDRUS software packages of Šimůnek et al. (1999b, 2005, 2008).
The HYDRUS Software Packages
HYDRUS-1D (Šimůnek et al. 2005), HYDRUS-2D (Šimůnek et al. 1999b), andHYDRUS (2D/3D) (Šimůnek et al. 2006c) are software packages (http://www.pcprogress.com/en/Default.aspx) that simulate the one- and two-dimensional movement of water, heat, and multiple contaminants in variably saturated porous
media, respectively. Both programs use finite elements to numerically solve
the Richards equation for saturated-unsaturated water flow and Fickian-based

834

D. Mallants et al.

advection-dispersion equations for both heat and contaminant transport. The
unsaturated soil hydraulic properties 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.
The HYDRUS-1D software package additionally includes modules for simulating carbon dioxide and major ion contaminant movement (Šimůnek et al. 1996;Šimůnek and Suarez 1993). Also included is a small catalog of unsaturated soil
hydraulic properties (Carsel and Parish 1988), as well as pedotransfer functions
based on neural network predictions (Schaap et al. 2001).
18.5.2.2 Biogeochemical Transport Models
Significant efforts have been made also in coupling physical flow and transportmodels with biogeochemical models to simulate increasingly more complex reactions, such as surface complexation, precipitation/dissolution, cation exchange,
and/or (micro)biological reactions. Reviews of the development of hydrogeochemical transport models involving reactive multiple components are given by Mangold
and Tsang (1991), Lichtner (1996), Steefel and MacQuarrie (1996), Šimůnek andValocchi (2002), and Bell and Binning (2004). Most modeling efforts involving multicomponent transport have thus far focused on the saturated zone, where changes in
the flow velocity, temperature and pH are often much more gradual and hence lessimportant than in the unsaturated zone. Consequently, most multicomponent transport models assumed one- or two-dimensional steady-state saturated water flow witha fixed value of the flow velocity, temperature and pH. Several multicomponent
transport models have been published also for variably-saturated flow problems.
These include DYNAMIX (Liu and Narasimhan 1989), HYDROGEOCHEM (Yehand Tripathi 1990), TOUGH-REACT (Pruess 1991), UNSATCHEM (Šimůnek
and Suarez 1994; Šimůnek 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 and Šimůnek 2005; Jacques et al. 2002, 2008a,
2008b), and HYDRUS-1D (Šimůnek et al. 2005).
Geochemical models can be divided into two major groups: those with specific chemistry and those characterized by more general chemistry (Šimůnek andValocchi 2002). Models with specific chemistry are limited in the number of species
they can handle, while their application is restricted to problems having a prescribedchemical system. They are, however, much easier to use and computationallycan 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 (Šimůnek and Suarez 1994;
Šimůnek et al. 1996), and HYDRUS-1D (Šimůnek 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 thesemodels are also much wider.

18

Leaching of Contaminants to Groundwater

835

HP1
HYDRUS-1D was coupled with the PHREEQC geochemical code (Parkhurst andAppelo 1999) to create a new comprehensive simulation tool, HP1 (acronym for
HYDRUS1D-PHREEQC) (Jacques and Šimůnek 2005; Jacques et al. 2003, 2008a,
2008b). The combined code contains modules simulating (1) transient water flow invariably-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 typeequations for heat and contaminant 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.
Jacques et al. (2003, 2008a, 2008b) and Jacques and Šimůnek (2005) demonstrated the versatility of the HP1 model on several examples such as:• the transport of heavy metals (Zn2+ , Pb2+ , and Cd2+ ) subject to multiple cationexchange reactions;• transport with mineral dissolution of amorphous SiO2 and gibbsite (Al(OH)3 );• heavy metal transport in a medium with a pH-dependent cation exchange
complex;• 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);• long-term transient flow and transport of major cations (Na+ , K+ , Ca2+ , andMg2+ ) and heavy metals (Cd2+ , Zn2+ , and Pb2+ ) in a soil profile• cadmium leaching in acid sandy soils;• radionuclide transport following phosphorus fertilization (U and its aqueous
complexes);• the fate and subsurface transport of explosives (TNT and its daughter products
2ADNT, 4ADNT, and TAT) (Šimůnek et al. 2006b).
To illustrate the capabilities of HP1 and the importance of coupling in an integrated manner all relevant physical and contaminant transport processes whenanalyzing the fate and transport of contaminants in variably saturated field soils, thenext section discusses the leaching of cadmium from a podsol soil under transientboundary conditions.
Applications to Unsaturated Flow and Geochemical Transport ModelingIn this example (Jacques et al. 2008a) we discuss a hypothetical HP1 application involving the transport of major cations and heavy metals in a soil duringtransient flow over a period of 30 years. Results will show that variations in water

836

D. Mallants et al.

contents and water fluxes can significantly influence the speciation, and thus themobility and availability, of elements. Decreasing water contents near the soil surface, furthermore, may lower the pH of the soil liquid phase and produce new cationexchange equilibrium conditions. The upward transport of Cl during summer due toincreased evapotranspiration, and subsequent accumulation of Cl near the soil surface, can cause an increase in the total aqueous cadmium concentration because of
the formation of cadmium-Cl complexes.
The HP1 multi-component transport simulator is used in this section to simulate
the leaching of cadmium and zinc in a dry Spodosol in a sandy region of Northern
Belgium using in-situ measured cadmium and zinc (and additional elements) concentration profiles. Soils in the region were contaminated by atmospheric depositionof cadmium and zinc from non-ferrous industry (Seuntjens 2000). Water flow andthe transport of major cations (Na, K, Mg, and Ca), heavy metals (cadmium andzinc), anions (Cl, Br) and Al were simulated for a 1-m deep multi-layered podsol
soil profile subject to atmospheric boundary conditions for a period of 30 years.
The main focus is on how processes affecting water contents and water fluxes also
influence the geochemical conditions in the soil. Specifically, the effect of cycles of
evaporation and infiltration on pH and cadmium speciation will be discussed.
Interactions between major cations and heavy metals with the soil solid phasewere simulated by means of cation exchange processes assuming local equilibriumon a single type of exchange sites. Voegelin (2001) earlier showed that this approachadequately describes various features of cadmium transport experiments. Althoughsorption on specific sites with a high affinity for cadmium may also occur in soils(Selim et al. 1992), this type of binding is unlikely in acid sandy soils (Voegelin
et al. 2001). Exchange parameters were calibrated using concentrations measured
in drainage water from a steady-state flow experiment on large, undisturbed soil
lysimeters (1 m long and 0.8 m diameter; Seuntjens et al. 2001). The initial composition of the cation exchange site was measured for each soil horizon (Jacques et al.
2008a). While the same log(K) parameters were used for all soil horizons, the sizeof the cation exchange complex was assumed to vary between horizons.
Daily values of precipitation (P) and potential evapotranspiration (Ep) during
the 30-year simulation period were imposed as climatic boundary condition, fromwhich the daily actual evaporation rate was calculated. A subset of the P-Ep data isshown in Fig. 18.18.
Figure 18.19 shows time series between 1972 and 1982 for the water content,
pH and total Cl and Cd concentrations in the liquid phase at two depths. The results
illustrate the relation between water flow and geochemical conditions in the soil.
The alternation between precipitation (wet conditions) and evaporation (dry conditions) as dictated by the atmospheric conditions clearly affected the dynamics of
the water content, with upward water flow during dry periods. The flow dynamics
in turn significantly influenced the geochemistry near the soil surface. As illustrated, the most mobile elements (anions such as Cl− and monovalent cations suchas Na+ ) move upwards during the evaporation periods, thus causing these ions toaccumulate near the soil surface. The decrease in water content near the soil surface
due to evaporation resulted in higher concentrations and a lower pH. PHREEQC

18

Leaching of Contaminants to Groundwater

837

6
5

P - Ep(cm day–1)

4
3
2
1
0
–1
1972

1973

1974

1975

1976

1977

1978

1979

1980

1981

1982

Year

Fig. 18.18 Potential precipitation surplus during years 1972–1982

calculations showed that removing only water from aqueous solutions results in a
pH decrease (for a similar geochemical system, i.e., only cation exchange and equilibrium with gibbsite). Another factor intensifying the pH decrease was the greater
mobility of anions (Cl− and Br− ) compared to cations. Thus, the physical factorsof having a decreasing water content and upward flow of water and contaminantscaused the pH to decrease when P − Ep is negative. It should be noted that the pH
in reality is also affected by other geochemical or biological processes not includedin the invoked conceptual model. For example, soil carbon dioxide concentrations
that usually change in response to biological activity and moisture status of the soil
can also affect soil pH (Šimůnek and Suarez 1993).
Although upward flow during the summer had almost no effect on the total
amount of heavy metals in the surface horizon due to the low mobility of theseelements, the aqueous concentrations of the metals did vary significantly during
the season. Several factors contributed to this. First, because of lower water contents, the concentration of all aqueous species increased during summer periods.
Changes in aqueous concentrations in turn caused changes in the cation exchange
equilibrium, thereby promoting monovalent cations to sorb onto the cation exchange
complex and bivalent cations to desorb into solution. This explains also the difference between TNa,aq and TCa,aq during the summer near the soil surface. The
aqueous concentration of Na was controlled more by the cation exchange complex
than that of Ca due to preferred adsorption of Na during dry soil conditions. This
process is further amplified by the increased supply of monovalent cations due toupward flow of water during summer, leading to relatively more sorption of themonovalent cations and higher concentrations of divalent cations and heavy metals
in the liquid phase. The complexation of cadmium with Cl contributed also to the838

D. Mallants et al.
0.3

5 cm
35 cm

0.25

θ

0.2
0.15
0.1
0.05

pH

4.8

4.6

TCl,aq (mol/1000 cm³)

4.4
4x10–5
3x10–5
2x10–5
1x10–5

TCd,aq (mol/1000 cm³)0x100

1.0x10–8

5.0x10–9

0.0x100
1972

1974

1976

1978

1980

1982

Year

Fig. 18.19 Time series of water content (θ), pH, Tcl,aq and TCd,aq (mol/1000 cm3 soil) between
1972 and 1982 at 5 cm (A-horizon) and 35 cm (Bh1-horizon)

18

Leaching of Contaminants to Groundwater839

increase in TCd,aq . At a depth of 3 cm, about 0.65 and up to 3.5 percent of the aqueous cadmium was in the form of CdCl+ during winter and summer, respectively.
Other cadmium complexes of the form CdCln (2−n) were present in much smaller
concentrations. Similar to the pH changes, changing water contents and upward
fluxes both affected the amount of cadmium in the liquid phase.
Our results above indicate that atmospheric boundary conditions can have a significant effect on the amount and transport of cadmium in a soil profile, and on its
bioavailability since uptake processes by plants and soil micro-organisms are often
concentration-dependent. Passive root uptake of contaminants together with water
increases with increasing contaminant concentrations. Similarly, active uptake as
described with Monod or Michealis-Menten kinetics will increase with increasing contaminant concentrations. Moreover, the high heavy metal concentrations
occurred during the summer months with the highest (micro)biological activity. In
addition, cadmium speciation may also play a role in uptake. For example, Smolders
and McLaughlin (1996) observed more cadmium uptake by chard (Beta vulgaris var.
cicla) when Cl concentrations increased while the Cd2+ activity was kept constant,likely due to phytoavailability of CdCl+ and other CdCln 2–n -species.
The example illustrates that simulators such as HP1 are potentially attractive
tools for studying reactive transport processes in the vadose zone during transientvariably-saturated flow. Geochemical conditions are an important factor since they
determine the speciation (both in the liquid phase and on the solid phase) of theelements involved, and thus their mobility and bioavailability. Small variations inprevailing geochemical conditions may alter significantly the speciation and mobility of heavy metals or other constituents. Soil systems, moreover, are subject tolarge transient variations since they are open to the atmosphere. Hence, changes
in the composition of rain water or atmospheric deposition, among other external
factors, can materially alter the geochemical conditions in a soil profile.
For the time being, only a limited number of test cases and experimental data
exist for the unsaturated zone to assess the full capabilities of a reactive transport
code (see also Davis et al. 2004). One elaborate recent application is a study by
Gonçalves et al. (2006), who successfully applied the chemistry-specific major ion
geochemistry module of HYDRUS-1D (version 3.0, Šimůnek et al. 2005) to a 4-year
experimental data set involving water flow and contaminant transport in lysimetersirrigated with waters of different quality and subjected to atmospheric conditions.
Additional studies of this type should provide more credibility to the use of coupledhydrogeochemical models for addressing flow and reactive transport problems inthe vadose zone.

18.6 Concluding Remarks
This chapter demonstrates the abundance of models and modeling approaches thatare currently available for simulating variably-saturated water flow and contaminanttransport at various levels of approximation and for different applications. Models
range from relatively simple analytical approaches for analyzing contaminant

840

D. Mallants et al.

transport problems during one-dimensional steady-state flow, to sophisticatednumerical models for addressing multi-dimensional variably-saturated flow andcontaminant transport problems at the field scale.
One may expect that unsaturated zone flow and transport models will be used
increasingly as tools for higher tier Groundwater-related Risk Assessments (seeChapter 17 by Swartjes and Grima, this book). Moreover, these models may support the development of cost-effective, yet technically sound strategies for resource
management, contamination remediation, and/or 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. Models can be helpful
tools also for designing, testing and implementing soil, water and crop management
practices that minimize soil and water contamination. Models are equally neededfor designing or remediating industrial waste disposal sites and landfills, for predicting contaminant transport from mining wastes, or for long-term stewardship ofnuclear waste repositories. A major challenge is to make the models as realistic as
possible for such complex applications. Yet another challenge is to obtain input dataappropriate for the time and spatial scale under consideration, accounting for spatial
variability and possibly even time-dependent parameters. Several approaches werediscussed that permit generation of input data from other more basic (soil) data, suchas pedotransfer functions. These efforts must continue to also include (bio)chemicalparameters.
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 continue 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 transportmodels 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 coupledphysico-chemical systems to account for possible changes in the water retention
and hydraulic conductivity functions. Better integration is also needed between
variably-saturated zone and existing larger-scale surface numerical models, which
in turn requires further research on such issues as spatial and temporal scaling ofhydrological, chemical and biological processes and properties, linking constitutive(soil hydraulic) relationships to measurements scales, preferential flow, and issuesof parameter and model uncertainty.
Many scientific questions related to colloid and colloid-facilitated transport arealso still largely unresolved. This is an area of research where our understanding lags
far behind current numerical capabilities. Much work is needed to better understandthe 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),

18

Leaching of Contaminants to Groundwater841

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 conditionsof transient flow.
Also, to the best of our knowledge, no models are currentely available that consider all of the many 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 andinorganic colloids and/or organic complexes, and solid and air phases of a soil,
including widely varying rates of these various reactions. Further integration of thedifferent types of numerical models is needed to address practical problems of contaminant transport (trace elements, radionuclides, organic contaminants) in complex
vadose zone environments.

References
Abriola LM, Pennell KD, Weber WJ Jr, Lang JR, Wilkins MD (1999) Persistence and interphasemass transfer of organic contaminants in the unsaturated zone: experimental observations andmathematical modeling. In: Parlange JY, Hopmans JW (eds) Vadose zone hydrology: cutting
across disciplines. Oxford University, New York, pp 210–234
Ahuja LR, Hebson C (1992) Root zone water quality model. GPSR Tech Rep No 2, USDA, ARS,
Fort Collins, CO
Allen RG (2001) REF-ET: reference evapotranspiration calculation software for FAO and ASCE
standardized equations, version 2.0 for windows. University of Idaho, USA
Allen RG, Pereira LS, Raes D, Smith M (1998) Crop evapotranspiration: guidelines for computing crop water requirements, Irrigation and Drainage Paper, no. 56, Food and Agriculture
Organization of the United Nations, RomeAnderson MP (1984) Movement of contaminants in groundwater: groundwater transport-advection
and dispersion. In: Groundwater contamination. Studies in geophysics. National Academy,
Washington, DC, pp 37–45
Babalola O (1978) Spatial variability of soil water properties in tropical soils of Nigeria. Soil Sci
126:269–279
Bear J (1972) Dynamics of fluids in porous media. Elsevier Science, NY
Bell LSJ, Binning PJ (2004) A split operator approach to reactive transport with the forward particle tracking Eulerian-Lagrangian localized adjoint method. Adv Water Resour 27:
323–334
Biggar JW, Nielsen DR (1967) Miscible displacement and leaching phenomena. Agronomy
11:254–274
Biggar JW, Nielsen DR (1976) Spatial variability of the leaching characteristics of a field soil.
Water Resour Res 12:78–84
Bouma J (1989) Using soil survey data for quantitative land evaluation. Adv Soil Sci 9:
177–213
Bradford SA, Šimůnek J, Bettahar M, Van Genuchten MTh, Yates SR (2003) Modeling colloid attachment, straining, and exclusion in saturated porous media. Environ Sci Technol 37:
2242–2250
Bradford SA, Yates SR, Bettahar M, Šimůnek J (2002) Physical factors affecting the transport and fate of colloids in saturated porous media. Water Resour Res 38(12):1327.
doi:10.1029/2002WR001340

842

D. Mallants et al.

Brooks RH, Corey AT (1964) Hydraulic properties of porous media. Hydrol Paper No. 3, Colorado
State University, Fort Collins, CO
Brusseau ML (1999) Nonideal transport of reactive solutes in porous media; Cutting across history
and disciplines. In: Parlange MB, Hopmans JW (eds) Vadose zone hydrology: cutting acrossdisciplines. Oxford University, New York, NY, pp 130–154
Burdine NT (1953) Relative permeability calculations from pore-size distribution data. Petr Trans
Am Inst Mining Metall Eng 198:71–77
Cameron DR, Nyborg M, Toogood JA, Laverty DH (1971) Accuracy of field sampling for soil
tests. Can J Soil Sci 51:165–175
Carsel RF, Parrish RS (1988) Developing joint probability distributions of soil water retention
characteristics. Water Resour Res 24:755–769
Casey FXM, Šimůnek J (2001) Inverse analyses of the transport of chlorinated hydrocarbons
subject to sequential transformation reactions. J Environ Qual 30(4):1354–1360
Casey FXM, Larsen GL, Hakk H, Šimůnek J (2003) Fate and transport of 17β-Estradiol in soilwater systems. Environ Sci Technol 37(11):2400–2409
Casey FXM, Larsen GL, Hakk H, Šimůnek J (2004) Fate and transport of testosterone in
agriculturally significant soils. Environ Sci Technol 38(3):790–798
Cassel DK (1983) Spatial and temporal variability of soil physical properties following tillage of
Norfolk loamy sand. Soil Sci Soc Am J 47:196–201
Cassel DK, Bauer, A (1975) Spatial variability in soils below depth of tillage. Bulk density andfifteen atmosphere percentage. Soil Sci Soc Am J 39:247–250
Cosby BJ, Hornberger GM, Clapp RB, Ginn TR (1984) A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils. Water Resour Res 6:
682–690
Courtin P, Feller MC, Klinka K (1983) Lateral variability in some properties of disturbed forest
soils in SW British Columbia. Can J Soil Sci 63:529–539
Dagan G (1989) Flow and transport in porous formations. Springer, Berlin
Dane JH, Hopmans JW (2002) Water retention and storage. In: Dane JH, Topp GC (eds)
Methods of soil analysis, part 1, Physical methods, chapter 3.3, 3rd edn. SSSA, Madison, WI,
pp 671–720
Dane JH, Topp GC (2002) Methods of soil analysis, part 1, physical methods, 3rd edn. SSSA,
Madison, WI
Davis JA, Yabusaki SB, Steefel CI, Zachara JM, Curtis GP, Redden GB, Criscenti LJ,
Honeyman BD (2004) Assessing conceptual models for subsurface reactive transport ofinorganic contaminants. EOS 85:449, 455
de Jonge LW, Kjaergaard C, Moldrup P (2004) Colloids and colloid-facilitated transport ofcontaminants in soils: an introduction. Vadose Zone J 3:321–325
Delhomme JP (1979) Spatial variability and uncertainty in groundwater flow parameters: ageostatistical approach. Water Resour Res 15:269–280
Durner W (1994) Hydraulic conductivity estimation for soils with heterogeneous pore structure,
Water Resour Res 32(9):211–183
El Abd H (1984) Spatial variability of the pesticide distribution coefficient. PhD thesis, Universityof California, Riverside, 100 pp
Fayer MJ (2000) UNSAT-H Version 3.0: unsaturated soil water and heat flow model. Theory, user
manual, and examples. Pacific Northwest National Laboratory 13249
Feddes RA, Kowalik PJ, Zaradny H (1978) Simulation of field water use and crop yield. Wiley,
New York, NY
Fetter CW (1999) Contaminant hydrology. Upper Saddle River, Prentice Hall, 458 pp
Flerchinger GN, Hanson CL, Wight JR (1996) Modeling evapotranspiration and surface energy
budgets across a watershed. Water Resour Res 32:2539–2548
Flury M, Fluhler H, Jury W, Leuenberger J (1994) Susceptibility of soils to preferential flow ofwater: a field study. Water Resour Res 30:1945–1954

18

Leaching of Contaminants to Groundwater

843

Flury M, Qiu H (2008) Modeling colloid-facilitated contaminant transport in the Vadose zone.
Vadose Zone J 7:682–697
Freundlich H (1909) Kapillarchemie; eine Darstellung der Chemie der Kolloide und verwandter
Gebiete, Akademische Verlagsgellschaft, Leipzig, Germany
Gajem YM, Warrick AW, Myers DE (1981) Spatial dependence of physical properties of a typic
torrifluvent soil. Soil Sci Soc Am J 45:709–715
Gelhar LW, Mantaglou A, Welty C, Rehfeldt KR (1985) A review of field scale physical solutetransport processes in unsaturated and saturated porous media. EPRI Topicl Report EA-4190,
Electric Power Research Institute, Palo Alto, CA
Gerke HH, Van Genuchten MTh (1993) A dual-porosity model for simulating the preferential
movement of water and solutes in structured porous media. Water Resour Res 29:305–319
Gerke HH, Van Genuchten MTh (1996) Macroscopic representation of structural geometry for
simulating water and solute movement in dual-porosity media. Adv Water Resour 19:343–357
Glotfelty DE, Schomburg CJ (1989) Volatilization of pesticides. In: Sawhney BL, Brown K (eds)
Reactions and movement of organic chemicals in soils, special publ. no. 18. Soil Science
Society of America, Madison, WI, pp 181–208
Gonçalves MC, Šimůnek J., Ramos TB, Martins JC, Neves MJ, Pires FP (2006) Multicomponent
solute transport in soil lysimeters irrigated with waters of different quality. Water Resour Res
42:W08401. doi:10.1029/2005WR004802
Grolimund D, Borkovec M, Barmettler K, Sticher H (1996) Colloid-facilitated transport of strongly
sorbing contaminants in natural porous media; a laboratory column study. Environ Sci Technol
30:3118–3123
Gumaa GA (1978) Spatial variability of in situ available water. PhD Dissertation, University of
Arizona, Tucson, 140 pp
Gwo JP, Jardine PM, Wilson GV, Yeh GT (1995) A multiple-pore-region concept to modeling masstransfer in subsurface media. J Hydrol 164:217–237
Hanselman TA, Graetz DA, Wilkie AC (2003) Manure-borne estrogens as potential environmental
contaminants: a review. Environ Sci Technol 37:5471–5478
Hansen S, Jensen HE, Nielsen NE, Svendsen H (1990) DAISY: soil plant atmosphere systemmodel. NPO report no. A 10. The National Agency for Environmental Protection, Copenhagen,272 ppHaws NW, Rao PSC, Šimůnek J (2005) Single-porosity and dual-porosity modeling of waterflow and transport in subsurface-drained fields using effective field-scale parameters. J Hydrol
313(3–4):257–273
Healy RW (1990) Simulation of solute transport in variably saturated porous media with supplemental information on modifications to the U.S. Geological Survey’s computer program
VS2DUS. Geological Survey, Water-Resources Investigation Report 90–4025:125
Hendrickx JMH, Dekker LW, Boersma OH (1993) Unstable wetting fronts in water repellent fieldsoils. J Environ Qual 18:109–118
Hendrickx JMH, Flury M (2001) Uniform and preferential flow, mechanisms in the Vadose zone.
In: Conceptual models of flow and transport in the fractured Vadose zone. National ResearchCouncil, National Academy, Washington, DC, pp 149–187
HydroGeoLogic, Inc (1996) MODFLOW-SURFACT, a comprehensive MODFLOW-based flow
and transport simulator, version 2.1, HydroGeoLogic, Inc.
Jacques D, Mouvet JC, Mohanty BP,Vereecken H, Feyen J (1999) Spatial variability of atrazine
sorption parameters and other soil properties in a podzoluvisol. J Contamin Hydrol 36:
31–52
Jacques D, Šimůnek J (2005) Multicomponent – variable saturated transport model, description,
manual, verification and examples, waste and disposal, SCK•CEN, BLG-998, Mol, Belgium,
79 pp
Jacques D, Šimůnek J, Mallants D, Van Genuchten MTh (2002) Multicomponent transport model
for variably-saturated porous media: application to the transport of heavy metals in soils. In:
Hassanizadeh SM, Schotting RJ, Gray WG, Pinder GF (eds) Computational methods in water

844

D. Mallants et al.

resources, XIVth international conference on computational methods in water resources, June
23–28, 2002, Delft, The Netherlands. Elsevier, Developments in Water Science, vol 47, pp
555–562
Jacques D, Šimůnek J, Mallants D, Van Genuchten MTh (2003) The HYDRUS-PHREEQC multicomponent transport model for variably-saturated porous media: code verification and application, MODFLOW and More 2003: Understanding through Modeling, Conference Proceedings,
Poeter E, Zheng Ch, Hill M, Doherty J (eds) International Ground Water Modeling Center,
Colorado School of Mines, pp 23–27
Jacques D, Šimůnek J, Mallants D, Van Genuchten MTh (2006) Operator-splitting errors in coupled reactive transport codes for transient variably saturated flow and contaminant transport inlayered soil profiles. J Contamin Hydrol 88:197–218
Jacques D, Šimůnek J, Mallants D, Van Genuchten MTh (2008a) Modelling coupled water flow,
solute transport and geochemical reactions affecting heavy metal migration in a podzol soil.
Geoderma 145:449–491
Jacques D, Šimůnek J, Mallants D, Van Genuchten MTh (2008b) Modeling coupled hydrological and chemical processes in the Vadose zone: a case study of long-term uranium transport
following mineral fertilization.Vadose Zone J 7(2):698–711
Jansson P-E, Karlberg L (2001) Coupled heat and mass transfer model for soil-plant-atmospheresystems. Royal institute of technology, Department of Civil and Environmental Engineering,
Stockholm, 325 pp
Jarvis NJ (1994) The MACRO model (version 3.1), Technical description and sample simulations.
Reports and dissertations 19. Department of Soil Science, Swedish University Agricutural
Science, Uppsala, Sweden, 51 pp
Jarvis NJ (1998) Modeling the impact of preferential flow on nonpoint source pollution. In: Selim
HM, Ma L (eds) Physical nonequilibrium in soils: modeling and application. Ann Arbor,
Chelsea, MI, pp 195–181
Jensen ME, Burman RD, Allen RG (1990) Evapotranspiration and irrigation water requirement.
ASCE manuals and reports on enigineering practice, vol 70, ASCE, New YorkJin Y, Jury WA (1996) Characterizing the dependence of gas diffusion coefficient on soil properties.
Soil Sci Soc Am J 60:66–71
Jury WA (1982) Simulation of solute transport using a transfer function model. Water Resour Res
18:363–368
Jury WA (1985) Spatial variability of soil physical parameters in solute migration: a critical
literature review. EPRI EA-4228, Interim report, EPRI
Jury WA, El Abd H, Collins TM (1983b) Field scale transport of nonadsorbing and adsorbing
chemicals applied to the soil surface. NWWA symposium, Las Vegas. National Water Well
Assoc, Worthington, OH, pp 203–222
Jury WA, Gardner WA, Gardner WH (1991) Soil physics, 5th edn. Wiley, New York, NY
Jury WA, Horton R (2004) Soil physics, 6th edn. Wiley, New York, NY
Jury WA, Spencer WF, Farmer WJ (1983a) Behavior assessment model for trace organics in soil:
I. Description of model. J Environ Qual 12:558–564
Jury WA, Spencer WF, Farmer WJ (1984) Behavior assessment model for trace organics in soil: I.
Chemical classification and parameter sensitivity J Environ Qual 13:567–572
Jury WA, Stolzy LH, Shouse P (1982) A field test of the transfer function model for predicting
solute transport. Water Resour Res 18:369–375
Kan AT, Tomson MB (1990) Ground water transport of hydrophobic organic compounds in thepresence of dissolved organic matter. Environ Toxicol Chem 9:253–263
Köhne S, Lennartz B, Köhne JM, Šimůnek J (2006) Bromide transport at a tile-drained field site:experiment, one- and two-dimensional equilibrium and non-equilibrium numerical modeling.
J Hydrol 321(1–4):390–408
Köhne JM, Köhne S, Šimůnek J (2009a) A review of model applications for structured soils: a)water flow and tracer transport. J Contamin Hydrol, Special issue “flow domains” 104(1–4):
4–35. doi:10.1016/j.jconhyd.2008.10.002

18

Leaching of Contaminants to Groundwater

845

Köhne JM, Köhne S, Šimůnek J (2009b) A review of model applications for structured soils:b) pesticide transport. J Contamin Hydrol, Special issue “flow domains” 104(1–4):36–60.
doi:10.1016/j.jconhyd.2008.10.003
Köhne JM, Mohanty B, Šimůnek J, Gerke HH (2004) Numerical evaluation of a second-orderwater transfer term for variably saturated dual-permeability models. Water Resour Res 40.
doi:10.1029/2004WR00385
Kosugi K (1996) Lognormal distribution model for unsaturated soil hydraulic properties. Water
Resour Res 32(9):2697–2703
Kung K-JS (1990) Preferential flow in a sandy vadose zone. 2. Mechanism and implications.
Geoderma 46:59–71
Langmuir I (1918) The adsorption of gases on plane surfaces of glass, mica, and platinum. J Am
Chem Soc 40:1361–1403
Leij FJ, Bradford SA (1994) 3DADE: a computer program for evaluating three-dimensional equilibrium solute transport in porous media. Research Report No. 134 US Salinity Laboratory,
USDA, ARS, Riverside, CA
Leij FJ, Toride N (1997) N3DADE: a computer program for evaluating nonequilibrium threedimensional equilibrium solute transport in porous media. Research Report No. 143 US Salinity
Laboratory, USDA, ARS, Riverside, CA
Leij FJ, Priesack E, Schaap MG (2000) Solute transport modeled with Green’s functions withapplication to persistent solute sources. J Contamin Hydrol 41:155–173
Leij FJ, Schaap MG, Arya LM (2002) Indirect methods. In: Dane JH, Topp GC (eds) Methodsof soil analysis, part 1, physical methods, chapter 3.6.3, 3rd edn. SSSA, Madison, WI, pp
1009–1045
Lenhard RJ, Oostrom M, Dane JH (2002) Multi-fluid flow. In: Dane JH, Topp GC (eds) Methodsof soil analysis, part 4, physical methods, SSSA book series 5, Soil Sci Soc Am, Madison, WI,
pp 1537–1619
Lichtner PC (1996) Continuum formulation of multicomponent-multiphase reactive transport. In:
Lichtner PC, Steefel CI, Oelkers EH (eds) Reactive transport in porous media, Reviews inmineralogy. Chapter 1. Mineralogical Society of America, vol 34, pp 1–81
Lichtner PC (2000) FLOWTRAN user’s manual. Los Alamos National Laboratory Document. NM
Lichtner PC, Seth MS (1996) 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
Lindqvist R, Enfield CG (1992) Biosorption of dichlorodiphenyltri-chloroethane and hexachlorobenzene in groundwater and its implications for facilitated transport. Appl Environ
Microbiol 58:1811–1818
Liu CW, Narasimhan TN (1989) Redox-controlled multiple-species reactive chemical transport. 1.
Model development. Water Resour Res 25(5):869–882
Liu C, Szecsody JE, Zachara JM, Ball WP (2000) Use of the generalized integral transform method
for solving equations of solute transport in porous media. Adv Water Resour 23(5):483–492
Mallants D, Jacques D, Vanclooster M, Diels J, Feyen J (1996b) A stochastic approach to simulate
water flow in a macroporous soil. Geoderma 70:299–324
Mallants D, Mohanty BP, Diederik J, Feyen J (1996c) Spatial variability of hydraulic properties ina multi-layered soil profile. Soil Sci 161(3):167–181
Mallants D, Tseng P-H, Toride NA (1997) Timmerman, and Feyen J. Evaluation of multimodal
hydraulic functions in characterizing a heterogeneous field soil. J Hydrol 195:172–199
Mallants D, Vanclooster M, Feyen J (1996a) A transect study on solute transport in a macroporous
soil. Hydrol Process 10:55–70
Mallants D, Volckaert G, Marivoet J (1999) Sensitivity of protective barrier performance tochanges in rainfall rate. Waste Manage 19:467–475
Magee BR, Lion LW, Lemley AT (1991) Transport of dissolved organic macromolecules and theireffect on the transport of phenanthrene in porous media. Environ Sci Technol 25:323–331

846

D. Mallants et al.

Mangold DC, Tsang CF (1991) A summary of subsurface hydrological and hydrochemical models.
Rev Geophys 29(1):51–79
Mansfeldt T, Leyer H, Barmettler K, Kretzschmar R (2004) Cyanide leaching from soil developed
from coking plant purifier waste as influenced by citrate. Vadose Zone J 3:471–479
Massoudieh A, Ginn TR (2009) Modeling colloid-enhanced contaminant transport in stormwater
infiltration basin best management practices. Vadose Zone J 7:1261–1268
Mayer U, Frind EO, Blowes DW (2002) Multicomponent reactive transport modeling in variably
saturated porous media using a generalized formulation for kinetically controlled reactions.
Water Resour Res 38:1174
MDH Engineered Solutions Corp (2003) 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. 07785-2493-0 (printed edition), ISBN No. 0-7785-2494-9 (on-line edition)
Millington RJ, Quirk JM (1961) Permeability of porous solids. Trans Faraday Soc 57:1200–1207
Moldrup O, Poulsen TG, Schjonning P, Olesen T, Yamaguchi T (1998) Gas permeability inundisturbed soils: measurements and predictive models. Soil Sci 163:180–189
Monteith JL (1965) Evaporation and environment. In: Fogg GE (ed) Symposium of the society forexperimental biology. The state and movement of water in living organisms, vol 19. Academic,
NY, pp 205–234
Mualem Y (1976) A new model for predicting the hydraulic conductivity of unsaturated porous
media. Water Resour Res 12(3):513–522
Narasimhan TN (2005) Buckingham, 1907; An appreciation. Vadose Zone J 4:434–441
Nielsen DR, Biggar JW, Erh KT (1973) Spatial variability of field measured soil water properties.
Hilgardia 42:215–259
Nielsen DR, Van Genuchten MTh, Biggar JW (1986) Water flow and solute transport processes in
the unsaturated zone. Water Resour Res 18(9):89S-108S
Noell AL, Thompson JL, Corapcioglu MY, Triay IR (1998) The role of silica colloids on facilitated cesium transport through glass bead columns and modeling. J Contamin Hydrol 31:
23–56
Pachepsky YA, Rawls WJ (eds) (2004) Development of pedotransfer functions in soil hydrology.
Elsevier, Amsterdam
Pachepsky YA, Rawls WJ, Lin HS (2006) Hydropedology and pedotransfer functions. Geoderma
131:308–316
Parker JC, Van Genuchten MTh (1984) Determining transport parameters from laboratory and field
tracer experiments. Bull 84–83, Va Agric Exp St, Blacksburg, VA
Parkhurst DL, Appelo CAJ (1999) User’s guide to PHREEQC (version 2) – a computer program forspeciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.
Water resources investigation, report 99-4259, Denver, CO, USA, 312 pp
Penman HL (1948) Natural evaporation from open water, bare soil, and grass. Proc Roy Soc Lon
Z193:120–156
Pot V, Šimůnek J, Benoit P, Coquet Y, Yra A, Martínez-Cordón M-J (2005) Impact of rainfall
intensity on the transport of two herbicides in undisturbed grassed filter strip soil cores. J.
Contamin Hydrol 8:63–88
Pruess K (1991) Tough2 – a general-purpose numerical simulator for multiphase fluid and heat
flow. Rep LBL-29400, Lawrence Berkeley Laboratory, Berkeley, CA
Pruess K (2004) The TOUGH codes -A family of simulation tools for multiphase flow and transport
processes in permeable media. Vadose Zone J 3:738–746
Pruess K, Wang JSY (1987) Numerical modeling of isothermal and non-isothermal flow in unsaturated fractured rock – a review, In: Evans DD, Nicholson TJ (eds) Flow and transport throughunsaturated fractured rock. Geophysics monograph, vol 42. American Geophysical Union,
Washington, DC, pp 11–18
Rathfelder KM, Lang JR, Abriola LM (2000) A numerical model (MISER) for the simulation of
coupled physical, chemical and biological processes in soil vapor extraction and bioventing
systems. J Contamin Hydrol 43(3–4):239–270
Richards LA (1931) Capillary conduction of fluid through porous mediums. Physics 1:318–333

18

Leaching of Contaminants to Groundwater

847

Ritsema CJ, Dekker LW (2000) Preferential flow in water repellent sandy soils: principles andmodeling implications. J Hydrol 231–232(29):308–319
Russo D, Bresler E (1981) Soil hydraulic properties as stochastic processes: 1. Analysis of field
spatial variability. Soil Sci Soc Am J 45:682–687
Scanlon BR, Christman M, Reedy R, Porro I, Šimůnek J, Flerchinger G (2002) Intercode comparisons for simulating water balance of surficial sediments in semiarid regions. Water Resour Res
38(12):1323. doi:10.1029/2001WR001233, 59.1–59.16
Schaap MG, Leij FJ (2000) Improved prediction of unsaturated hydraulic conductivity with theMualem-Van Genuchten model. Soil Sci Soc Am J 64:843–851
Schaap MG, Leij FJ, Van Genuchten MTh (2001) Rosetta: a computer program for estimating soil
hydraulic parameters with hierarchical pedotransfer functions. J Hydrol 251:163–176
Schaerlaekens J, Mallants D, Šimůnek J, Van Genuchten MTh, Feyen J (1999) Numerical simulation of transport and sequential biodegradation of chlorinated aliphatic hydrocarbons using
CHAIN_2D. J Hydrol Process 13(17):2847–2859
Schijven J, Šimůnek J (2002) Kinetic modeling of virus transport at field scale. J Contamin Hydrol
55(1-2):113-135
Selim HM, Büchter B, Hinz C, Ma L (1992) Modeling the transport and retention of cadmium insoils: multireaction and multicomponent approaches. Soil Sci Soc Am J 56:1004–1015
Selim HM, Davidson JM, Mansell RS (1976) Evaluation of a two-site adsorption-desorption modelfor describing solute transport in soils. Paper presented at proceedings. Summer Computer
Simulation Conference, Nat Sci Found, Washington, DC
Seuntjens P (2000) Reactive solute transport in heterogeneous porous medium: cadmium leachingin acid sandy soils. Ph.D., University of Antwerp, Belgium
Seuntjens P, Tirez K, Šimůnek J, Van Genuchten MTh, Cornelis C, Geuzens P (2001) Aging
effects on cadmium transport in undisturbed contaminated sandy soil columns. J Environ Qual
30:1040–1050Šimůnek J (2005) Models of water flow and solute transport in the unsaturated zone. In: Anderson
MG (ed) The encyclopedia of hydrological sciences. Wiley, Chichester, England, pp 1171–1180
Šimůnek J, Changming H, Pang JL, Bradford SA (2006a) Colloid-facilitated transport in variablysaturated porous media: numerical model and experimental verification. Vadose Zone J 5:1035–
1047
Šimůnek J, Jarvis NJ, Van Genuchten MTh, Gärdenäs A (2003) Nonequilibrium and preferential
flow and transport in the vadose zone: review and case study. J Hydrol 272:14–35
Šimůnek J, Jacques D, Van Genuchten MTh, Mallants D (2006b) Multicomponent geochemicaltransport modeling using the HYDRUS computer software packages. J Am Water Resour Assoc
42(6):1537–1547
Šimůnek J, Suarez DL (1993) Modeling of carbon dioxide transport and production in soil: 1.
Model development. Water Resour Res 29(2):487–497
Šimůnek J, Suarez DL (1994) Two-dimensional transport model for variably saturated porous
media with major ion chemistry. Water Resour Res 30(4):1115–1133Šimůnek J, Suarez DL (1997) Sodic soil reclamation using multicomponent transport modeling.
ASCE J Irrig Drain Eng 123(5):367–376
Šimůnek J, Šejna M, Van Genuchten MTh (1999b) The HYDRUS-2D software package forsimulating two-dimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0, IGWMC – TPS – 53, International Ground Water Modeling Center,
Colorado School of Mines, Golden, CO, 251 pp
Šimůnek J, Šejna M, Van Genuchten MTh (1998a) The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple in variably-saturated media.
Version 2.0, IGWMC – TPS – 70, International Ground Water Modeling Center, Colorado
School of Mines, Golden, CO, 202 pp
Šimůnek J, Suarez DL, Šejna M (1996) The UNSATCHEM software package for simulating
one-dimensional variably saturated water flow, heat transport, carbon dioxide production and848

D. Mallants et al.

transport, and multicomponent solute transport with major ion equilibrium and kinetic chemistry. Version 2.0. Research Report No. 141 US Salinity Laboratory, USDA, ARS, Riverside,CA, 186 ppŠimůnek J, Valocchi AJ (2002) Geochemical transport. In: Dane JH, Topp GC (eds) Methods ofsoil analysis, part 1, physical methods, chapter 6.9, 3rd edn. SSSA, Madison, WI, pp 1511–1536
Šimůnek J, Van Genuchten MTh (2006) Contaminant transport in the unsaturated zone: theory andmodeling, chapter 22. In: Delleur J (ed) The handbook of groundwater engineering, 2nd edn.
CRC, Boca Raton, FL, pp 22.1–22.46
Šimůnek J, Van Genuchten MTh (2008) Modeling nonequilibrium flow and transport with
HYDRUS. Vadose Zone J, Special issue “Vadose zone modeling” 7:782–797
Šimůnek J, Van Genuchten MTh, Šejna M (2005) 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, 270 pp
Šimůnek J, Van Genuchten MTh, Šejna M (2006c) The HYDRUS software package for simulating
two- and three-dimensional movement of water, heat, and multiple solutes in variably-saturated
media, Technical manual, version 1.0, PC Progress, Prague, Czech Republic, 241 pp
Šimůnek J, Van Genuchten MTh, Sejna M (2008) Development and applications of the HYDRUS
and STANMOD software packages and related codes. Vadose Zone J, Special issue “Vadosezone modelling” 7(2):587–600
Šimůnek J, Van Genuchten MTh, Šejna M, Toride N, Leij FJ (1999a) The STANMOD computer software for evaluating solute transport in porous media using analytical solutions ofconvection-dispersion equation. Versions 1.0 and 2.0, IGWMC – TPS – 71, International
Ground Water Modeling Center, Colorado School of Mines, Golden, CO, 32 pp
Šimůnek J, Wendroth O, Van Genuchten MTh (1998b) A parameter estimation analysis of theevaporation method for determining soil hydraulic properties. Soil Sci Soc Am J 62(4):894–905Šimůnek J, Wendroth O, Wypler N, Van Genuchten MTh (2001) Nonequilibrium water flow
characterized from an upward infiltration experiment. Eur J Soil Sci 52(1):13–24
Smolders E, McLaughlin M (1996) Effect of Cl on Cd uptake by Swiss Chard in nutrient solutions.
Plant Soil 179:57–64
Sposito G, Barry DA (1987) On the Dagan model of solute transport I groundwater: foundational
aspects. Water Resour Res 23:1867–1875
Steefel CI, MacQuarrie KTB (1996) Approaches to modeling of reactive transport in porous media.
In: Lichtner PC, Steefel CI, Oelkers EH (eds) Reactive transport in porous media. Reviews inmineralogy. Chapter 2, vol 34, Mineralogical Society of America, pp 83–129
Steefel CI, Yabusaki SB (1996) OS3D/GIMRT, software for multicomponent-multidimensional
reactive transport. User manual and programmer’s guide, PNL-11166, Pacific Northwest
National Laboratory, Richland, WA
Thibault DH, Sheppard MI, Smith PA (1990), A critical compilation and review of default soil
solid/liquid partition coefficients, Kd , for use in environmental assessments. AECL-10125
Thiele-Bruhn S (2003) Pharmaceutical antibiotic compounds in soils – a review. J Plant Nutr Soil
Sci 166:145–167
Tolls J (2001) Sorption of veterinary pharmaceuticals in soils: a review. Environ Sci Technol
35:3397–3406
Toride N, Leij FJ, Van Genuchten MTh (1993) A comprehensive set of analytical solutions fornonequilibrium transport with first-order decay and zero-order production. Water Resour Res
29(7):2167–2182
Toride N, Leij FJ, Van Genuchten MTh (1995) The CXTFIT code for estimating transport parameters from laboratory or field tracer experiments. Version 2.0, Research report no. 137 US
Salinity Laboratory, USDA, ARS, Riverside, CA
Totsche KU, Kögel-Knaber I (2004) Mobile organic sorbent affected contaminant transport in soil:
numerical case studies for enhanced and reduced mobility. Vadose Zone J 3:352–367

18Leaching of Contaminants to Groundwater

849

Ungs M, Wilkinson K, Summers K (1998) Users guide for MYGRT, Version 3.0: software forsimulating migration of organic and inorganic chemicals in grioundwater. EPRI Report TR111748, Electric Power Research Institute, Palo Alto, CA
van Dam JC, Huygen J, Wesseling JG, Feddes RA, Kabat P, van Valsum PEV, Groenendijk P, vanDiepen CA (1997) Theory of SWAP, version 2.0. Simulation of water flow, solute transport andplant growth in the soil- water- atmosphere- plant environment. Department Water Resources,
WAU, Report 71, DLO Winand Staring Centre, Wageningen, Technical Document 45
Vanderborght J, Kasteel R, Herbst M, Javaux M, Thiery D, Vanclooster M, Mouvet C, Vereecken
H (2005) A set of analytical benchmarks to test numerical models of flow and transport in soils.
Vadose Zone J 4:206–181
Van Genuchten MTh (1980a) A closed-form equation for predicting the hydraulic conductivity ofunsaturated soils. Soil Sci Soc Am J 44:892-898
Van Genuchten MTh (1980b) Determining transport parameters from solute displacement experiments. Research report no. 118 US Salinity Laboratory, USDA, ARS, Riverside, CA
Van Genuchten MTh (1981) Non-equilibrium transport parameters from miscible displacement
experiments. Research report no. 119 US Salinity Laboratory, USDA, ARS, Riverside, CA
Van Genuchten MTh (1987) A numerical model for water and solute movement in and below theroot zone. Research Report no. 121 US Salinity Laboratory, USDA, ARS, Riverside, CA
Van Genuchten MTh, Parker JC (1984) Boundary conditions for displacement experiments throughshort laboratory soil columns. Soil Sci Soc Am J 48(4):703-708
Van Genuchten MTh, Šimůnek J (1996) Evaluation of pollutant transport in the unsaturated zone.
In: Rijtema PE, Eliáš V (eds),Regional approaches to water pollution in the environment, NATOASI series: 2. Environment. Kluwer, Dordrecht, The Netherlands, pp 139–172
Van Genuchten MTh, Wagenet RJ (1989) Two-site/two-region models for pesticide transportand degradation: theoretical development and analytical solutions. Soil Sci Soc Am J 53:
1303–1310
Vinten AJA, Yaron B, Nye PH (1983) Vertical transport of pesticides into soil when adsorbed on
suspended particles. J Agric Food Chem 31:662–664
Voegelin A (2001) Competetive sorption and transport of heavy metal cations in soils. Diss ETH
No. 14101, Zürich, Switzerland
Voegelin A, Vulava V, Kretzschmar R (2001) Reaction-based model describing competitive
sorption and transport of Cd, Zn, and Ni in an acidic soil. Environ Sci Technol 35:
1651–1657
Von Gunten HR, Weber VE, Krahenbuhl U (1998) The reactor accident at Chernobyl: a possibility
to test colloid-controlled transport of radionuclides in a shallow aquifer. J Contamin Hydrol
2:237–247
Wagenet RJ, Hutson JL (1987) LEACHM: Leaching Estimation And CHemistry Model, A processbased model of water and movement, transformations, plant uptake and chemical reactions in
the unsaturated zone. Continuum 2. Water Resour Inst, Cornell University, Ithaca, New YorkWang HF, Anderson MP (1982) Introduction to groundwater modeling: finite difference and finite
element methods. Academic, San Diego
White N, Zelazny LW (1986) Charge properties in soil colloids. In: Sparks DL (ed) Soil physicalchemistry. CRC, Boca Raten, FL, pp 39–81
Wheatcraft SW, Cushman JH (1991) Hierarchical approaches to transport in heterogeneous porous
media. Rev Geophys 29(Suppl):261–267
Wild A, Babiker IA (1976) The asymmetric leaching pattern of NO3 and Cl in a loamy sand underfield conditions. J Soil Sci 27:460–466
Willardson LS, Hurst RL (1965) Sample size estimates in permeability studies. J Irrigation Drain
Div, ASCE, 91:1–9
Wilson GV, Albright WH, Gee GW, Fayer MJ, Ogan BD (1999) 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 underContract No. 44-0000-2038, 202 pp

850

D. Mallants et al.

Yeh GT (1981) Analytical transient one-, two-, and three-dimensional simulation of waste transport
in the aquifer system, ORNL-5602, Oak Ridge National Laboratory, Oak Ridge, TN
Yeh GT, Tripathi VS (1990) HYDROGEOCHEM: a coupled model of hydrological and geochemical equilibrium of multicomponent systems. Rep ORNL-6371, Oak Ridge Natl Lab, Oak
Ridge, TN
Yeh GT, Salvage KM, Gwo JP, Zachara JM, Szecsody JE (1998) HYDROBIOGEOCHEM: a coupled model of hydrological transport and mixed biochemical kinetic/equilibrium reactions insaturated-unsaturated media. Rep. ORNL/TM-13668, Oak Ridge Natl Lab, Oak Ridge, TN
Zyvoloski GA, Robinson BA, Dash ZV, Trease LL (1997) Summary of the models and methods
for the FEHM application – a finite element heat- and mass-transfer code. Los Alamos National
Laboratory Report LA-13307-MS, NM