78: Models of Water Flow and Solute Transport
in the Unsaturated Zone
JIRKA ŠIMŮNEK
Department of Environmental Sciences, University of California, Riverside, CA, US

A large number of models for simulating water flow and solute transport in the unsaturated zone are now
increasingly being used for a wide range of applications in both research and management. Modeling approaches
range from relatively simple analytical and semianalytical solutions, to complex numerical codes. While
analytical and semianalytical solutions are still popular for some applications, the ever-increasing power of
personal computers and the development of more accurate and numerically stable solution techniques havemotivated much wider use of numerical codes in recent decades. The wide use of numerical models is alsosignificantly enhanced by their availability in both the public and commercial domains, and by the developmentof sophisticated graphics-based interfaces that can tremendously simplify their use.
In this paper I focus mainly on numerical models, give a brief history of their development, and discuss someof the more often used numerical techniques including relatively efficient matrix solvers that now are available
for multidimensional models. Names and web addresses of some of the more popularly used numerical codes
simulating vadose zone processes are also provided. Finally, some typical problems in which the numericalcodes have been applied are identified.

INTRODUCTION
Many models of varying degree of complexity and dimensionality have been developed during the past several
decades to quantify the basic physical and chemical processes affecting water flow and pollutant transport in theunsaturated zone. Computer models based on analytical andnumerical solutions of the flow and solute transport equations are now increasingly being used for a wide range ofapplications in research and management of natural subsurface systems. Modeling approaches range from relatively
simple analytical and semianalytical models to more complex numerical codes that permit consideration of a large
number of simultaneous nonlinear processes. Whereas analytical and semianalytical solutions are still more popular
for most relatively simple applications, the ever-increasing
power of personal computers and the development of moreaccurate and numerically stable solution techniques havegiven rise to the much wider use of numerical models inrecent decades. The wide use of numerical models has also
been significantly enhanced by their availability in bothEncyclopedia of Hydrological Sciences. Edited by M G Anderson.
 2005 John Wiley & Sons, Ltd.

public and commercial domains, and by the developmentof sophisticated graphics-based interfaces that tremendously
simplify their use.
Analytical, semianalytical, and numerical models areusually based on the following three governing equations
for water flow, solute transport, and heat movement,
respectively:∂θ (h)
∂
∂h
=
K(h)
+1 −S
∂t
∂z
∂z

 

∂
∂θ Rc
∂c
=
θD
− qc − φ
∂t
∂z
∂z


 
∂
∂T
∂C(θ )T=λ(θ )− Cw qT
∂t
∂z
∂z

(1)
(2)
(3)

Suitable simplifications (mostly for analytical approaches)or extensions thereof (e.g. for two- and three-dimensional
systems) are also employed. In equation (1), often referred
to as the Richards equation, z is the vertical coordinate
positive upwards, t is time, h is the pressure head, θ is

1172

SOILS

the water content, S is a sink term representing root water
uptake or some other sources or sinks, and K(h) is theunsaturated hydraulic conductivity function, often given as
the product of the relative hydraulic conductivity, Kr , andthe saturated hydraulic conductivity, Ks . In equation (2),known as the convection-dispersion equation (CDE), c isthe solution concentration, R is the retardation factor that
accounts for adsorption, D is the dispersion coefficientaccounting for both molecular diffusion and hydrodynamic
dispersion, q is the volumetric fluid flux density, and φ isa sink/source term that accounts for various zero- and firstorder or other reactions. In equation (3), T is temperature, λis the apparent thermal conductivity, and C and Cw are thevolumetric heat capacities of the soil and the liquid phase,respectively. Solutions of the Richards equation (1) require
knowledge of the unsaturated soil hydraulic functions, thatis, the soil water retention curve, θ (h), describing therelationship between the water content θ and the pressure
head h, and the unsaturated hydraulic conductivity function,
K(h), defining the hydraulic conductivity K as a functionof h or θ . While under certain conditions (i.e. for linear
sorption, a concentration-independent sink term φ, and asteady flow field) equations (2) and (3) are linear equations,
equation (1) is generally highly nonlinear because of thenonlinearity of the soil hydraulic properties. Consequently,
many analytical solutions have been derived in the past
for equations (2) and (3) and these analytical solutions are
now widely used for analyzing solute and heat transport
under steady-state conditions. Although a large numberof analytical solutions of (1) exist, they can generally
be applied only to drastically simplified problems. The
majority of applications for water flow in the vadose zone
requires a numerical solution of the Richards equation.

ANALYTICAL MODELS
Analytical methods represent a classical mathematical
approach to solve differential equations, leading to an exact
solution for a particular problem. Analytical solutions areusually obtained by applying various transformations (e.g.
Laplace or Fourier transformations) to the governing equations, invoking a separation of variables, or the Green’s
function approach (e.g. Leij et al., 2000). Analytical models usually result in an explicit equation that, for example,
states that concentration (or the pressure head, water content, or temperature) is equal to a certain value at particulartime and location. One can therefore evaluate a particularvariable directly without time-stepping, typical of numerical methods. Many analytical solutions lead to relatively
complicated formulations that include infinite series and/or
integrals that need to be evaluated numerically, which suggests some ambiguity in the often-claimed advantage ofexactness of analytical methods over numerical techniques.
On the other hand, using analytical solutions one can often

more easily evaluate interrelationships among parameters,
and get better insight into how various processes control
the basic flow and transport processes (e.g. using dimensionless variables and parameters). Analytical solutions areoften also used to check the correctness and accuracy ofnumerical models, although numerical models can equallywell be used to check the correctness of the complex analytical solutions.
Analytical solutions can usually be derived only forsimplified transport systems involving linearized governing equations, homogeneous soils, simplified geometriesof the transport domain, and constant or highly simplified
initial and boundary conditions. Unfortunately, analyticalsolutions for more complex situations, such as for transient water flow or nonequilibrium solute transport withnonlinear reactions, are generally not available and/or cannot be derived, in which case numerical models must
be employed.
Solute Transport

