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

# 22

Contaminant Transport in the Unsaturated Zone Theory and Modeling

Ji ̆

rí Šim  ̊ unek University of California Riverside 

### Martinus Th. van Genuchten George E. Brown, Jr, Salinity Laboratory, USDA-ARS 22.1 Introduction

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

.......................... 22 -3 Mass Balance Equation • Uniform Flow • Preferential Flow 22.3 Solute Transport

......................................... 22 -8 Transport Processes • Advection–Dispersion Equations • Nonequilibrium Transport • Stochastic Models • Multicomponent Reactive Solute Transport • Multiphase Flow and Transport • Initial and Boundary Conditions 22.4 Analytical Models

........................................ 22 -22 Analytical Approaches • Existing Models 22.5 Numerical Models

....................................... 22 -27 Numerical Approaches • Existing Models 22.6 Concluding Remarks

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


## 22.1 Introduction
Human society during the past several centuries has created a large number of chemical substances that often find their way into the environment, either intentionally applied during agricultural practices, unintentionally released from leaking industrial and municipal waste disposal sites, or stemming from research or weapons production related activities. As many of these chemicals represent a significant health risk when they enter the food chain, contamination of both surface and subsurface water supplies has become a major issue. Modern agriculture uses an unprecedented number of chemicals, both in plant and animal production. A broad range of fertilizers, pesticides and fumigants are now routinely applied to agricultural lands, thus making agriculture one of the most important sources for non-point source pollution. The same is true for salts and toxic trace elements, which are often an unintended consequence of irrigation in arid and semiarid regions. While many agricultural chemicals are generally beneficial in surface soils, their leaching into the deeper vadose zone and groundwater may pose serious problems. Thus, management processes are being sought to keep fertilizers and pesticides in the root

22 -1JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 2 — #2 22 -2 The Handbook of Groundwater Engineering

zone and prevent their transport into underlying or down-gradient water resources. Agriculture also increasingly uses a variety of pharmaceuticals and hormones in animal production many of which, along with pathogenic microorganisms, are being released to the environment through animal waste. Animal waste and wash water effluent, in turn, is frequently applied to agricultural lands. Potential concerns about the presence of pharmaceuticals and hormones in the environment include: (1) abnormal physiological processes and reproductive impairment; (2) increased incidences of cancer; (3) development of antibiotic resistant bacteria; and (4) increased toxicity of chemical mixtures. While the emphasis above is mostly on non-point source pollution by agricultural chemicals, similar problems arise with point-source pollution from industrial and municipal waste disposal sites, leaking underground storage tanks, chemicals spills, nuclear waste repositories, and mine tailings, among other sources.

Mathematical models should be critical components of any effort to optimally understand and quantify site-specific subsurface water flow and solute transport processes. For example, models can be helpful tools for designing, testing and implementing soil, water, and crop management practices that minimize soil and water pollution. Models are equally needed for designing or remediating industrial waste disposal sites and landfills, or for long-term stewardship of nuclear waste repositories. A large number of specialized numerical models now exist to simulate the different processes at various levels of approximation and for different applications.

Increasing attention is being paid recently to the unsaturated or vadose zone where much of the subsurface contamination originates, passes through, or can be eliminated before it contaminates surface and subsurface water resources. Sources of contamination often can be more easily remediated in the vadose zone, before contaminants reach the underlying groundwater. Other chapters in this Handbook deal with water flow (Chapters xx) and solute transport (Chapters xx) in fully saturated (groundwater) systems and with water flow in the unsaturated zone (Chapter 5). The focus of this chapter thus will be on mathematical descriptions of transport processes in predominantly variably saturated media.

Soils are generally defined as the biologically active layer at the surface of the earth’s crust that is made up of a heterogeneous mixture of solid, liquid, and gaseous material, as well as containing a diverse community of living organisms (Jury and Horton, 2004). The vadose (unsaturated) zone is defined as the layer between the land surface and the permanent (seasonal) groundwater table. While pores between solid grains are fully filled with water in the saturated zone (groundwater), pores in the unsaturated zone are only partially filled with water, with the remaining part of the pore space occupied by the gaseous phase. The vadose zone is usually only partially saturated, although saturated regions may exist, such as when perched water is present above a low-permeable fine-textured (clay) layer or a saturated zone behind the infiltration front during or after a high-intensity rainfall event.

As the transport of contaminants is closely linked with the water flux in soils and rocks making up the vadose zone, any quantitative analysis of contaminant transport must first evaluate water fluxes into and through the vadose zone. Water typically enters the vadose zone in the form of precipitation or irrigation (Figure 22.1), or by means of industrial and municipal spills. Some of the rainfall or irrigation water may be intercepted on the leaves of vegetation. If the rainfall or irrigation intensity is larger than the infiltration capacity of the soil, water will be removed by surface runoff, or will accumulate at the soil surface until it evaporates back to the atmosphere or infiltrates into the soil. Part of this water is returned to the atmosphere by evaporation. Some of the water that infiltrates into the soil profile may be taken up by plant roots and eventually returned to the atmosphere by plant transpiration. The processes of evaporation and transpiration are often combined into the single process of evapotranspiration. Only water that is not returned to the atmosphere by evapotranspiration may percolate to the deeper vadose zone and eventually reach the groundwater table. If the water table is close enough to the soil surface, the process of capillary rise may move water from the groundwater table through the capillary fringe toward the root zone and the soil surface. Because of the close linkage between water flow and solute transport, we will first briefly focus on the physics and mathematical description of water flow in the vadose zone (Section 22.2). An overview is given of the governing equations for water flow in both uniform (Section 22.2.2) and structured (Section 22.2.3) media. This section is followed by a discussion of the governing solute transport equations (Section 22.3),JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 3 — #3 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -3

Transpiration

Rainfall Interception

Evaporation Irrigation

Runoff

Root water

uptake

Deep drainage Capillary rise

Groundwater recharge Water tableUnsaturated zonezoneSaturatedRoot zone Storage

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

again for both uniform (Section 22.3.2) and structured (fractured) (Section 22.3.3) media. We also briefly discuss alternative formulations for colloid (Section 22.3.3.2.2) and colloid-facilitated transport (Section 22.3.3.3), multicomponent geochemical transport (Section 22.3.5), and stochastic approaches for solute transport (Section 22.3.4). This is followed by a discussion of analytical (Section 22.4) and numerical (Section 22.5) approaches for solving the governing flow and transport equations, and an overview of computer models currently available for simulating vadose zone flow and transport processes (Sections 22.4.2 and 22.5.2). 

## 22.2 Variably Saturated Water Flow
In this section, we briefly present the equations governing variably saturated water flow in the subsurface. More details about this topic, including the description of the soil hydraulic properties and their constitutive relationship, are given in Chapter xx. Traditionally, descriptions of variably saturated flow in soils are based on the Richards (1931) equation, which combines the Darcy–Buckingham equation for the fluid flux with a mass balance equation. The Richards equation typically predicts a uniform flow process in the vadose zone, although possibly modified macroscopically by spatially variable soil hydraulic properties (e.g., as dictated by the presence of different soil horizons, but possibly also varying laterally).JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 4 — #4 22 -4 The Handbook of Groundwater Engineering

Unfortunately, the vadose zone can be extremely heterogeneous at a range of scales, from the microscopic (e.g., pore scale) to the macroscopic (e.g., field or larger scale). Some of these heterogeneities can lead to a preferential flow process that macroscopically is very difficult to capture with the standard Richards equation. One obvious example of preferential flow is the rapid movement of water and dissolved solutes through macropores (e.g., between soil aggregates, or created by earthworms or decayed root channels) or rock fractures, with much of the water bypassing (short-circuiting) the soil or rock matrix. However, many other causes of preferential flow exist, such as flow instabilities caused by soil textural changes or water repellency (Hendrickx and Flury, 2001; Šim  ̊ unek et al. 2003; Ritsema and Dekker 2005), and lateral funneling of water due to inclined or other textural boundaries (e.g., Kung 1990). Alternative ways of modeling preferential flow are discussed in a later section. Here we first focus on the traditional approach for uniform flow as described with the Richards equation. 

### 22.2.1 Mass Balance Equation Water flow in variably saturated rigid porous media (soils) is usually formulated in terms of a mass balance equation of the form: ∂θ ∂ t = − ∂ q i ∂ xi − S (22.1)

where θ is the volumetric water content [L 3 L − 3 ], t is time [T], xi is the spatial coordinate [L], q i is the volumetric flux density [LT

− 1 ], and S is a general sink/source term [L 3 L − 3 T
− 1 ], for example, to account for root water uptake (transpiration). Equation 22.1 is often referred to as the mass conservation equation or the continuity equation. The mass balance equation in general states that the change in the water content (storage) in a given volume is due to spatial changes in the water flux (i.e., fluxes in and out of some unit volume of soil ) and possible sinks or sources within that volume. The mass balance equation must be combined with one or several equations describing the volumetric flux density ( q ) to produce the governing equation for variably saturated flow. The formulations of the governing equations for different types of flow (uniform and preferential flow) are all based on this continuity equation.


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

q i = − K ( h )  KA ij ∂ h

∂ x

j + KA iz (22.2)

where K is the unsaturated hydraulic conductivity [LT

− 1 ], and KA ij are components of a dimensionless anisotropy tensor K A (which reduces to the unit matrix when the medium is isotropic). The Darcy– Buckingham equation is formally similar to Darcy’s equation, except that the proportionality constant (i.e., the unsaturated hydraulic conductivity) in the Darcy–Buckingham equation is a nonlinear function of the pressure head (or water content), while K ( h ) in Darcy’s equation is a constant equal to the saturated hydraulic conductivity, K s (e.g., see discussion by Narasimhan [2005]). Combining the mass balance Equation 22.1 with the Darcy–Buckingham Equation22.2 leads to the general Richards equation (Richards, 1931)
∂θ( h )

∂ t = ∂ ∂ xi  K ( h)  KA ij ∂ h

∂ x

j + KA iz 

− S ( h ) (22.3)

This partial differential equation is the equation governing variably saturated flow in the vadose zone. Because of its strongly nonlinear makeup, only a relatively few simplified analytical solutions can be derived. Most practical applications of Equation 22.3 require a numerical solution, which can be obtained using a variety of numerical methods such as finite differences or finite elements (Section 22.5a).JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 5 — #5 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -5

Equation 22.3 is generally referred to as the mixed form of the Richards equation as it contains two dependent variables, that is, the water content and the pressure head. Various other formulations of the Richards equation are possible.


### 22.2.3 Preferential Flow Increasing evidence exists that variably saturated flow in many field soils is not consistent with the uniform flow pattern typically predicted with the Richards equations (Flury et al., 1994; Hendrickx and Flury, 2001). This is due to the presence of macropores, fractures, or other structural voids or biological channels through which water and solutes may move preferentially, while bypassing a large part of the matrix pore-space. Preferential flow and transport processes are probably the most frustrating in terms of hampering accurate predictions of contaminant transport in soils and fractured rocks. Contrary to uniform flow, preferential flow results in irregular wetting of the soil profile as a direct consequence of water moving faster in certain parts of the soil profile than in others. Hendrickx and Flury (2001) defined preferential flow as constituting all phenomena where water and solutes move along certain pathways, while bypassing a fraction of the porous matrix. Water and solutes for these reasons can propagate quickly to far greater depths, and much faster, than would be predicted with the Richards equation describing uniform flow.

The most important causes of preferential flow are the presence of macropores and other structural features, development of flow instabilities (i.e., fingering) caused by profile heterogeneities or water repellency (Hendrickx et al., 1993), and funneling of flow due to the presence of sloping soil layers that redirect downward water flow. While the latter two processes (i.e., flow instability and funneling) are usually caused by textural differences and other factors at scales significantly larger than the pore scale, macropore flow and transport are usually generated at the pore or slightly larger scales, including scales where soil structure first manifests itself (i.e., the pedon scale) (Šim  ̊ unek et al., 2003). Uniform flow in granular soils and preferential flow in structured media (both macroporous soils and fractured rocks) can be described using a variety of single-porosity, dual-porosity, dual-permeability, multi-porosity, and multi-permeability models (Richards, 1931; Pruess and Wang, 1987; Gerke and van Genuchten, 1993a; Gwo et al., 1995; Jarvis, 1998; Šim  ̊ unek et al., 2003, 2005). While single-porositymodels assume that a single pore system exists that is fully accessible to both water and solute, dualporosity and dual-permeability models both assume that the porous medium consists of two interacting pore regions, one associated with the inter-aggregate, macropore, or fracture system, and one comprising the micropores (or intra-aggregate pores) inside soil aggregates or the rock matrix. Whereas dual-porosity models assume that water in the matrix is stagnant, dual-permeability models allow also for water flow within the soil or rock matrix. Figure 22.2 illustrates a hierarchy of conceptual formulations that can be used to model variablysaturated water flow and solute transport in soils. The simplest formulation (Figure 22.2a) is a singleporosity (equivalent porous medium) model applicable to uniform flow in soils. The other models apply in some form or another to preferential flow or transport. Of these, the dual-porosity model of Figure 22.2c assumes the presence of two pore regions, with water in one region being immobile and in the other region mobile. This model allows the exchange of both water and solute between the two regions (Šim  ̊ unek et al., 2003). Conceptually, this formulation views the soil as consisting of a soil matrix containing grains/aggregates with a certain internal microporosity (intra-aggregate porosity) and a macropore or fracture domain containing the larger pores (inter-aggregate porosity). While water and solutes are allowed to move through the larger pores and fractures, they can also flow in and out of aggregates. By comparison, the intra-aggregate pores represent immobile pockets that can exchange, retain, and store water and solutes, but do not contribute to advective (or convective) flow. Models thatassume mobile–immobile flow regions (Figure 22.2b) are conceptually somewhere in between the singleand dual-porosity models. While these models assume that water will move similarly as in the uniform flow models, the liquid phase for purposes of modeling solute transport is divided in terms of mobile and immobile fractions, with solutes allowed to move by advection and dispersion only in the mobile fractionJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 6 — #6 22 -6 The Handbook of Groundwater Engineering

(a) Uniform flow

Water Water Immob

Immob. Immob. Water Water

Solute Solute Solute Solute Mobile

Mobile Mobile Slow

Slow Fast

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

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

FIGURE 22.2 Conceptual models of water flow and solute transport ( θ is the water content, θmo and θ im in (b) and (c) are water contents in the mobile and immobile flow regions, respectively, and θm and θf in (d) are water contents in the matrix and macropore (fracture) regions, respectively).

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

Finally, dual-permeability models (Figure 22.2d) are those in which water can move in both the interand intra-aggregate pore regions (and matrix and fracture domains). These models in various forms are now also becoming increasingly popular (Pruess and Wang, 1987; Gerke and van Genuchten, 1993a; Jarvis, 1994; Pruess, 2004). Available dual-permeability models differ mainly in how they implement water flow in and between the two pore regions (Šim  ̊ unek et al., 2003). Approaches to calculating water flow in the macropores or inter-aggregate pores range from those invoking Poiseuille’s equation (Ahuja and Hebson, 1992), the Green and Ampt or Philip infiltration models (Ahuja and Hebson, 1992), the kinematic wave equation (Germann, 1985; Germann and Beven, 1985; Jarvis, 1994), and the Richards equation (Gerke and van Genuchten, 1993a). Multi-porosity and multi-permeability models (not shown in Figure 22.2) are based on the same concept as dual-porosity and dual-permeability models, but include additional interacting pore regions (e.g., Gwo et al., 1995; Hutson and Wagenet, 1995). These models can be readily simplified to the dual-porosity/permeability approaches. Recent reviews of preferential flow processes and available mathematical models are provided by Hendrickx and Flury (2001) and Šim  ̊ unek et al. (2003), respectively.
22.2.3.1 Dual-Porosity Models Dual-porosity models assume that water flow is restricted to macropores (or inter-aggregate pores and fractures), and that water in the matrix (intra-aggregate pores or the rock matrix) does not move at all. This conceptualization leads to two-region type flow and transport models (van Genuchten and Wierenga, 1976) that partition the liquid phase into mobile (flowing, inter-aggregate), θmo , and immobile (stagnant, intra-aggregate), θim , regions [L3 L − 3 ]: θ = θmo + θim (22.4) The dual-porosity formulation for water flow can be based on a mixed formulation of the Richards Equation 22.3 to describe water flow in the macropores (the preferential flow pathways) and a mass balance equation to describe moisture dynamics in the matrix as follows (Šim  ̊ unek et al., 2003):

∂θmo ( hmo ) ∂ t = ∂ ∂ xi  K ( hmo )  KA ij ∂ hmo ∂ x

j + KA iz 

− S mo ( hmo )

−  w

∂θim ( him ) ∂ t = − Sim ( him ) +  w (22.5)

## JACQ: “4316_c022” — 2006/5/26 — 19:26 — pag e 7 — #7 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -7

where Sim and S mo are sink terms for both regions [T

− 1 ], and  w is the transfer rate for water from the inter- to the intra-aggregate pores [T

− 1 ]. Several of the above dual-porosity features were recently included in the HYDRUS software packages (Šim  ̊ unek et al., 2003, 2005). Examples of their application to a range of laboratory and field data involving transient flow and solute transport are given by Šim  ̊ unek et al. (2001), Zhang et al. (2004), Köhne et al. (2004a, 2005), Kodešová et al. (2005), and Haws et al. (2005). 22.2.3.2 Dual-Permeability Models Different types of dual-permeability approaches may be used to describe flow and transport in structured media. While several models invoke similar governing equations for flow in the fracture and matrix regions, others use different formulations for the two regions. A typical example of the first approach is the work of Gerke and van Genuchten (1993a, 1996) who applied Richards equations to each of two pore regions. The flow equations for the macropore (fracture) (subscript f) and matrix (subscript m) pore systems in their approach are given by:

∂θf ( hf )

∂ t = ∂ ∂ xi  K f ( hf )  KA ij ∂ hf

∂ x

j + KA iz 

− S f ( hf )

−  w w (22.6)

and ∂θm ( hm ) ∂ t = ∂ ∂ xi  K m ( hm )  KA ij ∂ hm

∂ x

j + KA iz 

− S m ( hm ) +  w 1