Numerous analytical solutions of equations (2) and (3), or
their two- and three-dimensional equivalents, have beendeveloped in the last 40 years and are now widely used
for predictive purposes and/or analyzing laboratory and/or
field observed concentration distributions. The majority ofthese solutions pertains to equations (2) and (3) assuming
constant water content, θ , and flux, q, values (i.e. for steadystate water flow conditions). Since equation (2) and (3)have the same form, analytical solutions derived for solutetransport can often also be used immediately for many heattransport problems, and vice versa.
Some of the more popular one- and multidimensional
analytical transport models have been CXTFIT (Parker
and van Genuchten, 1984), AT123D (Yeh, 1981), and
3DADE (Leij and Bradford, 1994). A large number ofanalytical models for one-, two-, and three-dimensional
solute transport problems were recently incorporated into
the comprehensive software package STANMOD (STudio of ANalytical MODels) (Šimůnek et al., 1999b)(http://www.ussl.ars.usda.gov/models/stanmod/
stanmod.HTM). This Windows-based computer software
package includes not only programs for evaluating analytical solutions for equilibrium convective-dispersive solutetransport (i.e. the CFITM of van Genuchten (1980) for onedimensional transport and 3DADE for three-dimensional
problems), but also additional programs that solve more
complex problems. For example, it also incorporates the
CFITIM (van Genuchten, 1981) and N3DADE (Leij andToride, 1997) programs for nonequilibrium convectivedispersive transport (i.e. the two-region mobile-immobile
model for physical nonequilibrium and the two-site sorption model for chemical nonequilibrium) for one- andmultidimensional formulations, respectively. STANMOD
also includes CXTFIT2 (Toride et al., 1995), an updated

MODELS OF WATER FLOW AND SOLUTE TRANSPORT IN THE UNSATURATED ZONE

version of CXTFIT to solve both direct and inverse problems for three different one-dimensional transport models:(i) the conventional CDE (1); (ii) the chemical and physical nonequilibrium CDEs; and (iii) a stochastic stream tube
model based upon the local scale equilibrium or nonequilibrium CDE. All three models consider linear adsorption,
as well as the zero- and first-order decay/source terms.
In addition, STANMOD includes the CHAIN code of van
Genuchten (1985) for analyzing the convective-dispersive
transport of up to four solutes involved in sequential
first-order decay reactions. Examples are the migration ofradionuclides, in which the chain members form first-orderdecay reactions, and the simultaneous movement of various
interacting nitrogen or organic chemicals. The latest version
of STANMOD also includes the screening model of Juryet al. (1983) for describing transport and volatilization ofsoil-applied volatile organic chemicals.
Water Flow

The highly nonlinear Richards equation can be solved
analytically only for a very limited number of cases involving homogeneous soils, simplified initial and boundary
conditions, and relatively simple constitutive relationships
describing the unsaturated soil hydraulic properties. The
Richards equation for this purpose needs to be first linearized. This can be accomplished using several mathematical transformations, with the most common transformation
being the Kirchhoff integral transformation:= hK(h) dh(4)−∞

Linearization of the Richards equation often also involves
using an exponential law relating hydraulic conductivity
and the pressure head (Gardner, 1958):
K(h) = Ks exp(αh)(5)

where α is an empirical coefficient.
Development of analytical and semianalytical solutionsof the unsaturated flow equations has been geared mostly
towards infiltration problems (Wooding, 1968; Philip, 1969,1992). While some analytical solutions are widely used, forexample, for evaluating tension disc experiments (Wooding,1968), designing irrigation systems (Philip, 1992), or studying flow around buried objects (Warrick and Knight, 2002),
many others appear to be only of academic interest.

NUMERICAL METHODS
Numerical methods are superior to analytical methods interms of being able to solve practical problems. They allow
users to design complicated geometries that reflect complex

1173

natural geologic and hydrologic conditions, control parameters in space and time, prescribe realistic initial andboundary conditions, and to implement nonlinear constitutive relationships. Numerical methods usually subdivide
the time and spatial coordinates into smaller pieces, such as
finite differences, finite elements, and/or finite volumes, andreformulate the continuous form of governing partial differential equations in terms of a system of algebraic equations.
In order to obtain solutions at certain times, numericalmethods generally require intermediate simulations (timestepping) between the initial condition and the points intime for which the solution is needed. The following two
sections review the history of development of various
numerical techniques used in vadose zone flow and transport models. The review is based in part on an earlier review
by van Genuchten and Šimůnek (1996). After reviewing
various numerical techniques I will also discuss the need for
efficient matrix solvers for 2D and 3D models, and provide
a list of some of the most often used flow/transport models.
Numerical Solution of Richards Equation

A variety of numerical methods may be used to solve
the variably saturated flow equation. The popularity ofnumerical methods stems from the fact that the Richards
equation can be solved analytically only for a very limited
number of cases. Even so, the highly nonlinear nature of
the Richards equation also hampered the development ofcomputationally efficient numerical methods that are stable
under all conditions, particularly for infiltration in very dry
soils. Stable numerical solutions still require relatively finediscretizations of both the time and space domains, often
resulting in excessive CPU and simulation times, especiallyfor two- and three-dimensional problems and/or problems
involving highly transient boundary conditions.
Early applications of numerical methods for solving
variably saturated flow problems generally used classical finite differences. Integrated finite difference, finite
volumes, and especially, finite element methods (Neuman, 1973; Huyakorn et al., 1986; Šimůnek et al., 1999a)
became increasingly popular in the seventies and thereafter. While finite difference methods today are used in a
majority of one-dimensional models, finite volume methods
and/or finite element methods coupled with mass lumping of the mass balance term are usually used in twoand three-dimensional models. Finite element methods usedwith unstructured triangular and tetrahedral elements allow
for a more precise description of complex transport domains
compared to finite differences.
The Richards equation can be formulated in three different ways. A mixed formulation arises when both the watercontent and the pressure head variables appear simultaneously in the governing equation, such as in equation (1).
The h-based formulation is obtained when the time derivative of the water content (∂θ/∂t) on the left side of

1174

SOILS

(1) is rewritten using the soil water capacity C as follows:C∂h/∂t, where C is defined as the slope of the retention
curve, that is, dθ /dh. The θ -based formulation is obtained
when the product of the hydraulic conductivity K(h) andthe pressure head gradient (∂h/∂z) on the right side of (1) is
replaced with the product of the water diffusivity Dw (θ )and the water content gradient, ∂θ/∂z. The θ -formulation
allows for very efficient numerical solutions, even for infiltration into initially dry soils. The use of this formulation is,
however, straightforward only for homogeneous and unsaturated soils. This is because, the driving force for water
movement is not the water content gradient, but the pressure
head gradient; water in heterogeneous soils hence does notnecessarily flow from locations with a higher water content
to locations with a lower water content. Special provisions
must be taken for such heterogeneous systems (Hills et al.,1989), a likely reason why θ -formulations are rarely usedin numerical models.
Celia et al. (1990) suggested that numerical solutions
based on the standard h-based formulation of the Richards
equation often yield poor results, characterized by relatively large mass balance errors and incorrect predictionsof the pressure head distributions in the soil profile. They
solved the mixed formulation of the Richards equation
using a modified Picard iteration scheme that possesses
mass-conserving properties for both finite element andfinite difference spatial approximations. Milly (1985) earlier presented two mass-conservative schemes for computing nodal values of the water capacity in the h-based
formulation to force global mass balance. Several highly
efficient numerical schemes based on different types ofpressure head transformations were presented recently byHills et al. (1989), Ross (1990), and Kirkland et al. (1992).
Hills et al. (1989) showed that the θ -based form of the
Richards equation can yield fast and accurate solutions
for infiltration into very dry heterogeneous soil profiles.
However, the θ -based numerical scheme cannot be used
for soils having saturated regions. Kirkland et al. (1992)
expanded the work of Hills by combining the θ -based andh-based models to yield a transformation method applicable also to variably saturated systems. Their approach
involved a new dependent variable, being a linear function of the pressure head and the water content in thesaturated and unsaturated zones, respectively. Additionaltransformations of the Richards equation, with the common goal of decreasing its nonlinearity and increasing
the efficiency of the numerical solution, were reviewed
by Williams et al. (2000). Many of these transformations
were, however, used mostly by the authors themselves andhave not reached wide acceptance. Most popularly usedvadose zone flow models presently utilize the mixed formulation of Celia et al. (1990); these models include SWAP
(van Dam et al., 1997) and HYDRUS (Šimůnek et al.,1998, 1999a).

Time and space discretization of the Richards equation
generally leads to a nonlinear system of algebraic equations.
These equations are most often linearized and solved using
the Newton–Raphson or Picard iteration methods. Picard
iteration is widely used in unsaturated zone models becauseof its ease of implementation, and because this method
preserves symmetry of the final system of equations. The
Newton–Raphson iteration procedure is more complex andresults in nonsymmetric matrices, but often achieves afaster rate of convergence and can be more robust than
Picard iteration for certain types of problems (Paniconi and
Putti, 1994).
Numerical Solution of the Transport Equation

A large number of methods are also available to numerically solve the convection-dispersion solute transport
equation. These methods may be broadly classified into
three groups: (i) Eulerian, (ii) Lagrangian, and (iii) mixed
Lagrangian–Eulerian methods. In the Eulerian approach,the transport equation is discretized by means of a usual
finite difference or finite element method using a fixed
grid system. For the Lagrangian approach, the mesh moves
along with the flow or remains fixed in a deforming
coordinate system. A two-step procedure is followed fora mixed Lagrangian-Eulerian approach. First, convective
transport is considered using a Lagrangian approach inwhich Lagrangian concentrations are estimated from particle trajectories. Subsequently, all other processes including
sinks and sources are modeled with an Eulerian approach
using any finite element or finite differences method, leading to the final concentrations.
Standard finite difference and Galerkin-type finite element methods belong to the first group of Eulerian methods.
Finite differences and finite elements methods provided theearly tools for solving solute transport problems and still
are the most popular methods being used. Numerical experiments have shown that both methods give good results fortransport in which dispersion is a relatively dominant process (e.g. as indicated by the grid Peclet number, vx/D,
where v is pore-water velocity and x the grid size). However, both methods can lead to significant numerical oscillations and/or artificial dispersion for convection-dominatedtransport problems. Still, a relatively simple method may be
used to prevent or limit numerical oscillations. By selecting an appropriate combination of relatively small space
and time steps, it is possible to virtually eliminate all oscillations. Perrochet and Berod (1993) developed criteria forminimizing or eliminating numerical oscillations based ona “performance index”. They concluded that all oscillations
should be eliminated when the performance index, definedas the product of the local Peclet and Courant (vt/x)
numbers, is less than 2. When small oscillations in the solution can be tolerated, the performance index can be easily
increased to about 5 or 10 (Perrochet and Berod, 1993).

MODELS OF WATER FLOW AND SOLUTE TRANSPORT IN THE UNSATURATED ZONE