− w (22.7)
respectively, where w is the ratio of the volumes of the macropore (or fracture or inter-aggregrate) domain and the total soil system [

− ]. This approach is relatively complicated in that the model requires characterization of water retention and hydraulic conductivity functions (potentially of different form) for both pore regions, as well as the hydraulic conductivity function of the fracture–matrix interface. Note that the water contents θf and θm in Equation 22.6 and Equation 22.7 have different meanings than in Equation 22.5 where they represented water contents of the total pore space (i.e., θ = θmo + θim ) , while here they refer to water contents of the two separate (fracture or matrix) pore domains such that θ = wθf + ( 1

− w )θm . 22.2.3.3 Mass Transfer The rate of exchange of water between the macropore and matrix regions,  w , is a critical term in both the dual-porosity model (Equation 22.5) and the dual-permeability approach given by (Equation 22.6) and (Equation 22.7). Gerke and van Genuchten (1993a) assumed that the rate of exchange is proportional to the difference in pressure heads between the two pore regions:

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

in which αw is a first-order mass transfer coefficient [T

− 1 ]. For porous media with well-defined geometries, the first-order mass transfer coefficient, αw , can be defined as follows (Gerke and van Genuchten, 1993b):

αw = β d 2 K a γ

w (22.9)

where d is an effective diffusion path length [L] (i.e., half the aggregate width or half the fracture spacing), β is a shape factor that depends on the geometry [

− ], and γ w ( = 0.4) is a scaling factor [

− ] obtained by matching the results of the first-order approach at the half-time level of the cumulative infiltration curve to the numerical solution of the horizontal infiltration equation (Gerke and van Genuchten, 1993b). Several other approaches based on water content of relative saturation differences have also been used (Šim  ̊ unek et al., 2003).JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 8 — #8 22 -8 The Handbook of Groundwater Engineering 

## 22.3 Solute Transport
Similarly, as shown in Equation 22.1 for water flow, mathematical formulations for solute transport are usually based on a mass balance equation of the form:

∂ C T ∂ t = − ∂ J Ti ∂ xi − φ (22.10)

where C T is the total concentration of chemical in all forms [ML

− 3 ], J Ti is the total chemical 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 or other sources (negative) or sinks (positive) such as plant uptake [ML

− 3 T
− 1 ]. In its most general interpretation, Equation 22.10 allows the chemical to reside in all three phases of the soil (i.e., gaseous, liquid, and solid), permits a broad range of transport mechanisms (including advective transport, diffusion, and hydrodynamic dispersion in both the liquid and gaseous phases), and facilitates any type of chemical reaction that leads to losses or gains in the total concentration.

While the majority of chemicals are present only in the liquid and solid phases, and as such are transported in the vadose zone mostly only in and by water, some chemicals such as many organic contaminants, ammonium, and all fumigants, can have a significant portion of their mass in the gaseous phase and are hence subject to transport in the gaseous phase as well. The total chemical concentration can thus be defined as: C T = ρb s + θ c + ag (22.11)

where ρb is the bulk density [ML

− 3 ], θ is the volumetric water content [L3 L − 3 ], a is the volumetric air content [L3 L − 3 ], and s [MM

− 1 ], c [ML

− 3 ], and g [ML

− 3 ] are concentrations in the solid, liquid, and gaseous phases, respectively. The solid phase concentration represents solutes sorbed onto sorption sites of the solid phase, but can include solutes sorbed onto colloids attached to the solid phase or strained by the porous system, and solutes precipitated onto or into the solid phase. The reaction term φ of Equation 22.10 may represent various chemical or biological reactions that lead to a loss or gain of chemical in the soil system, such as radionuclide decay, biological degradation, anddissolution. In analytical and numerical models these reactions are most commonly expressed using zeroand first-order reaction rates as follows:

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

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

− 1 ], respectively, and γ s [T

− 1 ], γ w [ML

− 3 T
− 1 ], and γ g [ML

− 3 T
− 1 ] are zero-order production constants in the solid, liquid, and gaseous phases, respectively.


### 22.3.1 Transport Processes When a solute is present in both the liquid and gaseous phase, then various transport processes in both of these phases may contribute to the total chemical mass flux:

J T = J l + J g (22.13)

where J l and J g represent solute fluxes in the liquid and gaseous phases [ML

− 2 T
− 1 ], respectively. Note that in Equation 22.13, and further below, we omitted the subscript i accounting for the direction of flow. The three main processes that can be active in both the liquid and gaseous phase are molecular diffusion, hydrodynamic dispersion, and advection (often also called convection). The solute fluxes in the two phasesJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 9 — #9 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -9

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

J l = J lc + J ld + J lh J g = J gc + J gd + J gh (22.14)

where the additional subscripts c, d, and h denote convection (or advection), molecular diffusion, and hydrodynamic dispersion, respectively. 22.3.1.1 Diffusion Diffusion is a result of the random motion of chemical molecules. This process causes a solute to move from a location with a higher concentration to a location with a lower concentration. Diffusive transport can be described using Fick’s law:

J ld = − θξl (θ ) Dw l ∂ c

∂ z = − θ D s l ∂ c

∂ z J gd = − aξg (θ ) Dw g ∂ g

∂ z = − aD s g ∂ g

∂ z (22.15)

where Dw l and Dw g are binary diffusion coefficients of the solute in water and gas [L2 T
− 1 ], respectively; D s l and D s g are the effective diffusion coefficients in soil water and soil gas [L2 T
− 1 ], respectively; and ξl and ξg are tortuosity factors that account for the increased path lengths and decreased cross-sectional areas of the diffusing solute in both phases (Jury and Horton, 2004). As solute diffusion in soil water (air) is severely hampered by both air (water) and solid particles, the tortuosity factor increases strongly with water content (air content). Many empirical models have been suggested in the literature to account for the tortuosity (e.g., Moldrup et al., 1998). Among these, the most widely used model for the tortuosity factor is probably the equation of Millington and Quirk (1961) given by:

ξl (θ ) = θ 7 / 3 θ 2

s (22.16)

where θs is the saturated water content (porosity) [L3 L − 3 ]. A similar equation may be used for the tortuosity factor of the gaseous phase by replacing the water content with the air content. 22.3.1.2 Dispersion Dispersive transport of solutes results from the uneven distribution of water flow velocities within and between different soil pores (Figure 22.3). Dispersion can be derived from Newton’s law of viscosity which states that velocities within a single capillary tube follow a parabolic distribution, with the largest velocity in the middle of the pore and zero velocities at the walls (Figure 22.3a). Solutes in the middle of a pore, for this reason, will travel faster than solutes that are farther from the center. As the distribution of solute ions within a pore depends on their charge, as well as on the charge of pore walls, some solutes may move significantly faster than others. In some situations (i.e., for negatively charged anions in fine-textured soils, leading to anion exclusion), the solute may even travel faster than the average velocity of water (e.g., Nielsen et al., 1986). Using Poiseuille’s law, one can further show that velocities in a capillary tube depend strongly on the radius of the tube, and that the average velocity increases with the radius to the second power. As soils consist of pores of many different radii, solute fluxes in pores of different radii will be significantly different, with some solutes again traveling faster than others (Figure 22.3b).

The above pore-scale dispersion processes lead to an overall (macroscopic) hydrodynamic dispersion process that mathematically can be described using Fick’s law in the same way as molecular diffusion, that is, J lh = − θ D lh ∂ c

∂ z = − θλ v ∂ c ∂ z = − λ q ∂ c ∂ z (22.17)JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 10 — #10 22 -10 The Handbook of Groundwater Engineering

R

0

– R C = C 0 C = 0

Velocity (a) (b)

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

where D lh is the hydrodynamic dispersion coefficient [L2 T
− 1 ], v is the average pore-water velocity [LT

− 1 ], and λ is the dispersivity [L]. The dispersion coefficient in one-dimensional systems has been found to be approximately proportional to the average pore-water velocity, with the proportionality constant generallyreferred as the (longitudinal) dispersivity (Biggar and Nielsen, 1967). The discussion above holds for onedimensional transport; multi-dimensional applications require a more complicated dispersion tensor involving longitudinal and transverse dispersivities (e.g., Bear, 1972).

Dispersivity is a transport property that is relatively difficult to measure experimentally. Estimates are usually obtained by fitting measured breakthrough curves with analytical solutions of the advection– dispersion equation (discussed further below). The dispersivity often changes with the distance over which solute travels. Values of the longitudinal dispersivity typically range from about 1 cm for packed laboratory columns, to about 5 or 10 cm for field soils. Longitudinal dispersivities can be significantly larger (on the order of hundreds of meters) for regional groundwater transport problems (Gelhar et al., 1985). If no other information is available, a good first approximation is to use a value of one-tenth of the transport distance for the longitudinal dispersivity (e.g., Anderson, 1984), and a value of one-hundreds of the transport distance for the transverse dispersivity. 22.3.1.3 Advection Advective transport refers to solute being transported with the moving fluid, either in the liquid phase ( J

lc ) or the gas phase ( J gc ) , that is,

J lc = qc J gc = J g g (22.18)

where J g is the gaseous flux density [LT

− 1 ]. Advective transport in the gaseous phase is often neglected as its contribution in many applications is negligible compared to gaseous diffusion.

The total solute flux density in both the liquid and gaseous phases is obtained by incorporating contributions from the various transport processes into Equation 22.14 to obtain

J l = qc

− θ D s l ∂ c

∂ z − θ D lh ∂ c

∂ z = qc

− θ De ∂ c ∂ z J g = − aD s g ∂ g

∂ z (22.19)

where De is the effective dispersion coefficient [L2 T
− 1 ] that accounts for both diffusion and hydrodynamic dispersion. Dispersion in most subsurface transport problems dominates molecular diffusion in the liquid phase, except when the fluid velocity becomes relatively small or negligible. Notice that Equation 22.19 neglects advective and dispersive transport in the gaseous phase.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 11 — #11 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -11 

### 22.3.2 Advection–Dispersion Equations 22.3.2.1 Transport Equations The equation governing transport of dissolved solutes in the vadose zone is obtained by combining the solute mass balance (Equation 22.10) with equations defining the total concentration of the chemical (Equation 22.11) and the solute flux density (Equation 22.19) to give
∂(ρb s + θ c + ag) ∂ t = ∂ ∂ xi θ De ij ∂ c ∂ x

j + ∂ ∂ xi  aD s g ij ∂ g

∂ x

j − ∂  q i c ∂ xi − φ (22.20)

Notice that this equation is again written for multidimensional transport, and that De ij and D s ij are thus components of the effective dispersion tensor in the liquid phase and a diffusion tensor in the gaseous phase [L2 T
− 1 ], respectively.

Many different variants of Equation 22.20 can be found in the literature. For example, for one-dimensional transport of nonvolatile solutes, the equation simplifies to
∂(ρb s + θ c ) ∂ t = ∂(θ Rc ) ∂ t = ∂ ∂ z θ De ∂ c ∂ z − ∂( qc ) ∂ z − φ (22.21)

where q is the vertical water flux density [LT

− 1 ] and R is the retardation factor [

− ]

R = 1 + ρb θ d s ( c ) d c (22.22)
For transport of inert, nonadsorbing solutes during steady-state water flow we obtain

∂ c ∂ t = De ∂ 2 c ∂ z 2 − v ∂ c ∂ z (22.23)
The above equations are usually referred to as advection–dispersion equations (ADEs). 22.3.2.2 Linear and Nonlinear Sorption The ADE given by Equation 22.20 contains three unknown concentrations (those for the liquid, solid, and gaseous phases), while Equation 22.21 contains two unknowns. To be able to solve these equations, additional information is needed that somehow relates these concentrations to each other. The most common way is to assume instantaneous sorption and to use adsorption isotherms to relate the liquid and adsorbed concentrations. The simplest form of the adsorption isotherm is the linear isotherm given by

s = K d c (22.24)

where K d is the distribution coefficient [L3 M − 1 ]. One may verify that substitution of this equation into Equation 22.22 leads to a constant value for the retardation factor (i.e., R = 1 + ρb K d /θ ). While the use of a linear isotherm greatly simplifies the mathematical description of solute transport, sorption and exchange are generally nonlinear and most often depend also on the presence of competing species in the soil solution. The solute retardation factor for nonlinear adsorption is not constant, as is the case for linear adsorption, but changes as a function of concentration. Many models have been used in the past to describe nonlinear sorption. The most commonly used nonlinear sorption models are those by Freundlich (1909) and Langmuir (1918) given by:

s = K f c β (22.25)

s = K d c 1 + η c (22.26)JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 12 — #12 22 -12 The Handbook of Groundwater Engineering

0 0.2 0.4 0.6 0.8

1 1.2 1.4

Dissolved concentration [M/L 3] 0 0.5 1 1.5 0 0.5 1 1.5Sorbed concentration [M/M]

0.5 0.75 1 1.5 2 2.5 3 5 Increasing

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

Dissolved concentration [M/L 3]Sorbed concentration [M/M] 0 0.5 0.75 1 1.5 2 2.5 3 5 Increasing (a) (b)

FIGURE 22.4 Plots of the Freundlich adsorption isotherm given by (22.25), with K d = 1 and β given in the caption (a), and the Langmuir adsorption isotherm given by (22.26), with K d = 1 and η given in the caption (b).

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

s = k 1 c + k 2 Linear Lapidus and Amundson (1952) Lindstrom et al. (1967)

s = k 1 c k 3 1 + k 2 c k 3 Freundlich–Langmuir Sips (1950), Šim  ̊ unek et al. (1994, 2005)

s = k 1 c 1 + k 2 c + k 3 c 1 + k 4 c Double Langmuir Shapiro and Fried (1959)

s = k 1 c c k 2 / k3 Extended Freundlich Sibbesen (1981)

s = k 1 c 1 + k 2 c + k 3 √c Gunary Gunary (1970)

s = k 1 c k 2

− k 3 Fitter–Sutton Fitter and Sutton (1975)

s = k 1 { 1 − [ 1 + k 2 c k 3 ] k

4 } Barry Barry (1992)

s = RT k 1 ln ( k 2 c ) Temkin Bache and Williams (1971)

s = k 1 c exp (

− 2 k 2 s) Lindstrom et al. (1971) van Genuchten et al. (1974) s
s T = c [ c + k 1 ( cT − c ) exp { k 2 ( cT − 2 c ) }] − 1 Modified Kielland Lai and Jurinak (1971)

k 1 , k 2 , k 3 , k 4 : empirical constants; R : universal gas constant; T : absolute temperature; c T : maximum solute concentration; s T : maximum adsorbed concentration.

Source: Adapted from van Genuchten, M.Th. and Šimunek, J. In P.E. Rijtema and V. Eliáš (Eds.), Regional Approaches to Water Pollution in the Environment , NATO ASI Series: 2. Environment. Kluwer, Dordrecht, The Netherlands, pp. 139–172, 1996.

respectively, where K f [M

− β L − 3 β ] and β [

− ] are coefficients in the Freundlich isotherm, and η [L3 M − 1 ] is a coefficient in the Langmuir isotherm. Examples of linear, Freundlich and Langmuir adsorption isotherms are given in Figure 22.4. Table 22.1 lists a range of linear and other sorption models frequently used in solute transport studies.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 13 — #13 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -13

22.3.2.3 VolatilizationVolatilization is increasingly recognized as an important process affecting the fate of many organic chemicals, including pesticides, fumigants, and explosives in field soils (Jury et al., 1983, 1984; Glotfelty and Schomburg, 1989). While many organic pollutants dissipate by means of chemical and microbiological degradation, volatilization may be equally important for volatile substances, such as certain pesticides. The volatility of pesticides is influenced by many factors, including the physicochemical properties of the chemical itself as well-such environmental variables as temperature and solar energy. Even though only asmall fraction of a pesticide may exist in the gas phase, air-phase diffusion rates can sometimes be comparable to liquid-phase diffusion as gas-phase diffusion is about four orders of magnitude greater than liquid phase diffusion. The importance of gaseous diffusion relative to other transport processes depends also on the climate. For example, while transport of MTBE (gasoline oxygenate) is generally dominated by liquid advection in humid areas, gaseous diffusion may be equally or more important in arid climates; this even though only about 2% of MTBE may be in the gas phase. The general transport equation given by Equation 22.20 can be simplified considerably when assuming linear equilibrium sorption and volatilization such that the adsorbed ( s) and gaseous ( g ) concentrations are linearly related to the solution concentration ( c ) through the distribution coefficients, K d in (22.24) and K H , that is, g = K H c (22.27) respectively, where K H is the dimensionless Henry’s constant [

− ]. Equation 22.20 for one-dimensional transport then has the form:
∂(ρb K d + θ + aK H ) c ∂ t = ∂ ∂ z θ De ∂ c ∂ z + ∂ ∂ z  aD s g K H ∂ c

∂ z − ∂  qc ∂ x − φ (22.28)

or ∂θ Rc ∂ t = ∂ ∂ z θ D E ∂ c

∂ z − ∂  qc ∂ x − φ (22.29)

where the retardation factor R [

− ] and the effective dispersion coefficient D E [L2 T
− 1 ] are defined as follows:

R = 1 + ρb K d + aK H θ

D

E = De + aD s g K H θ (22.30)

Jury et al. (1983, 1984) provided for many organic chemicals their distribution coefficients K d , Henry’s constants K H , and calculated percent mass present in each phase. 

### 22.3.3 Nonequilibrium Transport As equilibrium solute transport models often fail to describe experimental data, a large number of diffusion-controlled physical nonequilibrium and chemical-kinetic models have been proposed and usedto describe the transport of both non-adsorbing and adsorbing chemicals. Attempts to model nonequilibrium transport usually involve relatively simple first-order rate equations. Nonequilibrium models have used the assumptions of two-region (dual-porosity) type transport involving solute exchange between mobile and immobile liquid transport regions, and one-, two- or multi-site sorption formulations (e.g., Nielsen et al., 1986; Brusseau, 1999). Models simulating the transport of particle-type pollutants, such as colloids, viruses, and bacteria, often also use first-order rate equations to describe such processes asattachment, detachment, and straining. Nonequilibrium models generally have resulted in better descriptions of observed laboratory and field transport data, in part by providing additional degrees of freedom for fitting observed concentration distributions.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 14 — #14 22 -14 The Handbook of Groundwater Engineering
22.3.3.1 Physical Nonequilibrium 22.3.3.1.1 Dual-Porosity and Mobile–Immobile Water Models The two-region transport model (Figures 22.2b and Figure 22.2c) assumes that the liquid phase can be partitioned into distinct mobile (flowing) and immobile (stagnant) liquid pore regions, and that solute exchange between the two liquid regions can be modeled as a first-order exchange process. Using the same notation as before, the two-region solute transport model is given by (van Genuchten and Wagenet, 1989; Toride et al., 1993):

∂θmo c mo ∂ t + ∂f ρ s mo ∂ t = ∂ ∂ z θmo D mo ∂ c mo ∂ z − ∂ qc mo ∂ z − φmo −  s
∂θim c im ∂ t + ∂( 1

− f )ρ s im ∂ t = − φim +  s (22.31)

for the mobile (macropores, subscript mo) and immobile (matrix, subscript im) domains, respectively, where f is the dimensionless fraction of sorption sites in contact with the mobile water [

− ], φmo and φim are reactions in the mobile and immobile domains [ML3 T
− 1 ], respectively, and  s is the solute transfer rate between the two regions [ML 3 T
− 1 ]. Notice that the same equations Equation 22.31 can be used to describe solute transport using both the mobile–immobile and dual-porosity models shown in Figures 22.2b and Figure 22.2c, respectively. 22.3.3.1.2 Dual-Permeability Model Analogous to Equations 22.6 and Equation 22.7 for water flow, the dual-permeability formulation for solute transport can be based on advection–dispersion type equations for transport in both the fracture and matrix regions as follows (Gerke and van Genuchten, 1993a):

∂θf c f ∂ t + ∂ρ s f ∂ t = ∂ ∂ z θf D f ∂ c f ∂ z − ∂ qf c f ∂ z − φf −  s w (22.32)

∂θm c m ∂ t + ∂ρ s m ∂ t = ∂ ∂ z θm D m ∂ c m ∂ z − ∂ qm c m ∂ z − φm −  s 1

− w (22.33) where the subscript f and m refer to the macroporous (fracture) and matrix pore systems, respectively; φf and φm represent sources or sinks in the macroporous and matrix domains [ML 3 T
− 1 ], respectively; and

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

− ]. Equation 22.32 and Equation 22.33 assume complete advective–dispersive type transport descriptions for both the fractures and the matrix. Several authors simplified transport in the macropore domain, for example, by considering only piston displacement of solutes (Ahuja and Hebson, 1992; Jarvis, 1994). 22.3.3.1.3 Mass Transfer The transfer rate,  s , in Equation 22.31 for solutes between the mobile and immobile domains in the dual-porosity models can be given as the sum of diffusive and advective fluxes, and can be written as

 s = α s ( c mo − c im ) +  w c ∗ (22.34)

where c ∗ is equal to c mo for  w > 0 and c im for  w < 0, and α s is the first-order solute mass transfer coefficient [T

− 1 ]. Notice that the advection term of Equation 22.34 is equal to zero for the mobile– immobile model (Figure 22.2b) as the immobile water content in this model is assumed to be constant. However,  w may have a nonzero value in the dual-porosity model depicted in Figure 22.2c.

The transfer rate,  s , in Equation 22.32 and Equation 22.33 for solutes between the fracture and matrix regions is also usually given as the sum of diffusive and advective fluxes as follows (e.g., Gerke and van Genuchten, 1996):  s = α s ( 1

− w m )( c f − c m ) +  w c ∗ (22.35)JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 15 — #15 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -15

in which the mass transfer coefficient, α s [T

− 1 ], is of the form:

α s = β d 2 Da (22.36)

where Da is an effective diffusion coefficient [L2 T
− 1 ] representing the diffusion properties of the fracture– matrix interface. Still more sophisticated models for physical nonequilibrium transport may be formulated. For example, Pot et al. (2005) and Köhne et al. (2006) considered a dual-permeability model that divides the matrix domain further into mobile and immobile subregions and used this model successfully to simulate bromide transport in laboratory soil columns at different flow rates or for transient flow conditions, respectively. 22.3.3.2 Chemical Nonequilibrium 22.3.3.2.1 Kinetic Sorption ModelsAn alternative to expressing sorption as an instantaneous process using algebraic equations (e.g., Equation 22.24, Equation 22.25 or Equation 22.26) is to describe the kinetics of the reaction using ordinary differential equations. The most popular and simplest formulation of a chemically controlled kinetic reaction arises when first-order linear kinetics is assumed:

∂ s ∂ t = αk ( K d c

− s) (22.37)

where αk is a first-order kinetic rate coefficient [T

− 1 ]. Several other nonequilibrium adsorption expressions were also used in the past (see Table 2 in van Genuchten and Šim  ̊ unek, 1996). Models based on this and other kinetic expressions are often referred to as one-site sorption models.

As transport models assuming chemically controlled nonequilibrium (one-site sorption) generally did not result in significant improvements in their predictive capabilities when used to analyze laboratory column experiments, the one-site first-order kinetic model was further expanded into a two-site sorption concept that divides the available sorption sites into two fractions (Selim et al., 1976; van Genuchten and Wagenet, 1989). In this approach, sorption on one fraction (type-1 sites) is assumed to be instantaneous while sorption on the remaining (type-2) sites is considered to be time-dependent. Assuming a linear sorption process, the two-site transport model is given by (van Genuchten and Wagenet, 1989)

∂(

f ρb K d + θ ) c ∂ t = ∂ ∂ z θ De ∂ c ∂ z − ∂  qc ∂ z − φe ∂ s

k ∂ t = αk [ ( 1

− f) K d c

− s k ] − φk (22.38)

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 nonequilibrium phases, respectively, and the subscript k refers to kinetic (type-2) sorption sites. Note that if f = 0, the two-site sorption model reduces to the one-site fully kinetic sorption model (i.e., when only type-2 kinetic sites are present). On the other hand, if f = 1, the two-site sorption model reduces to the equilibrium sorption model for which only type-1 equilibrium sites are present.
22.3.3.2.2 Attachment/Detachment Models Additionally, transport equations may include provisions for kinetic attachment/detachment of solutes to the solid phase, thus permitting simulations of the transport of colloids, viruses, and bacteria. The transport of these constituents is generally more complex than that of other solutes in that they are affected by such additional processes as filtration, straining, sedimentation, adsorption and desorption, growth, and inactivation. Virus, colloid, and bacteria transport and fate models commonly employ aJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 16 — #16 22 -16 The Handbook of Groundwater Engineering

modified form of the ADE, in which the kinetic sorption equations are replaced with equations describing kinetics of colloid attachment and detachment as follows:

ρ ∂ s

∂ t = θ k a ψ c

− k d ρ s (22.39)

where c is the (colloid, virus, bacteria) concentration in the aqueous phase [N c L − 3 ], s is the solid phase (colloid, virus, bacteria) concentration [N c M − 1 ], in which N c is a number of (colloid) particles, k a is the first-order deposition (attachment) coefficient [T

− 1 ], k d is the first-order entrainment (detachment) coefficient [T

− 1 ], and ψ is a dimensionless colloid retention function [

− ]. The attachment and detachment coefficients in Equation 22.39 have been found to strongly depend upon water content, with attachment significantly increasing as the water content decreases.

To simulate reductions in the attachment coefficient due to filling of favorable sorption sites, ψ is sometimes assumed to decrease with increasing colloid mass retention. A Langmuirian dynamics (Adamczyk et al., 1994) equation has been proposed for ψ to describe this blocking phenomenon:

ψ = s max − s

s max = 1

− s
s max (22.40)

in which s max is the maximum solid phase concentration [N c M − 1 ]. A similar equation as Equation 22.39 was used by Bradford et al. (2003, 2004) to simulate the process of pore straining. Bradford et al. (2003, 2004) hypothesized that the influence of straining and attachment processes on colloid retention should be separated into two distinct components. They suggested the following depth-dependent blocking coefficient for the straining process:

ψ =  d c + z

− z 0 d

c − β (22.41)

where d c is the diameter of the sand grains [L], z 0 is the coordinate of the location where the straining process starts [L] (the surface of the soil profile, or interface between soil layers), and β is an empirical factor (Bradford et al., 2003) [

− ].

The attachment coefficient is often calculated using filtration theory (Logan et al., 1995), a quasiempirical formulation in terms of the median grain diameter of the porous medium (often termed the collector), the pore-water velocity, and collector and collision (or sticking) efficiencies accounting for colloid removal due to diffusion, interception, and gravitational sedimentation (Rajagopalan and Tien, 1976; Logan et al., 1995): k a = 3 ( 1

− θ ) 2 d

c ηα v (22.42)
where d c is the diameter of the sand grains [L], α is the sticking efficiency (ratio of the rate of particles that stick to a collector to the rate they strike the collector) [

− ], v is the pore-water velocity [LT

− 1 ], and η is the single-collector efficiency [

− ]. In related studies, Schijven and Hassanizadeh (2000) and Schijven and Šim  ̊unek (2002) used a twosite sorption model based on two equations (22.39) to successfully describe virus transport at both the laboratory and field scale. Their model assumed that the sorption sites on the solid phase can be divided into two fractions with different properties and various attachment and detachment rate coefficients. 22.3.3.3 Colloid-Facilitated Solute Transport There is an increasing evidence that many contaminants, including radionuclides (Von Gunten et al., 1988; Noell et al., 1998), pesticides (Vinten et al., 1983; Kan and Tomson, 1990, Lindqvist and Enfield, 1992), heavy metals (Grolimund et al., 1996), viruses, pharmaceuticals (Tolls, 2001; Thiele-Bruhn, 2003), hormones (Hanselman et al., 2003), and other contaminants (Magee et al., 1991; Mansfeldt et al., 2004) are transported in the subsurface not only with moving water, but also sorbed to mobile colloids. As manyJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 17 — #17 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -17

colloids and microbes are negatively charged and thus electrostatically repelled by negatively charged solid surfaces, which may lead to an anion exclusion process, their transport may be slightly enhanced relative to fluid flow. Size exclusion may similarly enhance the advective transport of colloids by limiting their presence and mobility to the larger pores (e.g., Bradford et al., 2003). The transport of contaminants sorbed to mobile colloids can thus significantly accelerate their transport relative to more standard advection– transport descriptions. Colloid-facilitated transport is a relatively complicated process that requires knowledge of water flow, colloid transport, dissolved contaminant transport, and colloid-contaminant interaction. Transport and mass balance equations, hence, must be formulated not only for water flow and colloid transport, but also for the total contaminant, for contaminant sorbed kinetically or instantaneously to the solid phase, and for contaminant sorbed to mobile colloids, to colloids attached to the soil solid phase, and to colloids accumulating at the air–water interface. Development of such a model is beyond the scope of this chapter. We refer interested readers to several manuscripts dealing with this topic: Mills et al. (1991), Corapcioglu and Jiang (1993), Corapcioglu and Kim (1995), Jiang and Corapcioglu (1993), Noell et al. (1998), Saiers et al. (1996), Saiers and Hornberger (1996), van de Weerd et al. (1998), and van Genuchten and Šim  ̊ unek (2004).


### 22.3.4 Stochastic Models Much evidence suggests that solutions of classical solute transport models, no matter how refined to include the most relevant chemical and microbiological processes and soil properties, often still fail to accurately describe transport processes in most natural field soils. A major reason for this failure is the fact that the subsurface environment is overwhelmingly heterogeneous. Heterogeneity occurs at a hierarchy ofspatial and time scales (Wheatcraft and Cushman, 1991), ranging from microscopic scales involving timedependent chemical sorption and precipitation/dissolution reactions, to intermediate scales involving the preferential movement of water and chemicals through macropores or fractures, and to much larger scales involving the spatial variability of soils across the landscape. Subsurface heterogeneity can be addressed in terms of process-based descriptions which attempt to consider the effects of heterogeneity at one or several scales. It can also be addressed using stochastic approaches that incorporate certain assumptions about the transport process in the heterogeneous system (e.g., Sposito and Barry, 1987; Dagan, 1989). In this section we briefly review several stochastic transport approaches, notably those using stream tube models and the transfer function approach. 22.3.4.1 Stream Tube Models The downward movement of chemicals from the soil surface to an underlying aquifer may be described stochastically by viewing the field as a series of independent vertical columns, often referred to as “stream tubes” (Figure 22.5), while solute mixing between the stream tubes is assumed to be negligible. Transport

V

1 V 2 V 3 V n
– 1 V n

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

in each tube may be described deterministically with the standard ADE, or modifications thereof to include additional geochemical and microbiological processes. Transport at the field scale is then implemented by considering the column parameters as realizations of a stochastic process, having a random distribution (Toride et al., 1995). Early examples are by Dagan and Bresler (1979) and Bresler and Dagan (1979) who assumed that the saturated hydraulic conductivity had a lognormal distribution.

The stream tube model was implemented into the CXTFIT 2.0 program (Toride et al., 1995) for a variety of transport scenarios in which the pore-water velocity in combination with either the dispersion coefficient, De , the distribution coefficient for linear adsorption, K d , or the first-order rate coefficient for nonequilibrium adsorption, αk , are stochastic variables (Toride et al., 1995). 22.3.4.2 Transfer Function ModelsJury (1982) developed an alternative formulation for solute transport at the field scale, called the transfer function model. This model was developed based on two main assumptions about the soil system (a) the solute transport is a linear process, and (b) the solute travel time probabilities do not change over time. These two assumptions lead to the following transfer function equation that relates the solute concentration at the outflow end of the system with the time-dependent solute input into the system:

c out ( t ) = t 

0 c in ( t

− t )f ( t ) d t (22.43)

The outflow at time t , c out ( t ) [ML

− 3 ], consists of the superposition of solute added at all times less than t , c in ( t

− t ) [ML

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

− 1 ] (Jury and Horton, 2004). One important advantage of the transfer function approach is that it does not require knowledge of the various transport processes within the flow domain. Different model distribution functions can be used for the travel-time probability density function f ( t ) in Equation 22.43. Most commonly used (Jury and Sposito, 1985) are the Fickian probability density function

f ( t ) = L 2 √ π Dt 3 exp 

− ( L

− vt) 2 4 Dt  (22.44)

and the lognormal distribution

f ( t ) = 1 √2 πσ t exp  (ln t

− μ) 2 2 σ 2  (22.45)

where D [L2 T
− 1 ], v [LT

− 1 ], μ , and σ are model parameters, and L is the distance from the inflow boundary to the outflow boundary [L]. To accommodate conditions when the water flux through the soil system is not constant, Jury (1982) expressed the travel-time pdf as a function of the cumulative net applied water I :

I = t 

0 q ( t ) d t (22.46)

leading to the following transfer function equation:

cout ( I ) = I ( t ) 

0 c in ( I

− I )f ( I ) d I (22.47)JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 19 — #19 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -19 

### 22.3.5 Multicomponent Reactive Solute Transport The various mathematical descriptions of solute transport presented thus far all considered solutes thatwould move independently of other solutes in the subsurface. In reality, the transport of reactive contaminants is more often than not affected by many often interactive physico-chemical and biological processes. Simulating these processes requires a more comprehensive approach that couples the physical processes of water flow and advective–dispersive transport with a range of biogeochemical processes. The soil solution is always a mixture of many ions that may be involved in mutually dependent chemical processes, such as complexation reactions, cation exchange, precipitation–dissolution, sorption–desorption, volatilization, redox reactions, and degradation, among other reactions (Šim  ̊ unek and Valocchi, 2002). The transportand transformation of many chemical contaminants is further mediated by subsurface aerobic or anaerobic bacteria. Bacteria catalyze redox reactions in which organic compounds (e.g., hydrocarbons) act as the electron donor and inorganic substances (oxygen, nitrate, sulfate, or metal oxides) as the electron acceptor. By catalyzing such reactions, bacteria gain energy and organic carbon to produce new biomass. These and related processes can be simulated using integrated reactive transport codes that couple the physical processes of water flow and advective–dispersive solute transport with a range of biogeochemical processes. This section reviews various modeling approaches for such multicomponent transport systems. 22.3.5.1 Components and Reversible Chemical Reactions Multi-species chemical equilibrium systems are generally defined in terms of components. Components may be defined as a set of linearly independent chemical entities such that every species in the system can be uniquely represented as the product of a reaction involving only these components (Westall et al., 1976). As a typical example, the chemical species CaCO 0 3

Ca 2 + + CO 2

− 3 ⇔ CaCO 0 3 (22.48)

consists of the two components Ca2 + and CO 2

− 3 . Reversible chemical reaction processes in equilibrium systems are most often represented using mass action laws that relate thermodynamic equilibrium constants to activities (the thermodynamic effective concentration) of the reactants and products (Mangold and Tsang, 1991; Appelo and Postma, 1993; Bethke, 1996). For example, the reaction

bB + cC ⇔ dD + eE (22.49)

where b and c are the number of moles of substances B and C that react to yield d and e moles of products D and E , is represented at equilibrium by the law of mass action

K = a d D a e E
a b

B a c C (22.50)

where K is a temperature-dependent thermodynamic equilibrium constant, and ai is the ion activity, being defined as the product of the activity coefficient ( γ i ) and the ion molality ( mi ) , that is, ai = γ i mi . Single-ion activity coefficients may be calculated using either the Davies equation, an extended version of the Debye–Hückel equation (Truesdell and Jones, 1974), or by means of Pitzer expressions (Pitzer, 1979). Equation 22.50 can be used to describe all of the major chemical processes, such as aqueous complexation, sorption, precipitation–dissolution, and acid–base and redox reactions, provided that the local chemical equilibrium assumption is valid (Šim  ̊ unek and Valocchi, 2002).JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 20 — #20 22 -20 The Handbook of Groundwater Engineering
22.3.5.2 Complexation Equations for aqueous complexation reactions can be obtained using the law of mass action as follows (e.g., Yeh and Tripathi, 1990; Lichtner, 1996):

xi = K x i γ x

i N a 
k = 1 (γ a k c

k ) a x ik i = 1, 2,

... , M x (22.51)