One alternative is the use of upwind finite difference
methods. This method virtually eliminates numerical oscillations, even for purely convective transport. The disadvantage of this method is that it may create significant
and often unacceptable numerical dispersion. Similarly,
upstream weighting has been proposed for finite elements.
The method consists of using different weighting functions
for terms having spatial derivatives than for other terms in
the transport equation. This approach places greater weighton the upstream nodes within a particular element.
While Lagrangian methods (or method of characteristics)
virtually eliminate problems with numerical oscillations,they may introduce other problems, notably artificial dispersion and nonconservative solutions. Lagrangian methods
are also relatively difficult to implement in two and three
dimensions. Instabilities resulting from inappropriate spatial
discretizations may occur during longer simulations becauseof deformation of the stream function. Furthermore, nonrealistic distortions of the results may occur when modelingthe transport of solutes that are subjected to different sorption/exchange or precipitation reactions.
The Eulerian and Lagrangian approaches can also be
combined to the mixed Eulerian–Lagrangian method.
Because of the different mathematical nature of the diffusive (parabolic) and convective (hyperbolic) terms in
the convection-dispersion equation, the transport equation is best decomposed into a mixed problem consisting of a purely convective problem, followed by a pure
diffusion-only problem. Methods based on this approach arecalled operator-splitting or splitting-up methods. Convective transport is then solved with the Lagrangian approach,while all other terms of the transport equation are solved
using Eulerian methods.
Still other solutions exist, such as the use of a so-called“random walk” process or a combination of analytical andnumerical techniques. In a random walk approach, solutetransport is modeled using a large number of particles.
Displacement of each particle during each time step isgiven by a certain distance, this being the sum of two
velocity contributions – a deterministic and stochastic contribution. Studies with this type of method indicate thatit may be necessary to use many thousands of particles
in order to obtain relatively precise results. An example
of a combination of analytical and numerical techniques
is given by Sudicky (1989) who modeled solute transport
using Laplace transforms with respect to time, and Galerkinfinite elements for the spatial domain. The use of Laplacetransforms avoids the need for intermediate simulations(time-stepping) between the initial condition and the points
in time for which solutions are needed, while also less stringent requirements are needed for the spatial discretization.
This combination of analytical and numerical techniques,unfortunately, has however, one important limitation. Since

1175

Laplace transforms eliminate time as an independent variable in the governing transport equation, all coefficientssuch as water content, flow velocity, and retardation factors,
must be independent of time. This means that combination
methods can only solve solute transport problems during
steady-state water flow, and hence are inappropriate fortransient variably saturated flow situations typical of mostfield problems.
Many of the above methods for numerically solving
transport equation were developed primarily for saturated
conditions, for which coarse spatial discretizations andlarge flow velocities often produce large Peclet numbers,
which may lead to significant numerical oscillations. Flowvelocities in the vadose zone are usually much smaller.
Also, the nonlinearity of the Richards equation generally
requires numerical solutions on much finer spatial grids thanin groundwater studies. Consequently, the Peclet numbers
are significantly smaller for vadose zone applications thanin groundwater flow studies, thereby allowing the adoptionof more oscillation-prone methods. The majority of vadosezone models therefore can use relatively standard finite
element or finite difference methods, which, if combinedwith upwind or upstream weighting or proper self-adjustingor other time step scheme, should eliminate most or all
numerical oscillations.
Matrix Solvers

Discretization and subsequent linearization (as needed) of
the governing partial differential equations for water flow
and solute transport leads to a system of linear equations[A]{x} = {b}(6)in which {x } is an unknown solution vector, {b} is the
known right-hand-side vector of the matrix equation, and[A] is a sparse banded matrix, which is symmetric forwater flow if the modified Picard procedure is used,
but asymmetric for water flow if the Newton–Raphson
method is used. Matrix [A] is generally asymmetric forsolute transport, unless convection is not considered in theformulation. Matrix [A] is tridiagonal for one-dimensional
applications and thus can be solved very quickly andefficiently using simple Gaussian elimination. For higherdimensional application, matrix [A] is a sparse matrix
with a number of lines equal to the number of nodes inthe spatial discretization scheme. The number of nodes
is on the order of tens of thousands for a typical twodimensional application, and on the order of millions fora three-dimensional application. Since this matrix needs to
be inverted many times during a typical numerical run (ateach time step for solute transport, and additionally also ateach iteration for unsaturated water flow), the need for very
efficient solvers cannot be underestimated.

1176

SOILS

Traditionally, the matrix equations have been solved
using direct methods, such as Gaussian elimination or LU
decomposition. These methods usually take advantage ofthe banded nature of the coefficient matrices and, in thecase of water flow, of the symmetric properties of thematrix. Direct solution methods have several disadvantages
as compared to iterative methods. For example, they requirea fixed number of operations (depending upon the size ofthe matrix), which increases approximately by the square of
the number of nodes. Iterative methods, on the other hand,require a variable number of repeated steps, with the number increasing at a much smaller rate (about 1.5) with thesize of a problem (Mendoza et al., 1991). A similar reduction also holds for the memory requirement since iterative
methods do not require one to store nonzero matrix elements. Memory requirements, therefore, increase at a muchsmaller rate with the size of the problem when iterativesolvers are used. This memory requirement is associatedwith the need to minimize the size of the band (i.e. thelargest distance between two neighboring nodal numbers)of matrix [A] for direct methods. While minimization of
the matrix band is trivial for finite difference methods, it
can be rather complex when unstructured triangular finite
element meshes are used. Round-off errors also represent
less of a problem for iterative methods as compared withdirect methods. This is because round-off errors in iterative methods are self-correcting. Finally, for time-dependentproblems, a reasonable approximation of the solution (i.e.
the solution at the previous time step) exists for iterative
methods, but not for direct methods. In general, direct methods are more appropriate for relatively small problems andfor finite difference codes, while iterative methods are moresuitable for larger problems and codes using unstructured
finite element grids.
While many iterative methods have been used in thepast for handling large sparse matrix equations, a variety of increasingly powerful preconditioned accelerated
iterative methods, such as the preconditioned conjugate
gradient method (PCG), are now becoming available also.
Since the system of linear equations resulting from discretization of the solute transport equation is nonsymmetrical, it is necessary to either formulate the transport
problem in such a way that it leads to a symmetricmatrix, or to use an extension of PCG for nonsymmetrical matrices, such as ORTHOMIN (generalized conjugate
residual method) (Mendoza et al., 1991), GMRES (generalized minimal residual method), biconjugate gradients,
TFQMR (transpose-free quasi-minimal residual algorithm),
CGSTAB (conjugate gradient stabilized method), and conjugate gradient squared procedures. Both the preconditioned
conjugate gradient and ORTHOMIN methods consist of two
essential parts: initial preconditioning, and iterative solution
with either conjugate gradient, CGSTAB, or ORTHOMIN

acceleration (Mendoza et al., 1991). Incomplete lowerupper (ILU) factorization can be used as preconditioning ofmatrix [A], which is then factorized into lower and upper
triangular matrices by partial Gaussian elimination. The
preconditioned matrix is subsequently inverted repeatedlyusing updated estimates to provide a new approximation of
the solution.
Available Model for Unsaturated Zone

Most of the early models developed for studying processes
in the near-surface environment focused mainly on variably saturated water flow. They were used primarily inagricultural research for the purpose of optimizing moisture conditions to increase crop production. This focus
has increasingly shifted to environmental research, withthe primary concern now being the subsurface fate andtransport of various agricultural and other contaminants,
such as pesticides, nutrients, pathogens, pharmaceuticals,
viruses, bacteria, colloids, toxic trace elements, and/or fumigants, and also the evaluation of water recharge throughthe vadose zone. While the earlier models solved thegoverning equations (1) through (3) 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 systemdependent boundary conditions evaluating surface flow andenergy balances and accounting for the simultaneous movement of water, vapor, and heat. Examples are DAISY(Hansen et al., 1990), TOUGH2 (Pruess, 1991), SHAW
(Flerchinger et al., 1996), SWAP (van Dam et al., 1997),
HYDRUS-1D (Šimůnek et al., 1998), UNSATH (Fayer,2000), and COUP (Jansson and Karlberg, 2001). Several
models now account also 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 solutetransport 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 andchemical nonequilibrium transport, volatilization, gas diffusion, colloid attachment/detachment, decay chain reactions,and many other processes (e.g. the HYDRUS-1D and 2D codes of Šimůnek et al., 1998, 1999a). For example,
the general formulation of the transport equations in the
HYDRUS codes allows one to simulate not only nonadsorbing or linearly sorbing chemicals but also a variety ofother contaminants, such as viruses (Schijven and Šimůnek,2002), colloids (Bradford et al., 2002), cadmium (Seuntjens
et al., 2001), and hormones (Casey et al., 2003), or chemicals involved in the sequential biodegradation of chlorinatedMODELS OF WATER FLOW AND SOLUTE TRANSPORT IN THE UNSATURATED ZONE

aliphatic hydrocarbons (Schaerlaekens et al., 1999; Caseyand Šimůnek, 2001).
Options to simulate carbon and nitrogen cycles are also
becoming a standard feature of many environmental models, such as DAISY, LEACHN (Hutson and Wagenet,1992), RZWQM (Ahuja and Hebson, 1992), and COUP.
These models typically distribute organic matter, carbon,
and organic and mineral nitrogen over multiple computational pools, while allowing organic matter to be decomposed by multiple microbial biomass populations. They canaccount for most of the major reaction pathways, such as
mineralization-immobilization of crop residues, manure andother organic wastes, mineralization of the soil humus fractions, interpool transfer of carbon and nitrogen, nitrification(ammonium to nitrate-N), denitrification (leading to the production of N2 and N2 O), volatilization loss of ammonia(NH3 ), production and consumption of methane (CH4 ) andcarbon dioxide (CO2 ), changes in the carbon nitrogen ratioof organic matter, and microbial biomass growth and death(Ahuja and Hebson, 1992).
Efforts are also on to couple physical flow and transport models with geochemical models to simulate evenmore complex reactions, such as surface complexation,
precipitation/dissolution, cation exchange, and/or biological reactions (e.g. Ahuja and Hebson, 1992; Šimůnek andSuarez, 1994; Šimůnek and Valocchi, 2002; Jacques et al.,2002). Models considering these chemical reactions, including the ability to simulate the transport of multiple chemical
species and carbon dioxide, are required for studying water
management practices and irrigation techniques under arid
and semiarid conditions, evaluation of water suitability for
irrigation, and reclamation of sodic soils (Šimůnek andSuarez, 1997).
Another active area of research involves attempts toextend existing models that simulate uniform flow tosituations where nonequilibrium and/or preferential flow
occurs. Examples of this are the MACRO (Jarvis, 1994)
and HYDRUS-1D (Šimůnek et al., 2003) models. Possible approaches for simulating preferential flow differ interms of their underlying assumptions and complexity.
They range from relatively simplistic models to more complex, physically based, dual-porosity, dual-permeability,
and multiregion-type models. A relatively simple dualporosity flow model results when the Richards equationis combined with composite (double-hump type) equations
for the hydraulic properties to account for both soil textural(matrix) and soil structural (fractures, macropores, peds)
effects on flow. A more complex dual-porosity, mobileimmobile water flow model results when the Richards orkinematic wave equations are used for flow in the fractures,and immobile water is assumed to exist in the matrix. Even
more complex are various dual-permeability models suchas the formulations of Gerke and van Genuchten (1993)
and Pruess (1991), or the kinematic wave approach as used

1177

in the MACRO model of Jarvis (1994). These formulations all assume that water is mobile in both the matrix
and fracture domains, while invoking terms that accountfor the exchange of water and solutes between the matrix
and the fractures.
A large number of models are now available for simulating processes in the vadose zone. Some of these models arein the public domain, such as MACRO, SWAP, UNSATH(Fayer, 2000), VS2DI (Healy, 1990), and HYDRUS-1D(Šimůnek et al., 1998), while others are in the commercial domain, such as HYDRUS-2D (Šimůnek et al., 1999a).
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, withassociated difficulties of preparing input files and interpreting tabulated outputs, many others, especially thosein the commercial domain, are supported by sophisticatedgraphics-based interfaces that tremendously simplify theiruse (Šimůnek et al., 1998, 1999a).
Table 1 gives a summary of some of the more widelyused numerical models for simulating variably saturated
water flow and solute transport in soils. This table also
provides Internet addresses and references where additional
information about each model can be found. With the
exception of HYDRUS-2D, TOUGH2, and VS2DTI, allmodels given in Table 1 are one-dimensional models,perhaps reflecting the fact that the majority of applications
for unsaturated zone models is still only one dimensional.