where xi is the concentration of the ith complexed species, K x i is the thermodynamic equilibrium constant of the ith complexed species, γ x i is the activity coefficient of the ith complexed species, N a is the number of aqueous components, c k is the concentration of the kth aqueous component, γ a k is the activity coefficient of the kth aqueous component species, a x ik is the stoichiometric coefficient of the kth aqueous component in the ith complexed species, M x is the number of complexed species, and subscripts and superscripts x and a refer to complexed species and aqueous components, respectively. 22.3.5.3 Precipitation and Dissolution Equations describing precipitation–dissolution reactions are also obtained using the law of mass action, but contrary to the other processes are represented by inequalities rather than equalities, as follows (Šim  ̊ unek and Valocchi, 2002):

K p

i ≥ Q p i = N a 
k = 1 (γ a k c

k ) ap ik i = 1, 2,

... , M p (22.52)

where M p is the number of precipitated species, K p i is the thermodynamic equilibrium constant of the ith precipitated species, that is, the solubility product equilibrium constant, Q p i is the ion activity product of the ith precipitated species, and a p ik is the stoichiometric coefficient of the kth aqueous component in the ith precipitated species. The inequality in Equation 22.52 indicates that a particular precipate is formed only when the solution is supersaturated with respect to its aqueous components. If the solution is undersaturated then the precipitated species (if it exists) will dissolve to reach equilibrium conditions. 22.3.5.4 Cation Exchange Partitioning between the solid exchange phase and the solution phase can be described with the general exchange equation (White and Zelazny, 1986):

z j · A zi + X zi + zi · B z j + ⇔ zi · B z j + X z j + z j · A zi + (22.53)

where A and B are chemical formulas for particular cation (e.g., Ca2 + or Na + ) , X refers to an “exchanger” site on the soil, and zi is the valence of species. The mass action equation resulting from this exchange reaction is

K ij =   c z j + j γ a j c z j + j   zi   γ a i c z + i i

c z +

i i   z

j (22.54)

where K ij is the selectivity coefficient, and c i is the exchanger-phase concentration of the ith component (expressed in moles per mass of solid). 22.3.5.5 Coupled System of Equations Once the various chemical reactions are defined, the final system of governing equations usually consists of several partial differential equations for solute transport (i.e., ADEs for each component) plus a set of nonlinear algebraic and ordinary differential equations describing the equilibrium and kinetic reactions, respectively. Each chemical and biological reaction must be represented by the corresponding algebraic or ordinary differential equations depending upon the rate of the reaction. Since the reaction of oneJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 21 — #21 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -21

species depends upon the concentration of many other species, the final sets of equations typically are tightly coupled. For complex geochemical systems, consisting of many components and multidimensional transport, numerical solution of these coupled equations is challenging (Šim  ̊ unek and Valocchi, 2002). As an alternative, more general models have recently been developed that also more loosely couple transport and chemistry using a variety of sequential iterative or non-iterative operator-splitting approaches (e.g., Bell and Binning, 2004; Jacques and Simunek, 2005). Models based on these various approaches are further discussed in Section 22.5b. 

### 22.3.6 Multiphase Flow and Transport While the transport of solutes in variably saturated media generally involves two phases (i.e., the liquid and gaseous phases, with advection in the gaseous phase often being neglected), many contamination problems also increasingly involve nonaqueous phase liquids (NAPLs) that are often only slightly miscible with water. Nonaqueous phase liquids may consist of single organic compounds such as many industrial solvents, or of a mixture of organic compounds such as gasoline and diesel fuel. Some of these compounds can be denser than water (commonly referred to as DNAPLs) or lighter than water (LNAPLs). Their fate and dynamics in the subsurface is affected by a multitude of compound-specific flow and multicomponent transport processes, including interphase mass transfer and exchange (also with the solid phase). Multiphase fluid flow models generally require flow equations for each fluid phase (water, air, NAPL). Two-phase air–water systems hence could be modeled also using separate equations for air and water. This shows that the standard Richards Equation 22.3 is a simplification of a more complete multiphase (air–water) approach in that the air phase is assumed to have a negligible effect on variably saturated flow, and that the air pressure varies only little in space and time. This assumption appears adequate for most variably saturated flow problems. Similar assumptions, however, are generally not possible when NAPLs are present. Hence mathematical descriptions of multiphase flow and transport in general require flow equations for each of the three fluid phases, mass transport equations for all organic components (including those associated with the solid phase), and appropriate equations to account for interphase mass transfer processes. We refer readers to reviews by Abriola et al. (1999) and Rathfelder et al. (2000)for discussions of the complexities involved in modeling such systems subject to multiphase flow, multicomponent transport, and interphase mass transfer. An excellent overview of a variety of experimental approaches for measuring the physical and hydraulic properties of multi-fluid systems is given by Lenhard et al. (2002). 

### 22.3.7 Initial and Boundary Conditions 22.3.7.1 Initial Conditions The governing equations for solute transport can be solved analytically or numerically provided that the initial and boundary conditions are specified. Initial conditions need to be specified for each equilibrium phase concentration, that is,

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

where c i is the initial concentration [ML

− 3 ], as well as for all nonequilibrium phases such as concentrations in the immobile region, sorbed concentrations associated with kinetic sites, and initially attached or strained colloid concentrations. 22.3.7.2 Boundary Conditions Complex interactions between the transport domain and its environment often must be considered for the water flow part of the problem being considered since these interactions determine the magnitude of water fluxes across the domain boundaries. By comparison, the solute transport part of most analytical and numerical models usually considers only three types of boundary conditions. When the concentrationJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 22 — #22 22 -22 The Handbook of Groundwater Engineering

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

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

where c 0 is a prescribed concentration [ML

− 3 ] at or along the  D Dirichlet boundary segments. This boundary condition is often referred to as a concentration boundary condition. A third-type (Cauchy type) boundary condition may be used to prescribe the concentration flux at the boundary as follows:

− θ Dij ∂ c ∂ x

j n i + q i n i c = q i n i c 0 for (x , z ) ∈  C (22.57)

in which q i n i represents the outward fluid flux [LT

− 1 ], n i is the outward unit normal vector and c 0 is the concentration of the incoming fluid [ML

− 3 ]. In some cases, for example, when a boundary is impermeable ( q0 =0) or when water flow is directed out of the region, the Cauchy boundary condition reduces to a second-type (Neumann type) boundary condition of the form:

θ Dij ∂ c ∂ x

j n i = 0 for (x , z ) ∈  N (22.58)

Most applications require a Cauchy boundary condition rather than Dirichlet (or concentration) boundary condition. Since Cauchy boundary conditions define the solute flux across a boundary, the solute flux entering the transport domain will be known exactly (as specified). This specified solute flux is then in the transport domain divided into advective and dispersive components. On the other hand, Dirichlet boundary condition controls only the concentration on the boundary, and not the solute flux which, because of its advective and dispersive contributions, will be larger than for a Cauchy boundary condition. The incorrect use of Dirichlet rather than Cauchy boundary conditions may lead to significant mass balance errors at an earlier time, especially for relative short transport domains (van Genuchten and Parker, 1984).

A different type of boundary condition is sometimes used for volatile solutes when they are present in both the liquid and gas phases. This situation requires a third-type boundary condition, but modified to include an additional volatilization term accounting for gaseous diffusion through a stagnant boundary layer of thickness d [L] above the soil surface. The additional solute flux is often assumed to be proportional to the difference in gas concentrations above and below this boundary layer (e.g., Jury et al., 1983). This modified boundary condition has the form:

− θ Dij ∂ c ∂ x

j n i + q i n i c = q i n i c 0 + D g d ( k g c

− g atm ) for ( x , z ) ∈  C (22.59)

where D g is the molecular diffusion coefficient in the gas phase [L2 T
− 1 ] and g atm is the gas concentration above the stagnant boundary layer [ML

− 3 ]. We note that Jury et al. (1983) assumed g atm in Equation 22.59 to be zero. Still other types of boundary conditions can be used. One example is the use of the Bateman equations (Bateman, 1910) to account for a finite rate of release of multiple solutes that are subject to a first-order sequential decay (e.g., radionuclides, nitrogen species). These solutes are typically released from a waste site into the environment as a consequence of decay reactions in the waste site (e.g., van Genuchten, 1985). 

## 22.4 Analytical Models
A large number of computer models using both analytical and numerical solutions of the flow and solute transport equations are now available for a wide range of applications in research and management of natural subsurface systems (Šim  ̊ unek, 2005). Modeling approaches range from relatively simple analyticalJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 23 — #23 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -23

and semi-analytical models, to more complex numerical codes that permit consideration of a large number of simultaneous nonlinear processes. While for certain conditions (e.g., for linear sorption, a concentration independent sink term φ , and a steady flow field) the solute transport equations (e.g., Equation 22.20, Equation 22.21, Equation 22.23, Equation 22.31, and Equation 22.38) become linear, the governing flow equations (e.g., Equation 22.3) are generally highly nonlinear because of the nonlinear dependency of the soil hydraulic properties on the pressure head or water content. Consequently, many analytical solutions have been derived in the past for solute transport equations, and are now widely used for analyzing solute transport during steady-state flow. Although a large number of analytical solutions also exist for the unsaturated flow equation, they generally can be applied only to relatively simple flow problems. The majority of applications for water flow in the vadose zone require a numerical solution of the Richards equation.

Analytical methods are representative of the classical mathematical approach for solving differential equations to produce an exact solution for a particular problem. Analytical models usually lead to anexplicit equation for the concentration (or the pressure head, water content, or temperature) at a particular time and location. Hence, one can evaluate the concentration directly without the time stepping that is typical of numerical methods. While exceptions exist (e.g., Liu et al., 2000; Cotta et al., 2005),analytical solutions usually can be derived only for simplified transport systems involving linearized governing equations, homogeneous soils, simplified geometries of the transport domain, and constant or highly simplified initial and boundary conditions. Unfortunately, analytical solutions for more complex situations, such as for transient water flow or nonequilibrium solute transport with nonlinear reactions, are generally not available and cannot be derived, in which case numerical models must be employed (Šim  ̊ unek, 2005).


### 22.4.1 Analytical Approaches Analytical solutions are usually obtained by applying various transformations (e.g., Laplace, Fourier or other transforms) to the governing equations, invoking a separation of variables, and using the Green’s function approach (e.g., Leij et al., 2000; Cotta et al., 2005).
22.4.1.1 Laplace Transformation The Laplace transform, £, of the solute concentration with respect to time is defined as follows:

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

0 c ( x , t ) exp (

− st ) d t (22.60)

where s is the transform variable [T

− 1 ]. Laplace transforms can greatly simplify the governing solute transport equations by eliminating one independent variable, usually time. They also transform the original transport equation from a partial to an ordinary differential equation. Using Laplace transforms, one thus obtains a governing ADE in the Laplace domain that is much simpler to solve analytically than the original equation. The analytical solution in the Laplace space is subsequently inverted to the real space using either a table of Laplace transforms, by applying inversion theorems, or by using a numerical inversion program.

22.4.1.2 Fourier Transformation Two- and three-dimensional problems often require not only Laplace transforms with respect to time (transformation from t to s) and one spatial coordinate (usually x being transformed to r) , but often also a double Fourier transform with respect to two other spatial coordinates (y and z ) , which is given by (LeijJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 24 — #24 22 -24 The Handbook of Groundwater Engineering

and Toride, 1997):

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

−∞ ∞ 

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

where i 2 = − 1 and γ and κ are transformation variables pertaining to the y and z coordinates. Similarly as for Laplace transforms, Fourier transforms lead to an equation that can be solved analytically more easily. One also needs to transform this analytical solution from the Laplace and Fourier domains back to the original time and space domains. Several analytical solutions, including the derivation process, for different initial and boundary conditions, and different transport domain geometries are given by Leij and Toride (1997) in the N3DADE manual. 22.4.1.3 Method of Moments Statistical moments are often used to characterize the statistical distribution of the solute concentration versus time and space. For example, moments can be used to evaluate the average residence time or the width of the residence time distribution. These attributes can be quantified with the moments of the distribution (Skaggs and Leij, 2002). The pth time moment, m p , is defined as:

m

p = ∞ 

0 tp c ( t ) d t (p = 0, 1, 2, . . .) (22.62)

where c ( t ) is the breakthrough curve measured at a certain location. The zeroth moment is related to the mass of solute contained in the breakthrough curve, while the first moment is related to the mean residence time. Normalized and central moments are defined as:

M

p = m p
m0

M

p = 1

m0 ∞ 

0 ( 1

− M1 ) p c ( t ) d t (22.63)

The second central moment relates to the degree of solute spreading. Similar expressions can be written also for the spatial (depth) moments characterizing the spatial solute distribution. Moments for particular problems or processes can then be obtained by substituting appropriate analytical solutions in the above equations. 

### 22.4.2 Existing ModelsNumerous analytical solutions of the linear advection–dispersion solute transport equations (e.g., Equation 22.21, Equation 22.23, Equation 22.31, and Equation 22.38), or their two- and three-dimensional equivalents, have been developed during the last 40 years and are now widely used for predictive purposesand for analyzing laboratory and field observed concentration distributions. The majority of these solutions pertain to solute transport equations assuming constant water content, θ , and flux, q , values (i.e., for steady-state water flow conditions in a homogeneous medium). 22.4.2.1 One-Dimensional Models Some of the more popular one-dimensional analytical transport models have been CFITM (van Genuchten, 1980b), CFITIM (van Genuchten, 1981), CXTFIT (Parker and van Genuchten, 1984), andCXTFIT2 (Toride et al., 1995). While CFITM considers only one-dimensional equilibrium solute transport in both finite and infinite solute columns, CFITIM additionally considers physical and chemicalJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 25 — #25 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -25

nonequilibrium transport (i.e., the two-region mobile–immobile model for physical nonequilibrium and the two-site sorption model for chemical nonequilibrium). CXTFIT further expanded the capabilities of CFITIM by considering more general initial and boundary conditions, as well as degradation processes. CXTFIT2 (Toride et al., 1995), an updated version of CXTFIT, solves both direct and inverse problems for three different one-dimensional transport models (1) the conventional advection–dispersion equation, ADE; (2) the chemical and physical nonequilibrium ADEs; and (3) a stochastic stream tube model based upon the local-scale equilibrium or nonequilibrium ADE. These three types of models all consider linear adsorption, and include zero- and first-order decay/source terms. All of these models may be used also to solve the inverse problem by fitting the analytical ADE solutions to experimental results.

22.4.2.2 Multi-Dimensional Models Some of the more popular multi-dimensional analytical transport models have been AT123D (Yeh, 1981), 3DADE (Leij and Bradford, 1994), N3DADE (Leij and Toride, 1997), and MYGRT (Ungs et al., 1998). These programs provide analytical solutions to transport problems in two- and three-dimensional domains. 3DADE also included parameter estimation capabilities.