CONCLUSIONS
Much of the research in the field of soil science has focusedin recent decades upon understanding the fundamentals
of variably saturated water flow and pollutant fate andtransport processes. As society continues its rapid technological development, the types of pollution problems
and chemicals posing significant environmental threats,
have become increasingly complex. Problems such as thetransport of pesticides, colloids, bacteria, viruses, pharmaceuticals, reproductive hormones, nutrients, and toxic
trace elements, carbon sequestration, and bioremediation oforganic contaminants, all require a thorough understanding
and coupling of multiple hydrogeological, geochemical, andmicrobiological processes. It is the continually increasing
speed and power of modern computers that will enable suchmodels to become convenient tools for analysis of complex
geochemical systems. Although more and more complex
models are being constantly developed, currently available
models are still relatively specialized and no single model ispresently available that can describe the multiple problems
and chemicals mentioned above. Development of numericalmodels capable of describing unstable and/or preferential
flow, as well as models coupled with sophisticated geochemical models capable of describing both instantaneous

1178

SOILS

Table 1 Some of the widely used numerical models for simulating variably saturated water flow and solute transport insoils
Model name
COUP
DAISY
HYDRUS-1DHYDRUS-2D

MACRO
RZWQM
SHAW
SWAP
SWIM
TOUGH2UNSATH
VS2DTI

Internet address, reference, interesting/special features
http://www.lwr.kth.se/Vara%20 Datorprogram/CoupModel/ (Jansson and Karlberg, 2001)
Carbon/nitrogen cycle, thawing/freezing cycle, coupled water, vapor, and heat transport, crop
growthhttp://www.dina.dk/∼daisy/ (Hansen et al., 1990) Carbon/nitrogen cycle, crop growth,management practices, pesticide processes
http://www.hydrus2d.com (Šimůnek et al., 1998) Multiple soil hydraulic functions, neuralnetwork-based pedotransfer functions, nonlinear nonequilibrium solute transport, mobile-immobile
and two-site sorption concepts, chain reactions, volatilization, inverse option, intuitive sophisticatedgraphical interface
http://www.hydrus2d.com (Šimůnek et al., 1999a) Two-dimensional, multiple soil hydraulicfunctions, neural network-based pedotransfer functions, nonlinear nonequilibrium solute transport,
mobile-immobile and two-site sorption concepts, chain reactions, inverse option, unstructured
triangular finite element meshes, intuitive sophisticated graphical interface
http://www.mv.slu.se/BGF/Macrohtm/macro.htm (Jarvis, 1994) Preferential flow using
kinematic wave equation, snow accumulation, pesticide transport
http://gpsr.ars.usda.gov/products/rzwqm.htm (Ahuja and Hebson, 1992) Complex modular
program, crop growth, chemical equilibrium module, management practices, pesticide processes
http://www.nwrc.ars.usda.gov/models/shaw (Flerchinger et al., 1996) Thawing/freezing cycle,
coupled water, vapor, and heat transport, multispecies plant canopy
http://www.swap.alterra.nl/ (van Dam et al., 1997) A three-level drainage system at regional
scale, crop growth
http://www.clw.csiro.au/products/swim (Verburg et al., 1996) Bypass flow, flexible descriptionof hydraulic properties, hyperbolic sine transformation of the pressure head
http://www-esd.lbl.gov/TOUGH2/ (Pruess, 1991) Multidimensional multiphase fluid and heat
flow, dual permeability
http://hydrology.pnl.gov/resources/unsath/unsath−download.asp (Fayer, 2000) Coupledwater, vapor, and heat transport, no solute transport
http://water.usgs.gov/software/vs2di.html (Healy, 1990) Two-dimensional, finite differences

and kinetic chemical and biological reactions will undoubtedly remain a focus of research in the near future.
The accuracy of the obtained predictions depend toa large extent upon the accuracy of available model
input parameters and upon proper conceptualization of soil
heterogeneity and other system complexities, such as thepossible presence of nonequilibrium flow and transport,
including preferential flow. Processes are often described
and their parameters measured on a much smaller scale thanthose for which the model predictions are being sought.
Consequently, many model parameters often need to becalibrated so that they reflect the bulk behavior of theheterogeneous system, in which case they can be used
for larger scale predictions. New measuring techniques thatprovide model parameters on the scale at which predictions
are made are badly needed for successful applications of
unsaturated flow and transport models in a predictive modeat the larger scale.
One may expect that unsaturated zone flow and transportmodels will be used increasingly for integrating fundamental knowledge about the vadose zone to yield tools
for developing cost-effective, yet technically sound strategies for resource management and pollution remediationand prevention. Unsaturated zone transport models areindispensable tools for analyzing complex environmental

pollution problems, and for developing practical management strategies. Models can help guide field observations byidentifying which parameters and processes control systembehavior. Following Steefel and Van Cappellen (1998) andŠimůnek and Valocchi (2002), several specific key ways inwhich unsaturated flow and transport models can be used
are identified below:
1.

2.

3.

Physical, chemical, and biological processes are often
studied in isolation either in the laboratory or in
the field under controlled conditions. Mathematicalmodels can be used to investigate the impacts ofmultiple coupled biogeochemical reactions and otherinteractions in the presence of complex flow fields
and spatial heterogeneity. These models also enable
extrapolation to environmentally relevant temporal andspatial scales.
Numerical transport models provide a useful tool forinterpreting experimental results. Models can help
understand qualitative and quantitative trends andrelationships present in the data. Properly applied modeling to interpret results of field experiments can lead tomore effective quantitative understanding of underlying
biogeochemical processes.
One of the most powerful applications of numericalflow and transport models is conducting sensitivity

MODELS OF WATER FLOW AND SOLUTE TRANSPORT IN THE UNSATURATED ZONE

4.

analyses. Such analyses permit a systematic evaluationof the impact of model parameters (physical, chemical, and/or biological), initial conditions, and boundary
conditions upon the model output. The results of asensitivity analysis provide insight into the relative
importance of individual processes and reactions within
a complex biogeochemical system. Results can help
one identify the most important parameters and processes, and thereby provide guidance in allocation of
resources for laboratory and field investigations.
Numerical flow and transport models are tools for integrating all of our knowledge obtained from simulations,sensitivity analyses, and laboratory and field experimentation. This integration will often lead to more
coherent and rigorous conceptual models for the underlying coupled flow, transport, and reactions processes.

REFERENCES
Ahuja L.R. and Hebson C. (1992) Root Zone Water Quality Model ,GPSR Technical Report No. 2, USDA, ARS, Fort Collins.
Bradford S.A., Yates S.R., Bettehar M. and Šimůnek J. (2002)
Physical factors affecting the transport and fate of colloids
in saturated porous media. Water Resources Research, 38(12),1327, doi:10.1029/2002WR001340, 63.1 – 63.12.
Casey F.X.M., Larsen G.L., Hakk H. and Šimůnek J. (2003)
Fate and transport of 17β-Estradiol in soil-water systems.
Environmental Science and Technology, 37(11), 2400 – 2409.
Casey F.X.M. and Šimůnek J. (2001) Inverse analyses of thetransport of chlorinated hydrocarbons subject to sequential
transformation reactions. Journal of Environmental Quality,
30(4), 1354 – 1360.
Celia M.A., Bououtas E.T. and Zarba R.L. (1990) A generalmass-conservative numerical solution for the unsaturated flow
equation. Water Resources Research, 26, 1483 – 1496.
Fayer M.J. (2000) UNSAT-H Version 3.0: Unsaturated Soil Waterand Heat Flow Model . Theory, User Manual, and Examples.
Pacific Northwest National Laboratory 13249.
Flerchinger G.N., Hanson C.L. and Wight J.R. (1996) Modelingevapotranspiration and surface energy budgets across awatershed. Water Resources Research, 32, 2539 – 2548.
Gardner W.R. (1958) Some steady-state solutions of theunsaturated moisture flow equation with application toevaporation from a water table. Soil Science, 85, 228 – 232.
Gerke H.H. and van Genuchten M.T.h (1993) A dual-porosity
model for simulating the preferential movement of water andsolutes in structured porous media. Water Resources Research,
29, 305 – 319.
Hansen S., Jensen H.E., Nielsen N.E. and Svendsen, H. (1990)
DAISY: Soil Plant Atmosphere System Model , NPO Report
No. A 10, The National Agency for Environmental Protection,
Copenhagen, pp. 272.
Healy R.W. (1990) Simulation of Solute Transport in Variably Saturated Porous Media with Supplemental Information on Modifications to the U.S. Geological Survey’s Computer Program
VS2DI , Water-Resources Investigation Report 90 – 4025, U.S.
Geological Survey, p. 125.

1179

Hills R.G., Hudson D.B., Porro I. and Wierenga P.J. (1989)
Modeling one-dimensional infiltration into very dry soils, 1.
Model development and evaluation. Water Resources Research,
25(6), 1259 – 1269.
Hutson J.L. and Wagenet R.J. (1992) LEACHM:Leaching
Estimation and Chemistry Model, Research Series no. 92-3,
Cornell University: Ithaca.
Huyakorn P.S., Springer E.P., Guvanasen V. and Wadsworth T.D.
(1986) A three-dimensional finite-element model for simulatingwater flow in variably saturated porous media. Water Resources
Research, 22, 1790 – 1808.
Jacques D., Šimůnek J., Mallants D. and van Genuchten
M.T.h (2002) Multicomponent transport model for variablysaturated porous media: application to the transport ofheavy metals in soils. In Computational Methods in Water
Resources, Developments in Water Science 47 , Hassanizadeh
S.M., Schotting R.J., Gray W.G. and Pinder G.F. (Eds.), XIVthInternational Conference on Computational Methods in Water
Resources, June 23 – 28 2002, Elsevier: Delft, pp. 555 – 562.
Jansson P.-E. and Karlberg L. (2001) Coupled Heat and Mass
Transfer Model for Soil-Plant-Atmosphere Systems, RoyalInstitute of Technology, Department of Civil and Environmental
Engineering: Stockholm, p. 325.
Jarvis N.J. (1994) The MACRO Model (Version 3.1), Technical
Description and Sample Simulations, Reports and Dissertations
19, Department Soil Science Swedish University of Agricultural
Sciences, Uppsala, p. 51.
Jury W.A., Spencer W.F. and Farmer W.J. (1983) Behavior
assessment model for trace organics in soil: I. Modeldescription. Journal of Environmental Quality, 12(4), 558 – 564.
Kirkland M.R., Hills R.G. and Wierenga P.J. (1992) Algorithms
for solving Richards’ equation for variably saturated soils.
Water Resources Research, 28(8), 2049 – 2058.
Leij F.J. and Bradford S.A. (1994) 3DADE: A Computer Programfor Evaluating Three-Dimensional Equilibrium Solute Transport
in Porous Media, Research Report No. 134, U. S. Salinity
Laboratory, USDA, ARS, Riverside.
Leij F.J., Priesack E. and Schaap M.G. (2000) Solute transportmodeled with green’s functions with application to persistentsolute sources. Journal of Contaminant Hydrology, 41,
155 – 173.
Leij F.J. and Toride N. (1997) N3DADE: A Computer Programfor Evaluating Nonequilibrium Three-Dimensional Equilibrium
Solute Transport in Porous Media, Research Report No. 143,
U. S. Salinity Laboratory, USDA, ARS, Riverside.
Mendoza C.A., Therrien R. and Sudicky E.A. (1991)
ORTHOFEM User‘s Guide, Version 1.02. Waterloo Centre
for Groundwater Research, University of Waterloo: Waterloo,
Ontario.
Milly P.C.D. (1985) A mass-conservative procedure for timestepping in models of unsaturated flow. Advances in Water
Resources, 8, 32 – 36.
Neuman S.P. (1973) Saturated-unsaturated seepage by finite
elements. Proceedings of the ASCE, Journal of Hydraulic
Division, 99, 2233 – 2250.
Paniconi C. and Putti M. (1994) A comparison of Picard andNewton iteration in the numerical solution of multidimensionalvariably saturated flow problems. Water Resources Research,30(12), 3357 – 3374.

1180

SOILS