A large number of analytical models for one-, two-, and three-dimensional solute transport problems were recently incorporated into the public domain software package STANMOD (STudio of ANalytical MODels) (Šim  ̊ unek et al., 1999a) (http://www.hydrus2d.com). Figure 22.6 shows a STANMOD dialog window where users can make a selection of the analytical program to be used. This Windows-based computer software package includes not only programs for equilibrium advective-dispersive transport (i.e., the CFITM of van Genuchten 1980b) for one-dimensional transport and 3DADE (Leij and Bradford, 1994) for three-dimensional problems (Figure 22.7), but also programs for more complex problems. For example, STANMOD also incorporates the CFITIM (van Genuchten, 1981) and N3DADE (Leij and Toride, 1997) programs for nonequilibrium transport (i.e., the two-region mobile–immobile model for physical nonequilibrium and the two-site sorption model for chemical nonequilibrium) in one and multiple dimensions, respectively.

FIGURE 22.6 Selection of an analytical model in the STANMOD software package (Šim  ̊ unek et al., 1999a).

## JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 26 — #26 22-26 The Handbook of Groundwater Engineering

FIGURE 22.7 Selection of the geometry of the transport domain for the 3DADE code (Leij and Bradford, 1994) within the STANMOD software package (Šim  ̊ unek et al., 1999a).

Figure 22.7 shows a dialog window of STANMOD in which users can select the geometry of the transport domain for use in the 3DADE program. STANMOD also includes CXTFIT2 (Toride et al., 1995) as discussed earlier, and the CHAIN code of van Genuchten (1985) for analyzing the advective– dispersive transport of up to four solutes involved in sequential first-order decay reactions. Examples of the latter are the migration of radionuclides in which the chain members form first-order decay reactions, and the simultaneous movement of various interacting nitrogen or organic chemicals. The latest version of STANMOD also includes the screening model of Jury et al. (1983) for describing transport and volatilization of soil-applied volatile organic chemicals.

An application of CFITIM within of the STANMOD program is demonstrated in Figure 22.8, which shows measured and fitted breakthrough curves of a nonreactive ( 3 H 2 O) solute for transport through a saturated Glendale clay loam soil. In this example, a tritiated water pulse of 3.102 pore volumes wasapplied to a 30 cm long column, with the breakthrough curve being determined from the effluent concentrations. An analytical solution for two-region (mobile–immobile) physical nonequilibrium transport (van Genuchten, 1981) was used for the analysis. With the pore-water velocity known ( v = 37.54 cm d

− 1 ) and assuming that the retardation factor R for 3 H 2 O was equal to one (no sorption or anion exclusion), only three parameters were optimized against the breakthrough curve, that is, the dispersion coefficient D ( = 15.6 cm 2 d

− 1 ) , a dimensionless variable β ( = 0.823) for partitioning the liquid phase in mobile and immobile regions, and a dimensionless mass transfer coefficient ω ( = 0.870).

An example application of the 3DADE program (Leij and Bradford, 1994) in the STANMOD package is demonstrated in Figure 22.9, which shows calculated concentration contours for a transport problem in which solute is applied at a concentration C 0 = 1, assuming a third-type boundary condition, to a rectangular area of the soil surface (or within groundwater) of 15 × 15 cm 2 . Other transport parameters were as follows: v = 10 cm d

− 1
, D x = 100 and D y = D z = 10 cm 2 d

− 1 , and R = 1. Figure 22.9 shows concentration contours in the xy-plane at t = 1 days, and z = 0.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 27 — #27 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -27

0 0246 0.2 0.4 0.6 0.8

1 1.2

8 Pore volumes [

–]Concentration [–] Measured Fitted

FIGURE 22.8 Breakthrough curve for nonreactive tritium transport as analyzed with the two-region physical nonequilibrium program CXTFIT2 (Toride et al., 1995) within STANMOD (Šim  ̊ unek et al., 1999a).

20

0

y

–20

0.0 0.2 0.4 0.6 0.8 1.0 0 x 30

FIGURE 22.9 Solute distribution in the xy-plane for continuous application at x = 0 and

− 7.5 < y < 7.5 cm after 1 days, with v = 10 cm d

− 1
, D x = 100 and D y = D z = 10 cm 2 d

− 1 as calculated with 3DADE (Leij and Bradford, 1994) within STANMOD (Šim  ̊ unek et al., 1999a).


## 22.5 Numerical Models 22.5.1 Numerical ApproachesWhile analytical and semi-analytical solutions are still popularly used for many relatively simple applications, the ever-increasing power of personal computers, and the development of more accurate andJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 28 — #28 22 -28 The Handbook of Groundwater Engineering

stable numerical solution techniques, have led to the much wider use of numerical models in recent years. Numerical methods are in general superior to analytical methods in terms of being able to solve much more practical problems (Šim  ̊ unek, 2005). They allow users to design complicated geometries that reflect complex natural geologic and hydrologic conditions, control parameters in space and time, prescribe more realistic initial and boundary conditions, and permit the implementation of nonlinear constitutive relationships. Numerical methods usually subdivide the time and spatial coordinates into smaller pieces, such as finite differences, finite elements, and finite volumes, and reformulate the continuous form of governing partial differential equations in terms of a system of algebraic equations. To obtain solutions at certain times, numerical methods generally require intermediate simulations (time-stepping) between the initial condition and the points in time for which the solution is needed. Reviews of the history of development of various numerical techniques used in vadose zone flow and transport models are given by van Genuchten and Šim  ̊ unek (1996) and Šim  ̊ unek (2005). Here we will summarize only the main principles behind some of the more popular techniques. 22.5.1.1 Finite Differences The method of finite differences is generally very intuitive and relatively easy to implement. Time and space are both divided into small increments  t and  z (or  x and  z ) (Figure 22.10). Temporal and spatial derivatives in the governing equation are then replaced with finite differences (formally using Taylor series expansions). For example, the standard advection–dispersion equation for steady-state water flow:

∂ c ∂ t = − ∂ J ∂ z = D ∂ 2 c ∂ z 2 − v ∂ c ∂ z (22.64)
can be approximated as follows using an explicit (forward-in-time) finite difference scheme:

c j + 1 i − c j i  t = − J j i + 1 / 2 − J j i

− 1 / 2  z = D c j i + 1 − 2 c j i + c j i

− 1 (  z ) 2 − v c j i + 1 − c j i

− 1 2  z (22.65)

where subscripts refer to spatial discretization and superscript to temporal discretization (j and j +1 are for the previous and actual time levels, respectively; see Figure 22.10),  t is the time step, and  z is the spatial step (assumed to be constant). Notice that this equation contains only one unknown variable (i.e., the concentration c j + 1 i at the new time level), which hence can be evaluated directly (explicitly) by solving the equation.

z

i +1 c i +1

c

i c i c

i +1

c

i

– 1 c

i

– 1 j +1

j +1

j +1 j

j

j t i +1 ti

z

i

– 1 z

i

z

t ∆ t ∆z ∆z i

∆x ∆x

i

– 1 ∆x i ∆z ∆z j

x z j

z ∆z

i

– 1

z j

– 1 z j +1

∆z j

– 1 x

i

– 1

c

i

– 1, j +1

c

i

– 1, j

c

i

– 1, j

– 1 c i, j c i +1, j +1

c

i +1, j

c

i +1, j

– 1 c

i, j

– 1 c

i, j +1 x i + 1 x

i (a) (b)

FIGURE 22.10 Examples of the spatial and temporal finite difference discretization of a one-dimensional problem (a), and the finite difference discretization of a two-dimensional domain (b).JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 29 — #29 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -29

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

c j + 1 i − c j i  t = − J j + 1 i + 1 / 2 − J j + 1 i

− 1 / 2  z = D c j + 1 i + 1 − 2 c j + 1 i + c j + 1 i

− 1 (  z ) 2 − v c j + 1 i + 1 − c j + 1 i

− 1 2  z (22.66)

and an implicit (weighted) finite difference scheme as:

c j + 1 i − c j i  t = D ε( c j + 1 i + 1 − 2 c j + 1 i + c j + 1 i

− 1 ) + ( 1

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

− 1 ) ( z ) 2

− v ε( c j + 1 i + 1 − c j + 1 i

− 1 ) + ( 1

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

− 1 ) 2  z (22.67)

where ε is a temporal weighting coefficient. Different finite-difference schemes results depending upon the value of ε , that is, an explicit scheme when ε = 0, a Crank–Nicholson time-centered scheme when

ε = 0.5, and a fully implicit scheme when ε = 1. Of these, the latter two cases lead to several unknown concentrations in each equation at the new time level. Consequently, equations for all nodes of the transport domain need to be assembled into an algebraic system of linear equations to produce a matrix equation of the form [ P ] j + 1 { c }j + 1 = { F } (22.68) in which [ P ] is a tridiagonal matrix [ P ] given by:

[ P ] = d1 e 1 0 0 b2 d 2 e 2 0 0 0 b3 d 3 e 3 0 0

.

.

. .

.

. .

.

. 0 0 bN

− 2 d N

− 2 e N

− 2 0 0 0 bN

− 1 d N

− 1 e N

− 1 0 0 bN d N (22.69)

Equation 22.69 is subsequently solved using a matrix equation solver to yield the concentrations at the new time level.

The problem becomes more complex for problems involving transient flow. In that case one must first obtain estimates of the water contents and fluxes at the new time level using a numerical solution of the Richards equation. These water contents and fluxes are subsequently used for assembling the matrix equations for solute transport. While finite differences are relatively straightforward to implement for many transport problems, one major disadvantage is that this method cannot be easily applied to relatively complicated (irregular) transport domains in two and three dimensions. 22.5.1.2 Finite Elements Finite element methods can be implemented in very much the same way as finite differences for one-, two-, and three-dimensional problems. A major advantage of the finite elements is that they are much easier used to discretize complex two- and three-dimensional transport domains (Figure 22.11). As an example, Figure 22.11 shows triangular unstructured finite element grids for a regular rectangular and an irregular domain as generated with the automated MeshGen2D mesh generator of HYDRUS-2D (Šim  ̊ unek et al., 1999b). Note that even though the figure on the right (Figure 22.11) has an irregular soil surface, as well as a tile drain within the transport domain, MeshGen2D could easily discretize/accommodate this transport domain using an unstructured triangular finite element mesh.

The finite element method assumes that the dependent variable in the advection–dispersion equation (e.g., the concentration function c ( x , t ) in (22.64)) can be approximated by a finite series c (x , t ) of theJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 30 — #30 22 -30 The Handbook of Groundwater Engineering

Z

X (a) (b)

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

i -1 i i +11 2

∆x i

x c

c i– 1 c i +1

1

x c

i i

– 1 1

i + 1 x i (a)

(b)

FIGURE 22.12 Example of one-dimensional linear basis functions (a) and their use to approximate solution c ( x , t ) of the advection–dispersion equation (b).

form (Figure 22.12b):

c ( x , t ) = N 

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

where φm are the selected linear basis functions that fulfill the condition φm (xn ) = δ nm , δ nm is Kronecker delta ( δ nm = 1 for m = n , and δ nm = 0 for m = n ), c m are the unknown time-dependent coefficients which represent solutions of (22.64) at the finite element nodal points, and N is the total number of nodal points. For example, for the one-dimensional finite element xi < x < xi + 1 , linear basis functions haveJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 31 — #31 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -31

the following form (Figure 22.12a):

φ 1 = 1

− x

− xi  x φ2 = x

− xi  x 
  

   xi ≤ x ≤ xi + 1 (22.71)

where  x ( = xi + 1 − xi ) is the size of a finite element [L], that is, the distance between two neighboring nodal points. The approximate solution c ( x , t ) converges to the correct solution c ( x , t ) as the number of basis functions N increases.

Application of the Galerkin method, which postulates that the differential operator associated with the transport equation is orthogonal to each of the N basis functions (e.g., Pinder and Gray, 1977), leads to the following system of N time-dependent differential equations with N unknown values cn ( t ) :

L 

0 

− ∂ c ∂ t − D ∂ 2 c ∂ x 2 − v ∂ c ∂ x φ n d x = 0 (22.72)

where φ n are the selected linear weighting functions that can be the same as the basis functions φm . Integrating by parts the terms containing the spatial derivatives to remove the second-order derivatives, and substituting (22.70) for c ( x , t ) leads to the following equation:


e Le



0 

− ∂ cm ∂ t φm φ n dx

− 

e Le



0  Dcm ∂φm ∂ x − vcm φm ∂φ n ∂ x dx

− qsL φ n ( L) + qs 0 φ n ( 0 ) = 0

(22.73)

where qs 0 and qsL are solute fluxes across the left and right boundaries, respectively; e is the element index, Le is the size of an element e , and L is the size of the transport domain (length of the soil profile). The integrals in Equation 22.73 can be evaluated either numerically or analytically, depending upon the type of invoked basis and weighting functions ( φm and φ n , respectively). The equation can be rewritten in matrix form similarly as (Equation 22.68) and then solved using a matrix equation solver. For a detailed derivation of the matrix equation we refer to standard finite element textbooks (e.g., Pinder and Gray, 1977; Zienkiewicz, 1977; Huyakorn and Pinder, 1983) or manuals of some of the models (e.g., the HYDRUS models by Šim  ̊ unek et al., 1998a,b). 22.5.1.3 Stability and OscillationsNumerical solutions of the transport equation often exhibit oscillatory behavior and/or excessive numerical dispersion near relatively sharp concentration fronts (Huyakorn and Pinder, 1983; Šim  ̊ unek et al., 1998a). These problems can be especially serious for advection-dominated transport characterized by small dispersivities. One way to partially circumvent numerical oscillations is to use upstream weighing in which case the flux term of Equation 22.64 is not weighted using regular linear basis functions φ n , but instead with nonlinear quadratic functions φ u n . Undesired oscillations can often be prevented also by selecting an appropriate combination of the space and time discretizations. Two dimensionless numbersmay be used to characterize the discretizations in space and time. One of these is the grid Peclet number, Pe e , which defines the predominant type of solute transport (notably the ratio of the advective and dispersive transport terms) in relation to coarseness of the finite element grid:

Pe e = v  x D (22.74)JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 32 — #32 22 -32 The Handbook of Groundwater Engineering

where  x is the characteristic length of a finite element [L]. The Peclet number increases when the advective part of the transport equation dominates the dispersive part (i.e., when a relatively steep concentration front is present). To achieve acceptable numerical results, the spatial discretization must be kept relatively fine to maintain a low Peclet number. Numerical oscillations can be virtually eliminated when the local Peclet numbers do not exceed about 5. However, acceptably small oscillations may be obtained with local Peclet numbers as high as 10 (Huyakorn and Pinder, 1983). Undesired oscillations for higher Peclet numbers can be effectively eliminated by using upstream weighing.

A second dimensionless number, the Courant number, Cre , may be used to characterize the relative extent of numerical oscillations in the numerical solution. The Courant number is associated with the time discretization,  t [T], as follows:
Cr e = v  t  x (22.75) Given a certain spatial discretization, the time step must be selected such that the Courant number remains less than or equal to 1. Perrochet and Berod (1993) developed a criterion for minimizing or eliminating numerical oscillations that is based on the product of the Peclet and Courant numbers:

Pe

· Cr ≤ ωs ( = 2 ) (22.76)

where ωs is the performance index [

− ]. This criterion indicates that advection-dominated transport problems having large Pe numbers can be safely simulated provided Cr is reduced according to (22.76) (Perrochet and Berod, 1993). When small oscillations in the solution can be tolerated, ωs can be increased to about 5 or 10. 

### 22.5.2 Existing Models 22.5.2.1 Single-Species Solute Transport Models A large number of numerical models are now available for evaluating variably saturated water flow and solute transport processes in the subsurface. Some of these models are in the public domain, such as MACRO (Jarvis, 1994), SWAP (van Dam et al., 1997), UNSATH (Fayer, 2000), VS2DI (Healy, 1990), and HYDRUS-1D (Šim  ̊unek et al., 1998a, 2005), while others are in the commercial domain, such as HYDRUS2D (Šim  ̊ unek et al., 1999b) and MODFLOW-SURFACT (HydroGeoLogic, 1996). These models vary widely in terms of their complexity, sophistication, and ease of use. Although some models are still being run under the DOS operating system, with associated difficulties of preparing input files and interpreting tabulated outputs, many others, especially those in the commercial domain, are supported by sophisticated graphics-based interfaces that greatly simplify their use (Šim  ̊ unek et al., 1998a; 1999b). Several studies have recently reviewed and compared various numerical models for vadose zone applications (e.g., Wilson et al., 1999; Scanlon et al., 2002; MDH Engineered Solutions Corp., 2003; Vanderborght et al., 2005). These studies typically evaluated comparable precision, speed, and ease of use of the codes involved.

While many of the early numerical zone models were developed mostly for agricultural applications (e.g., optimizing soil moisture conditions for crop production), this focus has increasingly shifted toenvironmental research, with the primary concern now being the subsurface fate and transport of various agricultural and other contaminants, such as radionuclides, hydrocarbons, fumigants, pesticides, nutrients, pathogens, pharmaceuticals, viruses, bacteria, colloids, and/or toxic trace elements. The earlier models solved the governing flow and transport equations for relatively simplified system-independent boundary conditions (i.e., specified pressure heads or fluxes, and free drainage). Models developed recently can cope with much more complex system-dependent boundary conditions evaluating surface flow and energy balances and accounting for the simultaneous movement of water, vapor, and heat. Examples are DAISY (Hansen et al., 1990), TOUGH2 (Pruess, 1991), SHAW (Flerchinger et al., 1996), SWAP (van Dam et al., 1997), HYDRUS-1D (Šim  ̊ unek et al., 1998a, 2005), UNSATH (Fayer, 2000), and COUP (JanssonJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 33 — #33 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -33

and Karlberg, 2001). Several models now also account for the extremely nonlinear processes associated with the freezing and thawing cycle (e.g., DAISY, SHAW, and COUP). Models have recently also become increasingly sophisticated in terms of the type and complexity of solute transport processes that can be simulated. Transport models are no longer being limited to solutes undergoing relatively simple chemical reactions such as linear sorption and first-order decay, but now consider also a variety of nonlinear sorption and exchange processes, physical and chemical nonequilibrium transport, volatilization, gas diffusion, colloid attachment/ detachment, decay chain reactions, and many other processes (e.g., the HYDRUS-1D and -2D codes of Šim  ̊ unek et al., 1999b, 2005, or MODFLOW-SURFACT of HydroGeoLogic, Inc., 1996). For example, the general formulation of the transport equations in the HYDRUS codes permit simulations of not only non-adsorbing or linearly sorbing chemicals, but also of a variety of other contaminants, such a viruses (Schijven and Šim  ̊ unek, 2001), colloids (Bradford et al., 2002), cadmium (Seuntjens et al., 2001), and hormones (Casey et al., 2003, 2004), or chemicals involved in the sequential biodegradation of chlorinated aliphatic hydrocarbons (Schaerlaekens et al., 1999; Casey and Šim  ̊ unek, 2001). Much effort has recently been directed also toward improving models for purposes of simulating nonequilibrium and/or preferential flow. Examples of this are the TOUGH codes (Pruess, 1991, 2004), MACRO (Jarvis, 1994), and HYDRUS-1D (Šim  ̊ unek et al., 2003, 2005). These models typically assume the presence of dual-porosity and dual-permeability regions, with different fluxes possible in the two regions. As discussed in detail in Section 22.2b, dual-porosity, mobile–immobile water flow models result when the Richards or kinematic wave equations are used for flow in the fractures, and immobile water is assumed to exist in the matrix, whereas dual-permeability models assume that water is mobile in both the matrix and fracture domains, while invoking terms accounting for the exchange of water and solutes between the two regions. Example applications of these dual-porosity and dual-permeability models are given by Šim  ̊ unek et al. (2001), Haws et al. (2005), Köhne et al. (2004b, 2005), and Pot et al. (2005), among many others.

As an example of available vadose zone flow and transport models, we discuss here in somewhat greater detail the HYDRUS software packages of Šim  ̊ unek et al. (1999b, 2005) and the MODFLOW-SURFACT model developed by HydroGeoLogic (1996).
22.5.2.1.1 The HYDRUS Software Packages HYDRUS-1D (Šim  ̊ unek et al., 2005) and HYDRUS-2D (Šim  ̊ unek et al., 1999b) are software packages (www.hydrus2d.com) that simulate the one- and two-dimensional movement of water, heat, and multiple solutes in variably saturated media, respectively. Both programs use finite elements to numerically solve the Richards equation for saturatedunsaturated water flow and Fickian-based advectiondispersion equations for both heat and solute transport. The flow equation includes a sink term to account for water uptake by plant roots as a function of both water and salinity stress. The flow equation may also consider dual-porosity-type flow with a fraction of water content being mobile, and fraction immobile. The heat transport equation considers conduction as well as advection with flowing water. The solute transport equations assume advective-dispersive transport in the liquid phase, and diffusion in the gaseous phase. The transport equations also include provisions for nonlinear and/or nonequilibrium reactions between the solid and liquid phases, linear equilibrium reactions between the liquid and gaseous phases, zeroorder production, and two firstorder degradation reactions: one which is independent of other solutes, and one which provides the coupling between solutes involved in the sequential first-order decay reactions. In addition, physical nonequilibrium solute transport can be accounted for by assuming a two-region, dual-porosity type formulation which partitions the liquid phase into mobile and immobile regions. Alternatively, the transport equations include provisions for kinetic attachment/detachment of solute to the solid phase and thus can be used to simulate the transport of viruses, colloids, or bacteria (Schijven and Šim  ̊ unek, 2002; Bradford et al., 2003). Both programs may be used to analyze water and solute movement in unsaturated, partially saturated, or fully saturated media. The flow region itself may consist of nonuniform (layered) soils. HYDRUS-2D additionally can handle flow regions delineated by irregular boundaries. The flow region itself may beJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 34 — #34 22 -34 The Handbook of Groundwater Engineering

composed of nonuniform soils having an arbitrary degree of local anisotropy. Flow and transport can occur in the vertical plane, the horizontal plane, or in a three-dimensional region exhibiting radial symmetry about the vertical axis. The unsaturated soil hydraulic properties (the constitutive relationships) can be described using van Genuchten (1980a), Brooks and Corey (1964), Kosugi (1996), and Durner (1994) type analytical functions, or modified van Genuchten type functions that produce a better description of the hydraulic properties near saturation. HYDRUS-1D incorporates hysteresis by assuming that drying scanning curves are scaled from the main drying curve, and wetting scanning curves from the main wetting curve. Root growth is simulated by means of a logistic growth function, while root water uptake can be simulated as a function of both water and salinity stress.

The HYDRUS-1D software package additionally includes modules for simulating carbon dioxide and major ion solute movement (Šim  ̊unek et al., 1996). Diffusion in both the liquid and gas phases and advection in the liquid phase are considered as the main CO 2 transport mechanisms (Šim  ̊ unek and Suarez, 1993). The major variables of the chemical system are Ca, Mg, Na, K, SO 4 , Cl, NO 3 , H 4 SiO 4 , alkalinity, and CO 2 (Šim  ̊ unek and Suarez, 1994; Šim  ̊ unek et al., 1996). The model accounts for equilibriumchemical reactions between these components such as complexation, cation exchange and precipitationdissolution. For the precipitation-dissolution of calcite and dissolution of dolomite, either equilibrium or multicomponent kinetic expressions are used, involving both forward and back reactions. Otherdissolution-precipitation reactions considered include gypsum, hydromagnesite, nesquehonite, and sepiolite. Since the ionic strength of soil solutions can vary considerably with time and space and often reach high values, both modified Debye–Hückel and Pitzer expressions were incorporated into the model as options to calculate single ion activities.

In addition, both HYDRUS programs implement a Marquardt-Levenberg type parameter estimation technique (Marquardt, 1963; Šim  ̊ unek and Hopmans, 2002) for inverse estimation of soil hydraulic (Šim  ̊ unek et al., 1998b,c; Hopmans et al., 2002) and/or solute transport and reaction (Šim  ̊ unek et al., 2002) parameters from measured transient or steady-state flow and/or transport data.

The HYDRUS packages use a Microsoft Windows based Graphics User Interface (GUI) to manage the input data required to run the program, as well as for nodal discretization and editing, parameter allocation, problem execution, and visualization of results. All spatially distributed parameters, such as those for various soil horizons, the root water uptake distribution, and the initial conditions for water, heat and solute movement, are specified in a graphical environment. HYDRUS-2D includes the MeshGen2D mesh generator, specifically designed for variably-saturated subsurface flow transport problems, for defining more general domain geometries, and for discretizing the transport domain into an unstructured finite element mesh. The program offers graphs of the distributions of the pressure head, water content, waterand solute fluxes, root water uptake, temperature, and solute concentrations in the soil profile at preselected times. Also included is a small catalog of unsaturated soil hydraulic properties (Carsel and Parish, 1988), as well as pedotransfer functions based on neural networks (Schaap et al., 2001). Many different applications of both software packages were referenced throughout this book chapter.

22.5.2.1.2 MODFLOW-SURFACT MODFLOW-SURFACT (HydroGeoLogic, 1996) is a fully integrated flow and transport code, based on the MODFLOW groundwater modeling software package of the U.S. Geological Survey (McDonald and Harbaugh, 1988). While MODFLOW (see detailed discussion in Chapter 20) deals only with fully saturated groundwater flow, MODFLOW-SURFACT expands the applicability of the code to unsaturated domains. MODFLOW-SURFACT includes new flow packages that enhance schemes for performing unconfined simulations to rigorously model desaturation and resaturation of aquifers. It provides an option for discretizing the domain using an axi-symmetric geometry for efficient simulation of pumping tests, recovery tests, etc. In addition to normally allowed external stresses in MODFLOW (i.e., constant head, constant flux, areal recharge, evapotranspiration, drains, and streams), MODFLOW-SURFACT provides a rigorous well withdrawal package, unconfined recharge boundary conditions, and seepage face boundary conditions. Additionally, MODFLOW-SURFACT also includes options for adaptive time-stepping andJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 35 — #35 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -35

output control procedures, and a new solution package based on the preconditioned conjugate gradient method.

MODFLOW-SURFACT includes new Analysis of Contaminant Transport (ACT) modules that allow simulations of single-species and multi-component contaminant transport. The program can handle up to 5 chemical species, which may undergo linear or nonlinear retardation, first-order decay and biochemical degradation, as well as generate transformation products. MODFLOW-SURFACT also permits mass partitioning of a single or multicomponent contaminant with diffusive mass movement in the inactive phase, and mass partitioning of a single or multicomponent contaminant from a depleting immobile NAPL phase with advective and dispersive transport in the active phase and diffusive transport in the inactive phase. 22.5.2.2 Geochemical Transport Models Significant efforts have recently been recently carried out also in coupling physical flow and transport models with geochemical models to simulate ever more complex reactions, such as surface complexation, precipitation/dissolution, cation exchange, and biological reactions. Recent reviews of the development of hydrogeochemical transport models involving reactive multiple components are given by Mangold and Tsang (1991), Lichtner (1996) and Steefel and MacQuarrie (1996), Šim  ̊ unek and Valocchi (2002), and Bell and Binning (2004). Most modeling efforts involving multicomponent transport have thus far focused on the saturated zone, where changes in the water flow velocity, temperature and pH are often much more gradual and hence less important than in the unsaturated zone. Consequently, most multicomponent transport models assumed one- or two-dimensional steady-state saturated water flow with a fixed value of the flow velocity, temperature and p H. Only recently have several multicomponent transport models been published which also consider variably saturated flow; these include DYNAMIX (Liu and Narasimhan, 1989), HYDROGEOCHEM (Yeh and Tripathi, 1990), TOUGH-REACT (Pruess, 1991), UNSATCHEM (Šim  ̊ unek and Suarez, 1994; Šim  ̊ unek et al., 1996, 1997), FEHM (Zyvoloski et al., 1997), MULTIFLO (Lichtner and Seth, 1996), OS3D/GIMRT (Steefel and Yabusaki, 1996), HYDROBIOGEOCHEM (Yeh et al., 1998), FLOTRAN (Lichtner, 2000), MIN3P (Mayer et al., 2002), HP1 (Jacques et al., 2002; Jacques and Šim  ̊ unek, 2005), and HYDRUS-1D (Šim  ̊ unek et al., 2005). Geochemical models can be divided into two major groups: those with specific chemistry and general models (Šim  ̊ unek and Valocchi, 2002). Models with specific chemistry are limited in the number of species they can handle, while their application is restricted to problems having a prescribed chemical system. They are, however, much easier to use and computationally can be much more efficient than general models. Typical examples of models with specified chemistry are those simulating the transport of major ions, such as LEACHM (Wagenet and Hutson, 1987), UNSATCHEM (Šim  ̊ unek and Suarez, 1994; Šim  ̊ unek et al., 1996), and HYDRUS-1D (Šim  ̊ unek et al., 2005). Models with generalized chemistry (DYNAMIX, HYDROGEOCHEM, MULTIFLO, FLOTRAN, OS3D/GIMRT, and HP1, all referenced above) provide users with much more freedom in designing a particular chemical system; possible applications of these models are also much wider.

22.5.2.2.1 HP1 HYDRUS-1D was recently coupled with the PHREEQC geochemical code (Parkhurst and Appelo, 1999) to create a new comprehensive simulation tool, HP1 (acronym for HYDRUS1D-PHREEQC) (Jacques et al., 2003; Jacques and Šim  ̊ unek, 2005). The combined code contains modules simulating (1) transientwater flow in variably saturated media, (2) the transport of multiple components, (3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heat transport. HP1 is a significant expansion of the individual HYDRUS-1D and PHREEQC programs by preserving most of their original features and capabilities. The code still uses the Richards equation for simulating variably saturated water flow andadvection–dispersion type equations for heat and solute transport. However, the program can now simulate also a broad range of low-temperature biogeochemical reactions in water, the vadose zone and in ground water systems, including interactions with minerals, gases, exchangers, and sorption surfaces, based on thermodynamic equilibrium, kinetics, or mixed equilibrium-kinetic reactions.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 36 — #36 22 -36 The Handbook of Groundwater Engineering

0.0008

0.0006

0.0004

0.0002

0 0.012

0.01

0.008

0.006

0.004

0.002

0 0 3 6 9 12 15 0 3 6 9 12 15 Time [days] (a) (b)

Time [days] Na Ca

Cd PbZn Cl

AlConcentrations [mol/L]Concentrations [mol/L]

FIGURE 22.13 Major ions (a) and heavy metals (Zn, Pb, and Cd) (b) concentrations in the effluent of an 8-cm soil long column.
Jacques et al. (2003) and Jacques and Šim  ̊ unek (2005) demonstrated the versatility of the HP1 model on several examples such as (a) the transport of heavy metals (Zn 2 + , Pb 2 + , and Cd2 + ) subject to multiple cation exchange reactions, (b) transport with mineral dissolution of amorphous SiO 2 and gibbsite (Al(OH) 3 ) , (c) heavy metal transport in a medium with a pH-dependent cation exchange complex, (d) infiltration of a hyperalkaline solution in a clay sample (this example considers kinetic precipitation– dissolution of kaolinite, illite, quartz, calcite, dolomite, gypsum, hydrotalcite, and sepiolite), (e) long-term transient flow and transport of major cations (Na + , K + , Ca2 + , and Mg2 + ) and heavy metals (Cd2 + , Zn 2 + , and Pb 2 + ) in a soil profile, (f) cadmium leaching in acid sandy soils, (g) radionuclide transport (U and its aqueous complexes), and (h) the fate and subsurface transport of explosives (TNT and its daughter products 2ADNT, 4ADNT, and TAT) (Šim  ̊ unek et al., 2006). 22.5.2.2.2 Leaching of Heavy Metals from a Soil Column As an example application of HP1, Figure 22.13 shows calculated effluent concentrations of heavy metals (Zn 2 + , Pb 2 + , and Cd2 + ) leached from an 8-cm long soil column having an initial solution as defined in Table 22.2, and with its ion-exchange complex in equilibrium with this solution. Heavy metals initially present in the soil column are leached from the exchange complex by major ions (Ca2 + , Mg2 + , and Al 3 + )

. Water was applied to the top of the column at a steady flow rate of 2 cm day

− 1 and having a chemical composition as given in Table 22.2. Details are given by Jacques and Šim  ̊ unek (2005). 

## 22.6 Concluding Remarks
This chapter demonstrates the abundance of models and modeling approaches that are currently available for simulating variably saturated water flow and contaminant transport at various levels of approximation and for different applications. Models range from relatively simple analytical approaches for analyzing solute transport problems during one-dimensional steady-state flow, to sophisticated numerical models for addressing multi-dimensional variably saturated flow and contaminant transport problems at the field scale. One may expect that unsaturated zone flow and transport models will be used increasingly as tools fordeveloping cost-effective, yet technically sound strategies for resource management and pollution remediation and prevention. Improved understanding of the key underlying processes, continued advances innumerical methods, and the introduction of more and more powerful computers make such simulations increasingly practical for many field-scale problems. For example, models can be helpful tools for designing, testing and implementing soil, water and crop management practices that minimize soil and water pollution. Models are equally needed for designing or remediating industrial waste disposal sites

### JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 37 — #37 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -37

TABLE 22.2 Chemical Components (mmol/L) and Species Considered, and Elemental Compositions of Initial and Boundary Solutions Used in the Column Simulation. X Refers to Ion Exchanger
Solutions Components Species Boundary Initial

Al Al3 + , Al(OH) 2 + , Al(OH) + 2 , Al(OH) 0 3 , Al(OH) − 4 0.1 0.5 Br Br − 3.7 11.9 Cl Cl
− (and Cd, Pb, and Zn-species, see below) 10 0 Ca Ca2 + , Ca(OH) + 5 0 K K + , KOH 0 0 2 Na Na + , NaOH 0 0 6 Mg Mg2 + , Mg(OH) + 1 0.75 Cd Cd2 + , Cd(OH) + , Cd(OH) 0 2 , Cd(OH) − 3 , Cd(OH) 2

− 4 , CdCl +
, CdCl2 , CdCl − 3 0 0.9 Pb Pb 2 + , Pb(OH) + , Pb(OH) 0 2 , Pb(OH) − 3 , Pb(OH) 2

− 4 , PbCl + , PbCl0 2 , PbCl − 3 , PbCl2

− 4 0 0.1 Zn Zn2 + , Zn(OH) + , Zn(OH) 0 2 , Zn(OH) − 3 , Zn(OH) 2

− 4 , ZnCl + , ZnCl0 2 , ZnCl − 3 , ZnCl2

− 4 0 0.25 X AlX3 , AlOHX2 , CaX2 , CdX2 , KX, NaX, MgX2 , PbX2 , ZnX2 (mmol) NA 11

and landfills, for predicting contaminant transport from mining wastes, or for long-term stewardship of nuclear waste repositories. The main challenge is to make the models as realistic as possible for the various applications.Continued progress in subsurface flow and transport modeling requires equal advances in both numerical techniques as well as the underlying science. Addressing preferential flow phenomena, and the relatedproblems of subsurface heterogeneity, including the stochastic nature of boundary conditions (precipitation and/or evapotranspiration), will continute to pose formidable challenges. The same is true for improving multicomponent geochemical transport modeling for the vadose zone. For example, numerical algorithms and databases for multicomponent transport models must be extended to higher temperatures and ionic strengths, complex contaminant mixtures (including especially mixed organic and inorganic wastes), multiphase flow, redox disequilibria for low-temperature systems, and coupled physico-chemical systems to account for possible changes in the soil water retention and hydraulic conductivity functions. Better integration is also needed between variably satuarted subsurface and existing larger-scale surface numerical models, which in turn requires further research on such issues as spatial and temporal scaling of hydrological, chemical and biological processes and properties, linking constitutive (soil hydraulic) relationships to measurements scales, preferential flow, and issues of parameter and model uncertainty.

Many science questions related to colloid and colloid-facilitated transport are also still largely unresolved. This is an area of research where our understanding lags far behind current numerical capabilities. Much work is needed to better understand the processes of filtration, straining, size exclusion, colloid– colloid interactions, mobilization of colloids and microorganisms; accumulation at air–water interfaces, interactions between microorganisms and contaminants (including biodegradation), the effects of both physical factors (water content, flow velocity, textural interfaces) and chemical processes (ionic strength, solution composition, pH) on colloid retention and mobilization, and modeling colloid-facilitated transport during conditions of transient flow.

Also, to the best of our knowledge, no models are currentely available that consider all the various processes simultaneously and in their full complexity, including their mutual interactions. That is, no models exist that consider transient preferential flow and transport in structured soils or fractured rocks, while simultaneously considering complex biogeochemical reactions between contaminants, organic 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 the different types of numerical models isJACQ: “4316_c022” — 2006/5/26 — 19:26 — page 38 — #38 22 -38 The Handbook of Groundwater Engineering

needed to address practical problems of contaminant transport (trace elements, radionuclides, organic contaminants) in complex vadose zone environments.



### References Abriola, L.M., Pennell, K.D., Weber, Jr., W.J. Lang, J.R. and Wilkins, M.D., Persistence and interphase mass transfer of organic contaminants in the unsaturated zone: Experimental observations and mathematical modeling. In J.Y. Parlange and J.W. Hopmans (Eds.) Vadose Zone Hydrology: Cutting Across Disciplines. Oxford University Press, New York, pp. 210–234, 1999. Adamczyk, Z., Siwek, B., Zembala, M. and Belouschek, P., Kinetics of localized adsorption of colloid particles. Adv. Coll. Interf. Sci. , 48, 151–280, 1994. Ahuja, L.R. and Hebson, C., Root Zone Water Quality Model. GPSR Tech. Rep. No. 2, USDA, ARS, Fort Collins, CO, 1992. Anderson, M.P., Movement of contaminants in groundwater: Groundwater transport–advection and dispersion. In Groundwater Contamination, Studies in Geophysics, National Academy Press, Washington, D.C. pp. 37–45, 1984. Appelo, C.A.J. and Postma, D., Geochemistry, Groundwater and Pollution. A. A. Balkema. Rotterdam, 1993. Bache, B.W. and Williams, E.G., A phosphate sorption index for soils. J. Soil Sci. 27, 289–296, 1971. Barry, D.A., Modelling contaminant transport in the subsurface: theory and computer programs. In H. Ghadiri and C.W. Rose (Eds.), Modelling Chemical Transport in Soil: Natural and Applied Contaminants, Lewis Publishers, Chelsea, MI, pp. 105–144, 1992. Bateman, H., The solution of a system of differential equations occurring in the theory of radioactive transformation. Proc. Cambridge Philos. Soc. , 15 (pt V) 423–427, 1910. Bear, J., Dynamics of Fluids in Porous Media, Elsevier Science, New York, 1972. Bell, L.S.J. and Binning, P.J., A split operator approach to reactive transport with the forward particle tracking Eulerian–Lagrangian localized adjoint method. Adv. Water Resour . 27, 323–334, 2004. Bethke, C. M., Geochemical Reaction Modeling — Concepts and Applications. Oxford University Press, New York, p. 397, 1996. Biggar, J.W. and Nielsen, D.R., Miscible displacement and leaching phenomena. Agronomy 11, 254–274, 1967. Bradford, S.A., Yates, S.R., Bettahar, M. and Šim  ̊ unek, J., Physical factors affecting the transport and fate of colloids in saturated porous media. Water Resour. Res. 38, 1327, doi:10.1029/2002WR001340, 2002. Bradford, S. A., Šim  ̊ unek, J., Bettahar, M., van Genuchten, M. Th. and Yates, S. R., Modeling colloid attachment, straining, and exclusion in saturated porous media. Environ. Sci. Technol. , 37, 2242– 2250, 2003. Bradford, S. A., Bettahar, M., Šim  ̊ unek, J. and van Genuchten, M. Th., Straining and attachment of colloids in physically heterogeneous porous media. Vadose Zone J . 3, 384–394, 2004. Bresler, E. and Dagan G., Solute dispersion in unsaturated heterogeneous soil at field scale: I Applications. Soil Sci. Soc. Am. J. 43, 467–472, 1979. Brooks, R.H. and Corey, A. T., Hydraulic properties of porous media. Hydrol. Paper No. 3, Colorado State Univ., Fort Collins, CO, 1964. Brusseau, M.L., Nonideal transport of reactive solutes in porous media; Cutting across history and disciplines. In Parlange, M.B. and Hopmans, J.W. (Eds.), Vadose Zone Hydrology: Cutting Across Disciplines, Oxford University Press, New York, pp. 130–154, 1999. Carsel, R.F. and Parrish, R. S., Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24, 755–769, 1988. Casey, F.X.M. and Šim  ̊ unek, J., Inverse analyses of the transport of chlorinated hydrocarbons subject to sequential transformation reactions. J. Environ. Qual. , 30, 1354–1360, 2001. Casey, F.X.M., Larsen, G.L., Hakk, H. and Šim  ̊ unek, J., Fate and transport of 17 β -Estradiol in soil-water systems. Environ. Sci. Technol. 37, 2400–2409, 2003.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 39 — #39 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -39
Casey, F.X.M., Larsen, G.L., Hakk, H. and Šim  ̊ unek, J., Fate and transport of testosterone in agriculturally significant soils. Environ. Sci. Technol. 38, 790–798, 2004. Corapcioglu, M.Y. and Jiang, S., Colloid-facilitated ground-water contaminant transport. Water Resour. Res. 29, 2215–2226, 1993. Corapcioglu, M.Y. and Kim, S., Modeling facilitated contaminant transport by mobile bacteria. Water Resour. Res. 31, 2639–2644, 1995. AQ: Please update Cotta et al.Cotta, R.M., Goldstein, R., Heilbron, P.F.L. and Ungs, M.J. (Eds.), Hybrid Methods in Environmental Transport Phenomena, Center for Analysis and Simulation in Environmental Engineering, POLI/COPPE, Federal University of Rio de Janeiro, Brazil, 2005 (in press). Dagan G. and E. Bresler, Solute dispersion in unsaturated heterogeneous soil at field scale: I Theory. Soil Sci. Soc. Am. J. 43, 461–467, 1979. Dagan G., Flow and Transport in Porous Formations, Springer-Verlag, Berlin, 1989. Durner, W., Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 32, 211–223, 1994. Fayer, M.J., UNSAT-H Version 3.0: Unsaturated Soil Water and Heat Flow Model. Theory, User Manual, and Examples. Pacific Northwest National Laboratory 13249, 2000. Fitter, A.H. and Sutton, S.D., The use of the Freundlich isotherm for soil phosphate sorption data. J. Soil Sci. 26, 241–246, 1975. Flerchinger, G.N., Hanson, C.L. and Wight, J.R., Modeling evapotranspiration and surface energy budgets across a watershed. Water Resour. Res. 32, 2539–2548, 1996. Flury, M., Fluhler, H., Jury, W. and Leuenberger, J., Susceptibility of soils to preferential flow of water: a field study. Water Resour. Res. 30, 1945–1954, 1994. Freundlich, H., Kapillarchemie; eine Darstellung der Chemie der Kolloide und verwandter Gebiete, Akademische Verlagsgellschaft, Leipzig, Germany, 1909. Gelhar, L.W., Mantaglou, A., Welty, C. and Rehfeldt, K.R., A review of Field Scale Physical Solute Transport Processes in Unsaturated and Saturated Porous Media, EPRI Topicl Report EA-4190, Electric Power Research Institute, Palo Alto, CA, 1985. Gerke, H.H. and van Genuchten, M.Th., A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour. Res. 29, 305–319, 1993a. Gerke, H.H. and van Genuchten, M.Th., Evaluation of a first-order water transfer term for variably saturated dual-porosity flow models. Water Resour. Res. 29, 1225–1238, 1993b. Gerke, H.H. and van Genuchten, M.Th., Macroscopic representation of structural geometry for simulating water and solute movement in dual-porosity media. Adv. Water Resour. 19, 343–357, 1996. Germann, P.F., Kinematic wave approach to infiltration and drainage into and from soil macropores. Trans. ASAE 28, 745–749, 1985. Germann, P.F. and Beven, K., Kinematic wave approximation to infiltration into soils with sorbing macropores. Water Resour. Res. 21, 990–996, 1985. Glotfelty, D.E. and Schomburg, C.J., Volatilization of pesticides. In B.L. Sawhney and K. Brown (Eds.), Reactions and Movement of Organic Chemicals in Soils, Special Publ. No. 22, Soil Science Society of America, Madison, WI, pp. 181–208, 1989. Grolimund, D., Borkovec, M., Barmettler, K. and Sticher, H., Colloid-facilitated transport of strongly sorbing contaminants in natural porous media; a laboratory column study. Environ. Sci. Technol., 30, 3118–3123, 1996. Gunary, D., A new adsorption isotherm for phosphate in soil. J. Soil Sci. 21, 72–77, 1970. Gwo, J.P., Jardine, P.M., Wilson, G.V. and Yeh, G.T., A multiple-pore-region concept to modeling mass transfer in subsurface media. J. Hydrol. 164, 217–237, 1995. Hanselman, T.A., Graetz, D.A. and Wilkie, A.C., Manure-borne estrogens as potential environmental contaminants: A review. Environ. Sci. Technol. 37, 5471–5478, 2003. Hansen, S., Jensen, H.E., Nielsen, N.E. and Svendsen, H., DAISY: Soil Plant Atmosphere System Model. NPO Report No. A 10. The National Agency for Environmental Protection,Copenhagen, 272 pp, 1990.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 40 — #40 22 -40 The Handbook of Groundwater Engineering

Haws, N.W., Rao, P.S.C. and Šim  ̊ unek, J., Single-porosity and dual-porosity modeling of water flow and solute transport in subsurface-drained fields using effective field-scale parameters. J. Hydrol. 2005 (in press). AQ: Please update Haws et al. Healy, R.W., Simulation of solute transport in variably saturated porous media with supplemental information on modifications to the U.S. Geological Survey’s computer program VS2D. U.S. Geological Survey, Water-Resources Investigation Report 90-4025: 125, 1990. Hendrickx, J.M.H., Dekker, L.W. and Boersma, O.H., Unstable wetting fronts in water repellent field soils. J. Env. Qual. 22, 109–118, 1993. Hendrickx, J.M.H. and Flury, M., Uniform and preferential flow, Mechanisms in the vadose zone. In Conceptual Models of Flow and Transport in the Fractured Vadose zone. National Research Council, National Academy Press, Washington, D.C., pp. 149–187, 2001. Hopmans, J.W., Šim  ̊ unek, J., Romano, N. and Durner, W., Inverse modeling of transient water flow. In J.H. Dane and G.C. Topp (Eds.), Methods of Soil Analysis, Part 1, Physical Methods, Chapter 3.6.2, Third edition, SSSA, Madison, WI, pp. 963–1008, 2002. Hutson, J.L. and Wagenet, R.J., A multiregion model describing water flow and solute transport in heterogeneous soils. Soil Sci. Soc. Am. J . 59, 743–751, 1995. Huyakorn, P.S. and Pinder, G.F., Computational Methods in Subsurface Flow , Academic Press, London, 1983. HydroGeoLogic, Inc., MODFLOW-SURFACT, A Comprehensive MODFLOW-Based Flow and Transport Simulator, Version 2.1, HydroGeoLogic, Inc., 1996. Jacques, D., Šim  ̊ unek, J., Mallants, D. and van Genuchten, M. Th., Multicomponent transport model for variably-saturated porous media: Application to the transport of heavy metals in soils. In S.M., Hassanizadeh, R.J. Schotting, W.G. Gray, and G.F. Pinder (Eds.) Computational Methods in Water Resources, XIVth International Conference on Computational Methods in Water Resources, June 23–28, 2002, Delft, The Netherlands. Elsevier, Developments in Water Science, 47, 555–562, 2002. Jacques, D., Šim  ̊unek, J., Mallants, D. and van Genuchten, M. Th., The HYDRUS-PHREEQC multicomponent transport model for variably-saturated porous media: Code verification and application, MODFLOW and More 2003: Understanding through Modeling, Conference Proceedings, E. Poeter, Ch. Zheng, M. Hill, and J. Doherty, (Eds.), Int. Ground Water Modeling Center, Colorado School of Mines, pp. 23–27, 2003. Jacques, D. and Šim  ̊ unek, J., Multicomponent — Variable Saturated Transport Model, Description, Manual, Verification and Examples, Waste and Disposal, SCKCEN, BLG-998, Mol, Belgium, 79 pp, 2005. Jansson, P.-E. and Karlberg, L., Coupled heat and mass transfer model for soil-plant-atmosphere systems. Royal Institute of Technology, Dept of Civil and Environmental Engineering, Stockholm 325 pp, 2001. Jarvis, N.J., The MACRO model (Version 3.1), Technical description and sample simulations. Reports and Dissertations 19. Dept. Soil Sci., Swedish Univ. Agric. Sci., Uppsala, Sweden, 51 pp, 1994. Jarvis, N. J., Modeling the impact of preferential flow on nonpoint source pollution. In H.M. Selim and L. Ma (Eds.), Physical Nonequilibrium in Soils: Modeling and Application. Ann Arbor Press, Chelsea, MI, pp. 195–221, 1998. Jiang, S. and Corapcioglu, M.Y., A hybrid equilibrium model of solute transport in porous media in the presence of colloids. Coll. Surf. A: Physicochem. Eng. Aspects 73, 275–286, 1993. Jury, W.A., Simulation of solute transport using a transfer function model. Water Resour. Res. 18, 363–368, 1982. Jury, W.A., Spencer, W.F. and Farmer, W.J., Behavior assessment model for trace organics in soil: I. Description of model. J. Env. Qual. 12, 558–564, 1983. Jury, W.A., Spencer W.F. and Farmer, W.J., Behavior assessment model for trace organics in soil: I. Chemical classification and parameter sensitivity. J. Env. Qual. 13, 567–572, 1984.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 41 — #41 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -41

Jury W.A. and Sposito, G., Field calibration and validation of solute transport models for the unsaturated

zone. Soil Sci. Soc. Am. J. 49, 1331–1341, 1985. Jury, W.A. and Horton, R., Soil Physics, Sixth ed, John Wiley & Sons, Inc., New York, 2004. Kan, A.T. and Tomson, M.B., Ground water transport of hydrophobic organic compounds in the presence of dissolved organic matter. Env. Toxicol. Chem. 9, 253–263, 1990. Kodešová, R., Kozák, J., Šim  ̊ unek, J. and Vacek, O., Field and numerical study of chlorotoluron transport in the soil profile: Comparison of single and dual-permeability model. Plant, Soil Environ. 51, 310–315, 2005. Köhne, J.M., Köhne, S., Mohanty, B.P. and Šim  ̊ unek, J., Inverse mobile–immobile modeling of transport during transient flow: Effect of between-domain transfer and initial water content. Vadose Zone J. 3, 1309–1321, 2004a. Köhne, J.M., B. Mohanty, J. Šim  ̊ unek and H.H. Gerke, Numerical evaluation of a second-order water transfer term for variably saturated dual-permeability models, Water Resour. Res. 40, doi:10.1029/2004WR00385, 2004b.Köhne, J.M. and Mohanty, B.P., Water flow processes in a soil column with a cylindrical macropore: Experiment and hierarchical modeling. Water Resour. Res., 41, W03010, doi:10,1029/2004WR003303, 2005. Köhne, S., Lennartz, B., Köhne, J.M. and Šim  ̊ unek, J., Bromide transport at a tile-drained field site: Experiment, one- and two-dimensional equilibrium and non-equilibrium numerical modeling. J. Hydrol. 2005 (in press). AQ: Please update Köhne et al. Köhne, J.M., Köhne, S. and Šim  ̊ unek, J., Multi-process herbicide transport in structured soil columns: Experiment and model analysis. J. Contamin. Hydrol. (submitted). Kosugi, K., Lognormal distribution model for unsaturated soil hydraulic properties, Water Resour. Res., 32, 2697–2703, 1996. Kung, K.-J.S., Preferential flow in a sandy vadose zone. 2. Mechanism and implications. Geoderma 46, 59–71, 1990. Lai, Sung-Ho and Jurinak, J.J., Numerical approximation of cation exchange in miscible displacement through soil columns. Soil Sci. Soc. Am. Proc. 35, 894–899, 1971. Langmuir, I., The adsorption of gases on plane surfaces of glass, mica, and platinum. J. Am. Chem. Soc. 40, 1361–1403, 1918. Lapidus, L. and Amundson, N.R., Mathematics of adsorption in beds: VI. The effect of longitudinal diffusion in ion exchange and chromatographic columns. J. Phys. Chem. 56, 984–988, 1952. Leij, F.J. and Bradford, S.A., 3DADE: A computer program for evaluating three-dimensional equilibrium solute transport in porous media. Research Report No. 134, U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1994. Leij, F.J. and Toride, N., N3DADE: A computer program for evaluating nonequilibrium three-dimensional equilibrium solute transport in porous media. Research Report No. 143, U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1997. Leij, F.J., Priesack, E. and Schaap, M.G., Solute transport modeled with Green’s functions with application to persistent solute sources. J. Contam. Hydrol. , 41, 155–173, 2000. Lenhard, R.J., Oostrom, M. and Dane, J.H., Multi-fluid flow. In J.H. Dane and G.C. Topp (Eds.), Methods of Soil Analysis, Part 4, Physical Methods, SSSA Book Series 5, Soil Sci. Soc, Am., Madison, WI, pp. 1537–1619, 2002.Lichtner, P.C., Continuum formulation of multicomponent-multiphase reactive transport. In P.C. Lichtner, C.I. Steefel and E.H. Oelkers (Eds.), Reactive Transport in Porous Media, Reviews in Mineralogy. Chapter 1, Vol. 34. Mineralogical Society of America, pp. 1–81, 1996. Lichtner, P.C. and Seth, M.S., User’s manual for MULTIFLO: Part II MULTIFLOW 1.0 and GEM 1.0. Multicomponent-multiphase reactive transport model. CNWRA 96-010, Center for Nuclear Waste Regulatory Analyses, San Antonio, TX, 1996. Lichtner, P.C., FLOWTRAN User’s manual. Los Alamos National Laboratory Document. NM, 2000.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 42 — #42 22 -42 The Handbook of Groundwater Engineering

Lindstrom, F.T., Haque, R., Freed, V.H. and Boersma, L., Theory on the movement of some herbicides in soils. Environ. Sci. Technol. 1, 561–565, 1967. Lindstrom, F.T., Boersma, L. and Stockard, D., A theory on the mass transport of previously distributed chemicals in a water saturated sorbing porous medium. Isothermal cases, Soil Sci. 112, 291–300, 1971.Lindqvist, R. and Enfield, C.G., Biosorption of dichlorodiphenyltri-chloroethane and hexachlorobenzene in groundwater and its implications for facilitated transport. Appl. Environ. Microbiol. 58, 2211–2218, 1992. Liu, C.W. and Narasimhan, T.N., Redox-controlled multiple-species reactive chemical transport. 1. Model development. Water Resour. Res. 25:869–882, 1989. Liu, C., Szecsody, J.E., Zachara, J.M. and Ball, W.P., Use of the generalized integral transform method for solving equations of solute transport in porous media. Adv. Water. Resour. 23, 483–492, 2000. Logan, B.E., Jewett, D.G., Arnold, R.G., Bouwer, E.J. and O’Melia, C.R., Clarification of clean-bed filtration models. J. Environ. Eng. 121, 869–873, 1995. Magee, B.R., Lion, L.W. and Lemley, A.T., Transport of dissolved organic macromolecules and their effect on the transport of phenanthrene in porous media. Environ. Sci. Technol. 25, 323–331, 1991. Mangold, D.C. and Tsang, A summary of subsurface hydrological and hydrochemical models. Rev. Geophys. 29, 51–79, 1991. Mansfeldt, T., Leyer, H., Barmettler, K. and Kretzschmar, R., Cyanide leaching from soil developed from coking plant purifier waste as influenced by citrate. Vadose Zone J . 3, 471–479, 2004. Marquardt, D.W., An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math. 11, 431–441, 1963.Mayer, U., Frind, E.O. and Blowes, D.W., Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water Resour. Res. 38, 1174, 2002. McDonald, M.G. and Harbaugh, A.W., A modular three-dimensional finite-difference ground-water flow model, U.S. Geological Survey Techniques of Water-Resources Investigations Book 6, Chapter A1, 586 p, 1988. MDH Engineered Solutions Corp., Evaluation of Computer Models for Predicting the Fate and Transport of Salt in Soil and Groundwater, Phase I Report. Pub. No: T/403, ISBN No. 0-7785-2493-0 (Printed Edition), ISBN No. 0-7785-2494-9 (On-line Edition), 2003. Millington, R.J. and J.M. Quirk, Permeability of porous solids. Trans. Faraday Soc. 57, 1200–1207, 1961. Mills W.M., Liu S. and Fong, F.K., Literature review and model (COMET) for colloid/metals transport in porous media. Groundwater 29, 199–208, 1991. Moldrup, O., Poulsen, T.G., Schjonning, P., Olesen, T. and Yamaguchi, T., Gas permeability in undisturbed soils: measurements and predictive models. Soil Sci. 163, 180–189, 1998. Narasimhan, T.N., Buckingham, 1907; An appreciation. Vadose Zone J . 4, 434–441, 2005, Nielsen, D.R., van Genuchten, M. Th. and Biggar, J.W., Water flow and solute transport processes in the unsaturated zone. Water Resour. Res. 22, 89S–108S, 1986. Noell, A.L., Thompson, J.L., Corapcioglu, M.Y. and Triay, I.R., The role of silica colloids on facilitated cesium transport through glass bead columns and modeling. J. Contam. Hydrol. 31, 23–56, 1998. Parker, J.C. and van Genuchten, M.Th., Determining transport parameters from laboratory and field tracer experiments. Bull. 84-3, Va Agric. Exp. St., Blacksburg, VA, 1984. Parkhurst, D.L. and Appelo, C.A.J., User’s guide to PHREEQC (version 2) — A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Water Resources Investigation, Report 99-4259, Denver, Co, USA, 312 pp, 1999. Perrochet, P. and Berod, D., Stability of the standard Crank–Nicolson–Galerkin scheme applied to the diffusion-convection equation: some new insights. Water Resour. Res. 29, 3291–3297, 1993. Pinder, G.F. and Gray, W.G., Finite Element Simulation in Surface and Subsurface Hydrology , Academic Press, New York, 1977.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 43 — #43 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -43

Pitzer, K.S., Activity coefficients in electrolyte solutions. In R.M. Pitkowicz, Activity Coefficients in Electrolyte Solutions. CRC Press, Inc., Boca Raton, FL, pp. 157–208, 1979. AQ: Please update Pot, V., et al. Pot, V., Šim  ̊ unek, J., Benoit, P., Coquet, Y., Yra, A. and Martínez-Cordón, M.-J., Impact of rainfall intensity on the transport of two herbicides in undisturbed grassed filter strip soil cores. J. Contam. Hydrol. 2005 (in press). Pruess, K., Tough2 — A general-purpose numerical simulator for multiphase fluid and heat flow. Rep. LBL-29400, Lawrence Berkeley Laboratory, Berkeley, CA, 1991. Pruess, K., The TOUGH codes — A family of simulation tools for multiphase flow and transport processes in permeable media. Vadose Zone J . 3, 738–746, 2004. Pruess, K. and Wang, J.S.Y., Numerical modeling of isothermal and non-isothermal flow in unsaturated fractured rock — a review. In D.D. Evans and T.J. Nicholson (Eds.), Flow and Transport Through Unsaturated Fractured Rock, Geophysics Monograph 42, American Geophysical Union, Washington, D.C., pp. 11–22, 1987. Rajagopalan, R. and Tien, C., Trajectory analysis of deep-bed filtration with the sphere-in-cell porous media model. AIChE. J . 22, 523–533, 1976. Rathfelder, K.M., Lang, J.R. and Abriola, L.M., A numerical model (MISER) for the simulation of coupled physical, chemical and biological processes in soil vapor extraction and bioventing systems. J. Contam. Hydrol. 43, 239–270, 2000. Richards, L.A., Capillary conduction of fluid through porous mediums. Physics 1, 318–333, 1931. Saiers, J.E., Hornberger, G.M. and Harvey, C., Colloidal silicakaolinite in the transport of cesium through laboratory sand columns. Water Resour. Res. 32, 33–41, 1996. Saiers, J.E. and Hornberger, G.M., The role of colloidal kaolinite in the transport of cesium through laboratory sand columns. Water Resour. Res. 32, 33–41, 1996. Scanlon, B.R., Christman, M., Reedy, R., Porro, I., Šim  ̊ unek, J. and Flerchinger, G., Intercode comparisons for simulating water balance of surficial sediments in semiarid regions. Water Resour. Res. 38, 1323, doi:10.1029/2001WR001233, 59.1–59.16, 2002. Schaerlaekens, J., Mallants, D., Šim  ̊ unek, J. van Genuchten, M.Th. and Feyen, J. Numerical simulation of transport and sequential biodegradation of chlorinated aliphatic hydrocarbons using CHAIN_2D. J. Hydrologi. Proc. 13, 2847–2859, 1999. Schaap, M.G., Leij, F.J. and van Genuchten, M.Th., Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251, 163–176, 2001. Schijven, J.F. and Hassanizadeh, S.M., Removal of viruses by soil passage: Overview of modeling, processes and parameters. Crit. Rev. Environ. Sci.Technol. 30, 49–127, 2000. Schijven, J. and Šim  ̊ unek, J., Kinetic modeling of virus transport at field scale. J. Contam. Hydrol. 55, 113–135, 2002. Selim, H.M., Davidson, J.M. and Mansell, R.S., Evaluation of a two-site adsorption-desorption model for describing solute transport in soils, paper presented at Proceedings, Summer Computer Simulation Conference, Nat. Sci. Found. , Washington, DC, 1976. Seuntjens, P., Tirez, K., Šim  ̊ unek, J., van Genuchten, M.Th., Cornelis, C. and Geuzens, P., Aging effects on cadmium transport in undisturbed contaminated sandy soil columns. J. Environ. Qual. 30, 1040–1050, 2001. Shapiro, R.E. and Fried M., Relative release and retentiveness of soil phosphates. Soil Sci. Soc. Am. Proc. 23, 195–198, 1959. Sibbesen, E., Some new equations to describe phosphate sorption by soils. J. Soil Sci. 32, 67–74, 1981. Šim  ̊ unek, J. and Suarez, D.L., Modeling of carbon dioxide transport and production in soil: 1. Model

development. Water Resour. Res. 29, 487–497, 1993. Šim  ̊ unek, J. and Suarez, D.L., Two-dimensional transport model for variably saturated porous media with

major ion chemistry. Water Resour. Res. 30, 1115–1133, 1994. Šim  ̊ unek, J. and van Genuchten, M.Th., The CHAIN_2D code for simulating two-dimensional movement

of water flow, heat, and multiple solutes in variably-saturated porous media, Version 1.1, Research Report No 136 , U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 205 pp, 1994.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 44 — #44 22 -44 The Handbook of Groundwater Engineering

Šim  ̊unek, J., Suarez, D.L. and Šejna, M., The UNSATCHEM software package for simulating onedimensional variably saturated water flow, heat transport, carbon dioxide production and transport, and multicomponent solute transport with major ion equilibrium and kinetic chemistry. Version 2.0. Research Report No. 141, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 186 pp, 1996. Šim  ̊ unek, J. and Suarez, D.L., Sodic soil reclamation using multicomponent transport modeling. ASCE J.

Irrig. Drain. Eng. 123, 367–376, 1997. Šim  ̊ unek, J., Šejna, M. and van Genuchten, M.Th. The HYDRUS-1D software package for simulating

the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 2.0, IGWMC-TPS-70 , International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 202pp, 1998a. Šim  ̊ unek, J., Wendroth, O. and van Genuchten, M.Th., A parameter estimation analysis of the evaporation

method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62, 894–905, 1998b. Šim  ̊ unek, J., van Genuchten, M.Th., Gribb, M.M. and Hopmans, J.W., Parameter estimation of unsaturated

soil hydraulic properties from transient flow processes. Soil Tillage Res. 47, Special issue "State of the art in soil physics and in soil technology of antrophic soils", 27–36, 1998c. Šim  ̊ unek, J., van Genuchten, M.Th., Šejna, M., Toride, N. and Leij, F.J., The STANMOD computer software

for evaluating solute transport in porous media using analytical solutions of convection-dispersion equation. Versions 1.0 and 2.0, IGWMC-TPS-71 , International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 32 pp, 1999a. Šim  ̊unek, J., Šejna, M. and van Genuchten, M.Th., The HYDRUS-2D software package for simulating twodimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0, IGWMC-TPS-53, International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 251pp, 1999b. Šim  ̊ unek, J., Wendroth, O., Wypler, N. and van Genuchten, M.Th., Nonequilibrium water flow

characterized from an upward infiltration experiment. Eur. J. Soil Sci. 52, 13–24, 2001. Šim  ̊ unek, J. and Hopmans, J.W., Parameter optimization and nonlinear fitting. In: J.H. Dane and

G.C. Topp, (Eds.) Methods of Soil Analysis, Part 1, Physical Methods, Chapter 1.7, 3rd ed., SSSA, Madison, WI, pp. 139–157, 2002. Šim  ̊ unek, J., Jacques, D., Hopmans, J.W., Inoue, M., Flury, M. and van Genuchten, M.Th., Solute transport

during variably-saturated flow — inverse methods. In J.H. Dane and G.C. Topp, (Eds.), Methods of Soil Analysis, Part 1, Physical Methods, Chapter 6.6, 3rd ed., SSSA, Madison, WI, pp. 1435–1449, 2002. Šim  ̊ unek, J. and Valocchi, A.J., Geochemical transport, In J. H. Dane and G.C. Topp, (Eds.), Methods of

Soil Analysis, Part 1, Physical Methods, Chapter 6.9, 3rd ed., SSSA, Madison, WI, pp. 1511–1536, 2002. Šim  ̊ unek, J., Jarvis, N.J., van Genuchten, M.Th. and Gärdenäs, A., Nonequilibrium and preferential flow

and transport in the vadose zone: review and case study. J. Hydrol. 272, 14–35, 2003. Šim  ̊ unek, J., van Genuchten, M.Th. and Šejna, M., The HYDRUS-1D software package for simulating

the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 3.0, HYDRUS Software Series 1 , Department of Environmental Sciences, University of California Riverside, Riverside, CA, 270pp, 2005. AQ: Please update Šim  ̊ unek Šim  ̊ unek, J., Models of water flow and solute transport in the unsaturated zone. In: M.G. Anderson (Ed.),

The Encyclopedia of Hydrological Sciences, John Wiley & Sons, Ltd., Chichester, England, 25 pp, 2005 (in press). AQ: Please update Šim  ̊ unek et al. Šim  ̊ unek, J., Jacques, D., van Genuchten, M.Th. and Mallants, D., Multicomponent geochemical transport

modeling using the HYDRUS computer software packages. J. Am. Water Resour. Assoc. 2006 (in press). Sips, R., On the structure of the catalyst surface, II. J. Chem. Phys. 18, 1024–1026, 1950.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 45 — #45 Contaminant Transport in the Unsaturated Zone Theory and Modeling 22 -45

Skaggs, T.H. and Leij, F.J., Solute transport: Theoretical background. In J.H. Dane and G.C. Topp, (Eds.), Methods of Soil Analysis, Part 1, Physical Methods, Chapter 6.3, 3rd ed., SSSA, Madison, WI, pp. 1353–1380, 2002. Sposito, G. and Barry, D.A., On the Dagan model of solute transport I groundwater: Foundational aspects. Water Resour. Res. 23, 1867–1875, 1987.Steefel, C.I. and Yabusaki, S.B., OS3D/GIMRT, software for multicomponent–multidimensional reactive transport. User Manual and Programmer’s Guide, PNL-11166, Pacific Northwest National Laboratory, Richland, WA, 1996. Steefel, C.I. and MacQuarrie, K.T.B., Approaches to modeling of reactive transport in porous media. In P.C. Lichtner, C.I. Steefel and E.H. Oelkers (Eds.), Reactive Transport in Porous Media. Reviews in Mineralogy. Chapter 2, Vol. 34, Mineralogical Society of America, pp. 83–129, 1996. Thiele-Bruhn, S., Pharmaceutical antibiotic compounds in soils — a review. J. Plant Nutr. Soil Sci. 166, 145–167, 2003. Tolls, J., Sorption of veterinary pharmaceuticals in soils: A review. Environ. Sci. Technol. 35, 3397–3406. 2001.Toride, N., Leij, F.J. and van Genuchten, M.Th., A comprehensive set of analytical solutions for nonequilibrium solute transport with first-order decay and zero-order production. Water Resour. Res. 29, 2167–2182, 1993. Toride, N., Leij, F.J. and van Genuchten, M.Th., The CXTFIT code for estimating transport parameters from laboratory or field tracer experiments. Version 2.0, Research Report No. 137 , U. S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1995. Truesdell, A.H. and Jones, B.F., Wateq, a computer program for calculating chemical equilibria of natural waters. J. Res. U.S. Geol. Surv. 2, 233–248, 1974. Ungs, M., Wilkinson, K. and Summers, K., Users guide for MYGRT, Version 3.0: Software for simulating migration of organic and inorganic chemicals in grioundwater, EPRI Report TR-111748, Electric Power Research Institute, Palo Alto, CA, 1998. van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabat, P. van Valsum, P.E.V. and Groenendijk, P. and van Diepen, C.A., Theory of SWAP, version 2.0. Simulation of water flow, solute transport and plant growth in the Soil–Water–Atmosphere–Plant environment. Dept. Water Resources, WAU, Report 71, DLO Winand Staring Centre, Wageningen, Technical Document 45, 1997. Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiery, D.,Vanclooster, M., Mouvet, C. andVereecken, H., A set of analytical benchmarks to test numerical models of flow and transport in soils. Vadose Zone J. 4, 206–221, 2005. van de Weerd, H., Leijnse, A. and van Riemsdijk, W.H., Transport of reactive colloids and contaminants in groundwater: effect of nonlinear kinetic reactions. J. Contam. Hydrol. 32, 313–331, 1998. van Genuchten, M.Th. and Wierenga, P.J., Mass transfer studies in sorbing porous media, I. Analytical solutions. Soil Sci. Soc. Am. J . 40, 473–481, 1976. van Genuchten, M.Th., A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898, 1980a. van Genuchten, M.Th. and J. Šim  ̊ unek, Evaluation of pollutant transport in the unsaturated zone, In: P.E. Rijtema and V. Eliáš (Eds.), Regional Approaches to Water Pollution in the Environment , NATO ASI Series: 2. Environment, Kluwer, Dordrecht, The Netherlands, pp. 139–172, 1996. van Genuchten, M.Th., Davidson, J.M. and Wierenga, P.J., An evaluation of kinetic and equilibrium equations for the prediction of pesticide movement through porous media. Soil Sci. Soc. Am. Proc. 38, 29–35, 1974. van Genuchten, M.Th., Determining transport parameters from solute displacement experiments. Research Report No. 118, U. S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1980b. van Genuchten, M.Th., Non-equilibrium transport parameters from miscible displacement experiments. Research Report No. 119, U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1981.JACQ: “4316_c022” — 2006/5/26 — 19:26 — page 46 — #46 22 -46 The Handbook of Groundwater Engineering

van Genuchten, M.Th. and Wagenet, R.J., Two-site/two-region models for pesticide transport and degradation: Theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53, 1303–1310, 1989. van Genuchten, M.Th. and Parker, J.C., Boundary conditions for displacement experiments through short laboratory soil columns. Soil Sci. Soc. Am. J. 48, 703–708, 1984. van Genuchten, M.Th., Convective–dispersive transport of solutes involved in sequential first-order decay reactions. Comp. & Geosci. 11, 129–147, 1985. van Genuchten, M.Th. and Šim  ̊ unek, J., Integrated modeling of vadose zone flow and transport processes. Proc. Unsaturated Zone Modelling: Progress, Challenges and Applications, Eds. R. A. Feddes, G. H. de Rooij, and J. C. van Dam, Wageningen, The Netherlands, October 3–5, 2004, pp. 37–69, 2004. Vinten, A.J.A., Yaron, B. and Nye, P.H., Vertical transport of pesticides into soil when adsorbed on suspended particles. J. Agric. Food Chem. 31, 662–664, 1983. Von Gunten, H.R., Weber, V.E. and Krahenbuhl, U., The reactor accident at Chernobyl: A possibility to test colloid-controlled transport of radionuclides in a shallow aquifer. J. Contam. Hydrol. 2, 237–247, 1998. Wagenet, R.J. and Hutson, J.L., LEACHM: Leaching Estimation And CH emistry M odel, A process-based model of water and solute movement, transformations, plant uptake and chemical reactions in the unsaturated zone. Continuum 2. Water Resour. Inst., Cornell University, Ithaca, New York, 1987. Westall, J.C., Zachary, J.L. and Morel, F.M.M. MINEQL: A computer program for the calculation of chemical equilibrium composition of aqueous systems. Tech. Note 18, 91pp, Dep. of Civ. Eng., Mass. Inst. of Technol., Cambridge, MA, 1976. White, N. and Zelazny, L.W., Charge properties in soil colloids. In D.L. Sparks (Ed.) Soil Physical Chemistry. CRC Press, Boca Raten, Florida, pp. 39–81, 1986. Wheatcraft, S.W. and Cushman, J.H., Hierarchical approaches to transport in heterogeneous porous media. Rev. Geophys. 29 (Suppl.), 261–267, 1991. Wilson, G.V., Albright, W.H., Gee, G.W., Fayer, M.J. and Ogan, B.D., Alternative cover assessment project, Phase 1 Report. Final Technical Report prepared by Desert Research Institute, NV and Pacific Northwest National Laboratory, WA for U.S. Environmental Protection Agency under Contract No. 44-0000-2038, 202p, 1999. Yeh, G.T., Analytical transient one-, two-, and three-dimensional simulation of waste transport in the aquifer system, ORNL- 5602, Oak Ridge National Laboratory, Oak Ridge, TN, 1981. Yeh, G.T. and Tripathi, V.S., HYDROGEOCHEM: A coupled model of hydrological and geochemical equilibrium of multicomponent systems. Rep. ORNL-6371, Oak Ridge Natl. Lab., Oak Ridge, TN, 1990. Yeh, G.T., Salvage, K.M., Gwo, J.P., Zachara, J.M. and Szecsody, J.E., HYDROBIOGEOCHEM: A coupled model of hydrological transport and mixed biochemical kinetic/equilibrium reactions in saturated– unsaturated media. Rep. ORNL/TM-13668, Oak Ridge Natl. Lab., Oak Ridge, TN, 1998. Zhang, P., Šim  ̊ unek, J. and Bowman, R.S., Nonideal transport of solute and colloidal tracers through reactive zeolite/iron pellets. Water Resour. Res. 40, doi:10.1029/2003WR002445, 2004. Zienkiewicz, O.C., The Finite Element Method , 3rd ed., McGraw-Hill, London, 1977. Zyvoloski, G.A., Robinson, B.A., Dash, Z.V. and Trease, L.L., Summary of the models and methods for the FEHM applicatio n — a finite element heat- and mass-transfer code. Los Alamos National Laboratory Rept. LA-13307-MS, NM, 1997.