Parker J.C. and van Genuchten M.T.h (1984) Determining
Transport Parameters from Laboratory and Field Tracer
Experiments, Bulletin 84-3, Va Agricultural Experiment
Station: Blaacksburg.
Perrochet P. and Berod D. (1993) Stability of the standard CrankNicholson-Galerkin scheme applied to the diffusion-convection
equation: some new insights. Water Resources Research, 29(9),3291 – 3297.
Philip J.R. (1969) Theory of infiltration. Advances inHydroscience, 5, 215 – 296.
Philip J.R. (1992) What happens near a quasi-linear point source?
Water Resources Research, 28, 47 – 52.
Pruess K. (1991) Tough2 – A General-Purpose Numerical
Simulator for Multiphase Fluid and Heat Flow , Report LBL29400, Lawrence Berkeley Laboratory, Berkeley.
Ross P.J. (1990) Efficient numerical methods for infiltration
using Richards’ equation. Water Resources Research, 26(2),279 – 290.
Schaerlaekens J., Mallants D., Šimůnek J., van Genuchten M.T.hand Feyen J. (1999) Numerical simulation of transport andsequential biodegradation of chlorinated aliphatic hydrocarbons
using CHAIN− 2D. Journal of Hydrological Processes, 13(17),2847 – 2859.
Schijven J. and Šimůnek J. (2002) Kinetic modeling of virustransport at field scale. Journal of Contaminant Hydrology,55(1 – 2), 113 – 135.
Seuntjens P., Tirez K., Šimůnek J., van Genuchten M.T.h,
Cornelis C. and Geuzens P. (2001) Aging effects on cadmiumtransport in undisturbed contaminated sandy soil columns.
Journal of Environmental Quality, 30, 1040 – 1050.
Šimůnek J., Jarvis N.J., van Genuchten M.T.h and Gärdenäs A.
(2003) Nonequilibrium and preferential flow and transport in
the vadose zone: review and case study. Journal of Hydrology,
272, 14 – 35.
Šimůnek J. and Suarez D.L. (1994) Two-dimensional transportmodel for variably saturated porous media with major ion
chemistry. Water Resources Research, 30(4), 1115 – 1133.
Šimůnek J. and Suarez D.L. (1997) Sodic soil reclamation using
multicomponent transport modeling. Journal of Irrigation andDrainage Engineering-ASCE, 123(5), 367 – 376.
Šimůnek J., Šejna M. and van Genuchten M.T.h (1998)
The HYDRUS-1D Software Package for Simulating the OneDimensional Movement of Water, Heat, and Multiple Solutes
in Variably-Saturated Media, Version 2.0, IGWMC – TPS-70,
International Ground Water Modeling Center, Colorado Schoolof Mines: Golden, p. 202.
Šimůnek J., Šejna M. and van Genuchten M.T.h (1999a)
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 Schoolof Mines: Golden, p. 251.
Šimůnek J., van Genuchten M.T.h, Šejna M., Toride N.
and Leij F.J. (1999b) 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, p. 32.
Šimůnek J. and Valocchi A.J. (2002) Geochemical transport.
In Methods of Soil Analysis, Part 1, Physical Methods, Third
edition, Chap. 6.9, Dane J.H. and Topp G.C. (Eds.), Soil
Science Society of America, Madison, pp. 1511 – 1536.
Steefel C.I. and Van Cappellen P. (1998) Reactive transportmodeling of natural systems. Journal of Hydrology, 209, 1 – 7.
Sudicky E.A. (1989) The Laplace transform Galerkin technique:a time-continuous finite element theory and application to masstransport in groundwater. Water Resources Research, 25(8),1833 – 1846.
Toride N., Leij F.J. and van Genuchten M.T.h (1995) TheCXTFIT Code for Estimating Transport Parameters fromLaboratory or Field Tracer Experiments, Version 2.0, Research
Report No. 137, U. S. Salinity Laboratory, USDA, ARS,
Riverside.
van Dam J.C., Huygen J., Wesseling J.G., Feddes R.A., Kabat P.,
van Valsum P.E.V., Groenendijk P. and van Diepen C.A.
(1997) Theory of SWAP, Version 2.0. Simulation of Water
Flow, Solute Transport and Plant Growth in the Soil- WaterAtmosphere- Plant Environment, Department Water Resources,
WAU, Report 71, DLO Winand Staring Centre, TechnicalDocument 45, Wageningen.
van Genuchten M.T.h (1980) Determining Transport Parametersfrom Solute Displacement Experiments, Research Report No.
118, U. S. Salinity Laboratory, USDA, ARS, Riverside.
van Genuchten M.T.h (1981) Non-Equilibrium TransportParameters from Miscible Displacement Experiments, Research
Report No. 119, U. S. Salinity Laboratory, USDA, ARS,
Riverside.
van Genuchten M.T.h (1985) Convective-dispersive transportof solutes involved in sequential first-order decay reactions.
Computers & Geosciences, 11(2), 129 – 147.
van Genuchten M.T.h and Šimůnek J. (1996) Evaluation ofpollutant transport in the unsaturated zone. In Regional
Approaches to Water Pollution in the Environment, NATO ASI
Series: 2. Environment, Rijtema P.E. and Eliáš V. (Eds.),
Kluwer: Dordrecht, 139 – 172.
Verburg K., Ross P.J. and Bristow K.L. (1996) SWIMv2.1 User
Manual , Divisional Report 130, CSIRO.
Warrick A.W. and Knight J.H. (2002) Two-dimensional
unsaturated flow through a circular inclusion. Water Resources
Research, 38(7), 18.1 – 18.6, 10.1029/2001WR001041.
Williams G.A., Miller C.T. and Kelley C.T. (2000) Transformation approaches for simulating flow in variably saturated porous
media. Water Resources Research, 36(4), 923 – 934.
Wooding R.A. (1968) Steady infiltration from large shallow
circular pond. Water Resources Research, 4, 1259 – 1273.
Yeh G.T. (1981) Analytical Transient One-, Two-, and ThreeDimensional Simulation of Waste Transport in the Aquifer
System, ORNL-5602, Oak Ridge National Laboratory: Oak
Ridge.