STANMOD: MODEL USE, CALIBRATION,
AND VALIDATION
M. Th. van Genuchten, J. Šimůnek, F. J. Leij, N. Toride, M. Šejna

ABSTRACT. This article provides an overview of STANMOD, a Windows-based computer software package for evaluating
solute transport in soils and groundwater using analytical solutions of the advection-dispersion equation. The software integrates seven separate codes that have been popularly used over the years for a broad range of one-dimensional andmulti-dimensional solute transport applications: the CFITM, CFITIM, CXTFIT, CHAIN, and SCREEN models for onedimensional transport, and the 3DADE and N3DADE models for multi-dimensional transport. All of the models can be
run for direct (forward) problems, and several (CFITM, CFITM, CXTFIT, and 3DADE) can also be run for inverse problems. CXTFIT further includes a stochastic stream tube model assuming local-scale equilibrium or nonequilibriumtransport conditions. The 3DADE and N3DADE models apply to two- and three-dimensional transport during steady unidirectional water flow assuming equilibrium and nonequilibrium transport, respectively. Nonequilibrium transport can be
simulated using the assumption of either physical nonequilibrium (two-region or mobile-immobile type transport) or
chemical nonequilibrium (two-site partial equilibrium, partial kinetic sorption). The STANMOD software comes with alarge number of example applications illustrating the utility of the different codes for a variety of laboratory and fieldscale solute transport problems.
Keywords. Analytical solutions, Contaminant transport, Nonequilibrium transport, Parameter estimation, STANMOD.

T

he solution of many subsurface contaminanttransport problems requires the use of appropriate
modeling tools consistent with the application.
While many field problems require comprehensive numerical models simulating transient water, heat,
and/or solute transport (e.g., Šimůnek et al., 2012), somesolute transport problems may well be addressed using
simplified one-dimensional or multi-dimensional analyticalmodels. As pointed out by Javandel et al. (1984), Leij et al.
(1991), Bauer et al. (2001), Vanderborght et al. (2005), Perez-Guerrero et al. (2009), and many others, analyticalmodels are useful for a variety of applications, such as forproviding initial or approximate analyses of alternative pollution scenarios, conducting sensitivity analyses to investigate the effects of various parameters or processes on contaminant transport, extrapolating results over large times
and spatial scales where numerical solutions become impractical, serving as screening models, estimating transport

Submitted for review in September 2011 as manuscript number SW
9428; approved for publication by the Soil & Water Division of ASABE in
May 2012.
The authors are Martinus Th. van Genuchten, Professor, Departmentof Mechanical Engineering, Federal University of Rio de Janeiro, Brazil;Jiří Šimůnek, Professor, Department of Environment Sciences, University of California, Riverside, California; Feike J. Leij, Professor, Department of Civil Engineering and Construction Engineering Management,
California State University, Long Beach, California; Nobuo Toride, Professor, Faculty of Bioresources, Mie University, Tsu, Japan; and Miroslav
Šejna, Director, PC-Progress, Ltd., Prague, Czech Republic. Corresponding author: M. Th. van Genuchten, Department of Mechanical Engineering COPPE/LTTC, Federal University of Rio de Janeiro, UFRJ, Rio deJaneiro, RJ, 21945-970, Brazil; phone: +55-21-3264-7025; e-mail: rvangenuchten@hotmail.com.

parameters from laboratory or well-defined field experiments, providing benchmark solutions for more complex
transport processes that cannot be solved analytically, andfor validating more comprehensive numerical solutions of
the governing transport equations. At the same time, it is also important to recognize the limitations of analyticaltransport modeling. For example, analytical models are restricted to linear or linearized problems involving homogeneous profiles subject to a uniform unidirectional flow fieldthat is constant in time and space (e.g., no layered profiles
or any root water uptake can be considered). In addition,model parameters generally must be constant, while only
relatively simple initial and boundary conditions are permitted.
To have flexibility in optimally addressing different
transport problems, one may need a toolbox containing a
range of computer programs of varying degrees of complexity and dimensionality. Here we describe the Windowsbased STANMOD computer software package, which
combines a large number of analytical solute transportmodels that were developed during the past 30 years or
more. Most of the models within STANMOD were developed jointly by the U.S. Salinity Laboratory and the University of California, both in Riverside, California, andplaced in the public domain (Šimůnek et al., 1999; webbased updates available at www.pc-progress.com/
en/Default.aspx?stanmod). The objective of this article is tobriefly describe the various codes integrated within
STANMOD, which stands for “STudio of ANalytical
MODels,” and to illustrate typical laboratory- and fieldscale applications.

Transactions of the ASABE
Vol. 55(4): 1353-1366

2012 American Society of Agricultural and Biological Engineers ISSN 2151-0032

1353

Model
CFITM

CFITIMCXTFIT

CHAIN
SCREEN3DADE

N3DADE

Table 1. Overview of computer codes included in the STANMOD software package
Transport
InitialBoundary
Parameter
DomainGeneral Description
Condition
Condition
Estimation?
1D
BTCs only,
Constant
Constant orYes
Finite or
Equilibrium transport,
finite pulse
semi-infinite
No production or decay
profile1D
BTCs only,
Constant
Constant orYes
Semi-infinite
Equilibrium or
finite pulse
profilenonequilibrium transport,
No production or decay
1D
Equilibrium orConstant, Dirac,
Constant, DiracYes
Semi-infinite
nonequilibrium transport,
exponential, or
exponential, or
profileStream tube models,or multiple-step
multiple-step
Zero-order production,functionfunctionsFirst-order degradation
1D
Solute decay chains,
ZeroConstant orNo
Semi-infinite
No production,
initial condition
finite pulse
profileFirst-order degradation
1D
Pesticide screening,
Single-step function
VolatilizationNo
Semi-infinite
First-order degradation,versus distanceboundary
profileVolatilization
condition
2,3D
Equilibrium transport,
Constant, finite
Constant orYes
Semi-infinite
Zero-order production,parallelepipedal
finite pulse
Cartesian orFirst-order degradationor cylindricalcylindricalinitial condition
coordinates
2,3D
Nonequilibrium transport,
Constant, finite
Constant orNo
Semi-infinite
Zero-order production,parallelepipedal,
finite pulse
Cartesian orFirst-order degradationor cylindricalcylindricalinitial condition
coordinates

STANMOD DESCRIPTION
Table 1 gives an overview of the different codes incorporated thus far in STANMOD, along with references forthe original codes. The models include the CFITM code forpredicting or analyzing measured solute breakthrough data(concentration curves versus time) in terms of the onedimensional equilibrium advection-dispersion equation(ADE), the CFITIM code for similar problems permitting
both physical or chemical nonequilibrium transport, theCXTFIT code allowing predictions and analyses of equilibrium and nonequilibrium transport in time and/or space (including probabilistic analyses of solute transport using
stream tube models), the CHAIN code for predicting thetransport of solutes subject to consecutive first-order decay
chain reactions, the SCREEN code for environmental assessments (transport, degradation, adsorption, and volatilization) of soil-applied organic chemicals, and the 3DADE
and N3DADE codes for analyses of multi-dimensionaltransport problems assuming equilibrium and nonequilibrium conditions, respectively.
All programs within STANMOD are for forward problems in that solute concentrations can be predicted for aprescribed set of transport parameters, such as the porewater velocity, the dispersion coefficient and, if applicable,
zero-order production and first-order degradation coefficients. Several codes (CFITM, CFITIM, CXTFIT, and
3DADE) also permit inverse analyses, in which selectedtransport parameters can be estimated from observed concentration distributions versus time and/or distance. Parameter estimation in STANMOD is accomplished using a
Marquardt-Levenberg type weighted nonlinear leastsquares optimization approach (Marquart, 1963) that minimizes the objective function O:1354

O ( b) =

n


i =1

Key
Reference
van Genuchten(1980)
van Genuchten(1981a)Toride et al.
(1999)

van Genuchten(1985)
Jury et al.
(1983)
Leij and Bradford(1994)

Leij and Toride(1997)

2wi  ci* ( x, t ) − ci ( x, t; b) 



(1)

where n is the number of concentration measurements;
ci*(x, t) are observed concentrations at time t and location x(in one, two, or three dimensions); ci(x, t; b) represent corresponding model predictions for the vector b of unknown
transport parameters; and wi are weights associated with aparticular concentration data point.
The Marquardt-Levenberg parameter estimation approach as implemented in the STANMOD codes assumes
that the variance-covariance (weighting) matrices, which
provide information about the measurement accuracy, arediagonal (Šimůnek and Hopmans, 2002). The method uses
a local optimization gradient procedure that requires initial
estimates of the parameters to be optimized. The behaviorof objective function O in the neighborhood of this initial
estimate is used to select a direction vector, from which updated values of the unknown parameter vector b in equation 1 are determined. Depending upon the problem being
considered (i.e., the magnitude of the measurement errors
and the number of parameters being optimized), the objective function may sometimes lack a well-defined global
minimum, or may have several local minima in parameter
space. The optimization may then become sensitive to theinitial values of the optimized parameters. Depending on
the initial estimate, the final solution of the calibration insuch cases may then not be the global minimum but insteada local minimum. Consequently, we generally recommend
repeating the minimization problem with different initial
estimates of the optimized parameters, and then selecting
those parameter values among the different runs that pro-TRANSACTIONS OF THE ASABE

vide the lowest value of the objective function O(b). Generally, however, the parameter estimation approach in
STANMOD has proven to be very robust for most transport
problems, unless a large number of parameters are determined simultaneously (e.g., for nonequilibrium transport)
from data that are noisy or that are not providing good resolution of the expected concentration distribution, or when
the inverse solution shows strong correlation between someof the estimated parameters (such as sometimes is the casebetween the pore-water velocity and the retardation factor).
Not further discussed in this article is the general issueof model validation as applied to analytical solutions. Manyreferences are provided showing a wide range of applications of the various codes within STANMOD, including
comparisons with observed data versus distance and/or
time. Still, as mentioned earlier, the codes in most cases are
restricted to very specific situations, such as homogeneous
soils and constant fluid flow regimes typical of analyticalsolutions. Validation of simplified or approximate models isinherently problematic in that results cannot be extrapolated to more general transient flow conditions generally encountered in the field. We also refer to discussions by
Konikow and Bredehoeft (1992) and Oreskes et al. (1994),
among many others later, about the presumed difficulties ofmodel validation in general. Still, we note that the mathematical accuracy of all solutions and codes within
STANMOD has been tested thoroughly (“mathematicallyverified”) through comparisons with other codes, against
numerical software such as the HYDRUS codes (Šimůnek
et al., 2012), or otherwise as documented in the various
manuals.
We next briefly describe the different analytical solutetransport models included in STANMOD and discuss sometypical examples, most of which are installed with thesoftware. The software itself can be downloaded freely
from the STANMOD homepage (www.pc-progress.com/
en/Default.aspx?stanmod). Detailed descriptions of most or
all models are given in the original manuals, which can be
downloaded from several websites, including the
STANMOD site. The manuals provide complete descriptions of the governing transport equations, applicable initialand boundary conditions, the derived analytical solutions,listings of the older FORTRAN programs, as well as detailed descriptions of illustrative examples. The graphicsbased user-interface of STANMOD is for MS Windows environments and is mostly based on libraries developed for
the HYDRUS-1D and HYDRUS (2D/3D) software packages (Šimůnek et al., 2008a, 2008b; Šejna et al., 2011). All
computational programs within STANMOD were written in
FORTRAN (essentially identical to the original codes) andthe graphical user interfaces in MS Visual C++.

THE CFITM CODE
The CFITM code (van Genuchten, 1980) may be used toanalyze observed solute concentration distributions versus
time using analytical solutions of the one-dimensional equilibrium advection-dispersion equation given by:

55(4): 1353-1366

R

∂c
∂ 2c
∂c
=D
−v
2
∂t
∂x
∂x

(2)

where c is the solution concentration, x is distance, t is
time, D is the dispersion coefficient, v is the average porewater velocity (water flux q divided by the water content
θ), and R is the retardation factor, defined as:
R = 1+ρkθ

(3)where ρ is the dry soil bulk density, and k (sometimes denoted as Kd) is a linear partitioning or distribution coefficient of the solute between the liquid and solid phases.
The CFITM code assumes the use of dimensionless parameters, in which the initial and input (boundary) concentrations are scaled to zero and one, respectively. Initial applications of the code focused on the analyses of theobserved effluent curves from finite laboratory columns oflength L. Introduction of the column Peclet number, P, given by:
vL
D

P=

(4)

leads then to the dimensionless ADE as follows:R∂c 1 ∂ 2 c ∂c=−
∂T P ∂X 2 ∂X

(5)where dimensionless time T (often referred to as pore volume) and distance X are given by:
T=

vL
L

(6a)

x
L

(6b)

X =

and c is the normalized concentration, in contrast withequation 2.
CFITM provides analytical solutions of equation 5 subject to both first-type (Dirichlet or concentration type) andthird-type (Cauchy or flux type) boundary conditions at the
inlet (x = 0), and either a zero concentration gradient at theoutlet of the finite column (0 ≤ x ≤ L) or assuming that thefinite column is part of a semi-infinite system (0 ≤ x < ∞).
The concentration at the outlet is then assumed not to be affected by the exit boundary. As shown elsewhere (e.g., vanGenuchten and Parker, 1984; Skaggs and Leij, 2002), theproper mass-conservative analytical solution to use is theclassical solution of Lapidus and Amundson (1952) for theflux-averaged concentration leaving the column at x = L:c ( X ,T ) ≡ A( X ,T ) = P ( RX − T ) 1
erfc 

2
 2 RT P ( RX − T ) 1+ e PX erfc 
2
 2 RT


(7)

1355

THE CFITIM CODE

Figure 1. Calculated breakthrough curves for different values of theretardation factor (R = 1, 2, 3, 4, and 5) assuming P = 50.

The CFITIM code (van Genuchten, 1981a) is similar to
CFITM but extended to analytical solutions for physical orchemical nonequilibrium transport problems. As opposed
to CFITM, however, CFITIM only considers solutions forsemi-infinite media. Physical nonequilibrium refers to situations in which physical phenomena such as the presenceof immobile water pockets in the medium (e.g., as in aggregated soils or fractured rock) are responsible for thenonequilibrium situation. The governing equations in thiscase are (van Genuchten and Wierenga, 1976):(θ m + fρk )∂cm∂ 2 cm∂c− θ m vm m − α(cm − cim )= θ m Dm
2
∂t
∂x
∂x

[θim + (1 − f )ρk ]Figure 2. Observed (circles) and fitted (solid line) breakthroughcurves for Cr6+ transport through a 5 cm long soil column (P = 18.6, R= 1.348). Unpublished data from Wierenga (1980).

where erfc is the complementary error function.
Here we give two example applications using CFITM.
Figure 1 is for a relatively simple forward problem inwhich solute breakthrough curves (BTCs) are calculatedusing equation 7 for different values of the retardation factor (R) and assuming a value of 50 for the Peclet number(P).
Figure 2 is for an inverse problem in which equation 7 isfitted to observed Cr6+ data (Wierenga, 1980, unpublisheddata; van Genuchten, 1980). This is actually one of the
CFITM
examples(EXAMPLE2)
installed
with
STANMOD. The analytical solution in this case was for a
finite medium assuming a third-type input boundary condition and a zero concentration gradient at the exit boundary(solution not further given here), in which case the resident
and flux concentrations at the outlet become identical (vanGenuchten and Parker, 1984). The example illustrates theconvenient way in which transport parameters can be determined from concentration distributions by fitting analytical solutions to the observed data. The unknown parameters in this case were the column Peclet number (P) and theretardation factor (R) but could also include as needed thedimensionless pulse time (To = vto/L) when a finite solute
pulse of duration to is applied, rather than a single-step input function. The solution for the flux concentration leaving the column is then given by:
 A( X , T )c( X , T ) =  A( X , T ) − A( X , T − To )0 < t ≤ tot > to(8)(9b)

where the subscripts m and im refer to the mobile and immobile regions of the soil, f is the fraction of sorption sites
located in the mobile region, α is a first-order mass transfer
coefficient, and vm is the pore-water velocity for the mobile
phase (q/θm). Models based on equations 9a and 9b are alternatively referred to as mobile-immobile transport models, dual-porosity transport models, or two-region models.
Similar to equation 5, equations 9a and 9b can be put in
dimensionless form as follows:βR∂cm
1 ∂ 2 cm ∂cm−− ω(cm − cim )=∂T
Pm ∂X 2 ∂X
(1 − β) R∂cim= ω(cm − cim )
∂t

(10a)

(10b)

where R is the same as before, and additionally:
vm L
D

(11a)

θ + fρk
β= m
θR

(11b)

Pm =

ω=

αL
q

(11c)

The analytical solution of equations 10a and 10b, again
in terms of the flux-averaged concentration applicable tocolumn effluent curves, is given by equation 8 with A(x,t)
now defined as (Skaggs and Leij, 2002, among others):c( X , T ) =

T X

0 τ

βPm R
 P (βRX − τ) 
exp  − m

4πτ
4βRτ



× [1 − J (b, a) ] dτ

where A(x,t) is the solution for the constant inlet concentration (van Genuchten, 1980).

1356

∂cim
= α(cm − cim )
∂t

(9a)

(12)

where J(b,a) is Goldstein’s J function, and:
ωτ
a=
βR

(13a)

TRANSACTIONS OF THE ASABE

Figure 3. CFITIM window for selecting one of three possible
transport models, the desired inlet boundary condition (BC), the typeof problem being considered (direct or inverse), the maximum number of iterations for the iterative inverse solution, and the number ofdata points to be fitted.

b=ω (T − τ )(1 − β ) R(13b)Mathematically, very similar equations can also be formulated for chemical nonequilibrium based on the assumption that adsorption does not proceed at an equal rate for all
sorption sites on the solid phase. The one-site and two-site
nonequilibrium transport equations are examples of thisapproach. The dimensionless equations are then exactly thesame as given by equation 10 (van Genuchten and
Wierenga, 1976) but with a different definition of the dimensionless parameter β.
The dimensionless transport model given by equations
10a and 10b contains four dimensionless parameters thatcan be fitted to observed data: the column Peclet number(Pm), the retardation factor (R), a dimensionless nonequilibrium partitioning coefficient (β) having values between 0(all nonequilibrium) and 1 (all equilibrium), and a dimensionless mass transfer coefficient (ω). Additionally, again,
the dimensionless pulse time (To) could be fitted to the data
if the amount of solute mass entering the column must also
be estimated. Refer to the manual (van Genuchten, 1981a)
for a detailed discussion of the governing equations.
Here we show an example CFITIM application, based on
figure 7.9b of van Genuchten (1981a), but with a first-type
input boundary condition appropriate for flux-averaged concentrations typical of column effluent curves. In this case, thedimensionless nonequilibrium transport model was fitted toobserved boron effluent data from a 30 cm long soil column(van Genuchten, 1974). Figure 3 shows the CFITIM window
where the current project and its boundary conditions are defined, including the option of running an inverse problem(lower left corner of fig. 3), the maximum number of permitted iterations in the inverse problem (20 in this case), and thenumber of data points that are fitted (30 in this example).
Figure 4 shows the results of this application, including values of the fitted parameters.

THE CXTFIT CODE
To our knowledge, CFITM and CFITIM were the first

55(4): 1353-1366

Figure 4. Observed (circles) and fitted (solid line) breakthroughcurves for boron transport through a 30 cm long soil column filledwith Glendale clay loam (van Genuchten, 1974). The project involved
fitting five dimensionless parameters of the physical nonequilibrium
model to 30 data points (P = 37.9, R = 3.89, β = 0.625, ω = 0,650, andTo = 6.187).

computerized parameter estimation codes that consideredinverse analyses of laboratory soil column breakthroughdata in terms of equilibrium and nonequilibrium transportmodels. The models were updated later in very significant
ways, first by Parker and van Genuchten (1984) and subsequently by Toride et al. (1999), to produce the widely usedCXTFIT code. CXTFIT allows analyses of concentration
distributions versus time as well as versus depth (both forward and inverse problems), permits the use of the standard
equilibrium and nonequilibrium transport formulations, andincludes a variety of stochastic stream tube models thatconsider the effects of areal variations in the pore-water velocity on field-scale transport. The code can be applied to abroad range of laboratory and field data due to the implementation of very flexible initial and boundary conditions,the inclusion of general zero-order production and firstorder degradation processes, and by considering both equilibrium and nonequilibrium transport in the stochastic
stream tube models.
Examples included in the CXTFIT code are divided into
three separate groups (workspaces). The first group (Direct)
contains direct problems, the second group (Inverse) involves inverse problems, and the third group (Stochast)contains stochastic stream tube applications. Most examples were taken from the most recent CXTFIT manual (Toride et al., 1999), and book chapters by van Genuchten and
Cleary (1979) and Leij and Toride (1998). Here we give
one example for each of the three workspaces. Hundreds ofother examples can be retrieved readily from the peerreviewed literature, with applications involving thetransport of broad range of chemicals, including nonadsorbing tracers, pesticides, heavy metals, radionuclides, andother organic and inorganic chemicals. We list here rather
arbitrarily a few publications in which the CXTFIT codes
have been used successfully to estimate transport parameters from observed concentration distributions versus time
and/or distance: Gamerdinger et al. (1990), Jacobsen et al.
(1992a, 1992b), Pivetz and Steenhuis (1995), Spurlock etal. (1995), Huang et al. (1995), Jensen et al. (1996), Seuntjens et al. (2001), Kay and Conklin (2001), Kamra et al.
(2001), Antoniadis and McKinley (2003), Pace et al.
(2003), Gaur et al. (2003), Tilahun et al. (2005), Duke et al.
(2007), Wellman et al. (2008), Mayer et al. (2008), Srivastava et al. (2009), Tao et al. (2009), Köhne et al. (2011),

1357

and Koestel et al. (2011). Similar applications can be found
for the older CFITM and CFITIM codes but are not further
listed here.

DIRECT CXTFIT APPLICATIONThe CXTFIT code comes with many preprogrammed
forward problems, such as showing the effects of various
parameters, including of productions and decay, on distributions versus distance or time. Here we consider a direct
problem showing the effect of first-order degradation on
concentration distributions versus distance. The governingtransport equation for this case is given by:2R∂c∂ c ∂c=D- v − μc + γ∂t∂x 2 ∂x(14)

where μ is a first-order degradation rate constant, and γ is azero-order production term. Assuming a pulse type solute
application with concentration Co to a profile having an initial concentration Ci, the analytical solution for the volumeaveraged (resident) concentration is given by (van Genuchten and Alves, 1982; Toride et al., 1999):
Co B( x, t ) + Ci A( x, t )
0 < t ≤ to
 + γ [1 − A( x, t ) − B( x, t )] μ
c ( x, t ) = C B( x, t ) − B( x, t − t ) + C A( x, t )
o ]
i
 o[
 γ
t > to + (1 − A( x, t ) − B( x, t ) ) 
μ

(15)

where
 Rx − vt  μt   1
A( x, t ) = exp  −  1 − erfc 

2
R


 2 DRt − ( Rx − vt )2 
v 2t
exp  −πDR
4 DRt 1  vx v 2 t  Rx + vt  vx 
+ 1 + + exp   erfc 
2 
D DR 
D
 2 DRt  

B ( x, t ) =v
 Rx − ut  (v − u ) x 
exp 
erfc 

u+v
 2D 
 2 DRt −v
 Rx + ut  (v + u ) x exp  erfc  2 DRt u −v
 2D 
+ Rx + vt v2
 vx μt 
exp  −  erfc 

2μD
D R
 2 DRt 

u = v 2 + 4μD(16)

Co A( x, t ) + Ci [1 − A( x, t )]

0 < t ≤ to
 + γ B ( x, t ) Rc ( x, t ) = 
C A( x, t ) − A( x, t − t )o ] o[t > toγ+Ci [1 − A( x, t ) ] + B( x, t )
R


(19)

where
A( x, t ) =+1
 Rx − vt erfc 

2
 2 DRt  ( Rx − vt ) 2 
v 2t
exp  −πDR
4 DRt 

 Rx + vt 1  vx v 2t  vx 
− 1 + + exp   erfc D DR 
2 
D
 2 DRt 

(20)

DR 

 Rx − vt + v  Rx − vt B ( x, t ) = t +  erfc 2v
 2 DRt  ( Rx − vt )2 
2 DR t − Rx + vt + exp  −4πDR 
4 DRt v 
 t DR ( Rx + vt )2  Rx + vt  vx 
+ −+ exp   erfc 
2
4 DR 
D
 2 DRt  2 2v(21)

The above solutions are given in a structure that reflects thederivation of the analytical solutions using superposition principles (Toride et al., 1999) as well as entry of the model input
parameters into CXTFIT using its graphical interface. They
consist of the sum of solutions for a boundary value problem(the terms containing the input concentration Co in eqs. 15 and19), an initial value problem (the terms containing the initialcondition Ci), and a production value problem (the terms containing the zero-order production coefficient γ).
Figure 5 shows results for one example of a forward
problem with CXTFIT (FIG71 of the Direct workspace).
Calculated resident concentrations are presented versus distance 7.5 days after starting the application of a 5-day solute
pulse (to = 5 d) with a concentration Co = 1 mg L-3 to an initially solute-free medium (Ci = 0) and assuming degradation

(17)
(18)

Equation 15 holds for nonzero values of the first-orderdegradation coefficient (μ). When μ = 0, the solution reduces to (van Genuchten and Alves, 1982):
Figure 5. Calculated concentration distributions versus depth showing
the effect of degradation. The degradation coefficient increased fromμ = 0 (top curve) to μ = 1 d-1 (bottom curve).

1358

TRANSACTIONS OF THE ASABE

coefficients (μ) of 0, 0.25, 0.5, and 1 d-1. Other parameters
for this example are v = 25 cm d-1, D = 35.5 cm2 d-1, R = 3
and γ = 0.5 mg L-3 d-1. As expected, results indicate decreasing concentrations as the degradation coefficient increases.

INVERSE CXTFIT APPLICATIONIn this example, the pore-water velocity (v) and dispersion
coefficient (D) are estimated from observed concentrations atthree different depths (11, 17, and 23 cm) measured with
four-electrode conductivity sensors (Shiozawa, 1988, unpublished data). The experiments involved the continuous
application of a 0.01 M NaCl solution to an initially solutefree saturated sand (θ = 0.30), followed by leaching with solute free-water during unsaturated conditions (θ = 0.12). The
analytical solution for the resident concentration in this caseis given by equations 19 through 21 with γ = 0 (Lindstrom etal., 1967). Results for this example (projects FIG73A andFIG73B of the Inverse workspace of CXTFIT) are shown infigure 6, along with the fitted parameter values for v and D(R was assumed to be unity). Fitting all three curves simultaneously (either the breakthrough curves or the leaching
curves) with CXTFIT yielded about the same results as fitting v and D separately for each curve at the three depths.
However, the dispersivity (D/v) for the unsaturated sand wasfound to be about 3 times larger than the dispersivity for thesaturated experiment (Toride et al., 1999, 2003).
STOCHASTIC CXTFIT APPLICATIONCXTFIT also contains several examples demonstrating
the use of stochastic stream tube models (Toride et al.,1999). The modeling approach assumes that a field may beviewed as a collection of independent vertical soil columns,
often referred to as “stream tubes” (Dagan, 1993; Torideand Leij, 1996a; Vanderborght et al., 2006). Local-scale

Figure 6. Observed (open circles) and fitted (solid lines) breakthroughcurves for transport into a saturated sand (top) or leaching from unsaturated sand (bottom). In this example the three breakthroughcurves, measured at x = 11, 17, and 23 cm inside the column, were fitted simultaneously (yielding v = 2.50 cm min-1 and D = 0.130 cm2 min1), as were the three leaching curves (v = 0.253 cm min-1 and D =0.0391 cm2 min-1).

55(4): 1353-1366

transport in each stream tube can be described deterministically using one of several transport models. In CXTFIT,
these are the traditional equilibrium and nonequilibriumtransport models, as discussed earlier. Transport at the field
scale is then modeled by assuming that selected parametersin these transport models for each tube are the realization ofa stochastic process, with the mean solute concentration of
the entire field given by the ensemble average of the local
concentrations of all stream tubes. Stochastic variables in
the CXTFIT stream tube models are the pore-water velocity(v) in combination with either the dispersion coefficient(D), the partitioning coefficient for linear sorption (k), or the
first-order rate coefficient for nonequilibrium adsorption (α).
The three different pairs of random variables are describedusing bivariate lognormal joint probability density functions
denoted as f(v,D), f(v,k), and f(v,α) (Toride and Leij, 1996a,
1996b), Two additional stochastic parameters can be included, provided they are perfectly correlated with v.
We show here one forward and one inverse problem.
The forward problem demonstrates the effect of correlation(ρvk = -1, 0, +1) between the pore water velocity (v) and the
distribution coefficient (k) on calculated field-scale resident
concentration (cr) profiles (project FIG713 of the Stochastworkspace; figure 7.13 of Toride et al., 1999). Field-scaleconcentrations at t = 5 d resulting from a Dirac delta input
at t = 0 were calculated versus depth for perfect negative
correlation, no correlation, and perfect positive correlation
between v and k. Values of the other parameters in this example were <v> = 50 cm d-1, σv = 0.2, D = 20 cm2 d-1, <k>= 1 g-1 cm3, σk = 0.2, <R> = 5, and ρb/θ = 4 g cm-3, where <> indicates an ensemble average, σv is the standard deviation of the log transform of v, and σk is the standard deviation of the log transform of k. Results, shown in figure 7,
indicate that negative correlation between v and k (which
implies that v and R are inversely correlated) leads to increased spreading in the field-scale concentrations. Moredetails of this example are given by Toride and Leij (1996a)
and Toride et al. (1999).
The stochastic option of CXTFIT, together with parameter estimation, is demonstrated with the FIG712 project of
the Inverse workspace. The example (figure 7.12 of Toride
et al., 1999) concerns the analyses of observed resident
concentrations in a 0.64 ha field to which a bromide pulse
was applied for 1.69 d, followed by leaching with solutefree water (Jury et al., 1982). The stream tube model wasused to estimate the mean pore-water velocity <v>, theFigure 7. Effect of perfect positive and negative correlations and nocorrelation between v and k on field-scale resident concentration profiles. Negative correlation caused the most spreading.

1359

ry condition is of the form:
vfi (t )∂c= − D + vc ∂x x = 0 0

0 < t ≤ tot > to(23)

where
f1 (t ) = B1e − λ1tFigure 8. Observed (circles) and fitted values (solid line) of the resident concentration for field-scale bromide transport (after Toride etal., 1999).

mean dispersion coefficient <D>, and the standard deviation (σv) of the log transform of v from areally averaged
concentrations at a depth of 30 cm. The approach assumed
perfect positive correlation between v and D such that ρvD =+1 and σv = σD, which implies that the local-scale dispersivity, λ (= D/v), is the same for all stream tubes (Torideand Leij, 1996a). Results are shown in figure 8. Solution of
the inverse problem yielded <v> = 30.5 mm d-1, σv = σD =0.8, and <D> = 2.5 mm2 d-1. The example shows that localscale dispersion will have a relatively minor effect on fieldscale transport when the stream tube model is used.

THE CHAIN CODE
The CHAIN code (van Genuchten, 1985) may be used
to analyze the advective-dispersive transport of solutes subject to sequential first-order decay reactions. Examples are
the subsurface transport of various interacting nitrogen species, pesticides, chlorinated hydrocarbons, radionuclides,
pharmaceuticals, and explosives. The CHAIN code considers the transport of up to four species involved in thesetypes of decay chains. The governing transport equations
for such a system are given by:

R1

Ri

∂c1
∂ 2c
∂c
= D 21 − v 1 − μ1c1
∂t
∂x
∂x

∂ci ∂ 2 ci∂c=
− v i + μ i −1ci −1 − μ i ci
2
∂t
∂x
∂x
(i = 2, 4)(22a)f 3 (t ) = B4 e − λ1t + B5e − λ 2t + B6 e − λ 3t

(24)

f 4 (t ) = B7 e − λ1t + B8e − λ 2t + B9 e − λ 3t + B10 e − λ 4t

in which the coefficients λi (i = 1,4) and Bi (i = 1,10) are all
constants. The multiple terms of equation 24 are a consequence of decay reactions within the waste site itself (e.g.,
an industrial waste site, landfill, or nuclear waste repository), and the slow release of the solutes from the waste into
the soil profile (van Genuchten, 1985).
Many scenarios can be accounted for with the above
boundary conditions. For one particular release process, theconstants Bj are related to each other through the Bateman
equations (Bateman, 1910; Higashi and Pigford, 1980). Refer to the article by van Genuchten (1985) for a detailed
discussion of the decay chain boundary conditions, as wellas the resulting analytical solutions. Because of their complexity, these solutions are not repeated here. Here we provide only one relatively simple example involving thethree-species nitrification chain NH4+ → NO2- → NO3 (project NITROG in the CHAIN workspace of STANMOD).
Results are shown in figure 9 for the following parametervalues: v = 1 cm h-1, D = 0.18 cm2 h-1, t = 100 d, R1 = 2,R2 = R3 = 1, μ1 = 0.01 h-1, μ2 = 0.1 h-1, μ3 = 0, B1 = 1, Bi = 0(i = 2 through 10), λi = 0 (i = 1 through 4), and t = 200 h.
The above example is a classic application of consecutive decay chains, first published by Cho (1971). A few
other studies in which the CHAIN model within
STANMOD, or some of its predecessors, have been usedfor various solutes and applications are described by
Wernberg (1998), Pontedeiro et al. (2007), Srinivasan andClement (2008), Perez-Guerrero et al. (2009, 2010), andKasteel et al. (2010).

(22b)

where the subscript i refers to the ith member of the decay
chain. Like equation 14, equations 22a and 22b hold forsituations where degradation occurs only in the liquid
phase. When radionuclides are considered, decay occurs
equally in the liquid and adsorbed phases, which means that
retardation factors must be added to the first-order terms,
i.e., each μi in equations 22a and 22b must be replaced by
μiRi (van Genuchten, 1981b, 1985). This option is directly
available in CHAIN. For situations with unequal degradation rates in the liquid and solid phases, the value of each μi
must be adjusted and defined by the user upon input.
The analytical solutions in CHAIN hold for very general
first- and third-type inlet boundary conditions and a solutefree initial semi-infinite soil profile. For a third-type inlet
condition applicable to resident concentrations, the bounda-1360

f 2 (t ) = B2 e − λ1t + B3e − λ 2t

THE SCREEN CODEThe SCREEN code is based on a behavior assessment

Figure 9. Calculated concentration distributions versus depth of thenitrogen nitrification decay chain NH4+ → NO2- → NO3 (after Cho,
1971; van Genuchten, 1985).

TRANSACTIONS OF THE ASABE

model developed by Jury et al. (1983) to describe the fate
and transport of soil-applied organic chemicals. The model
assumes linear, equilibrium partitioning between the liquid,
solid, and gaseous phases, first-order degradation, andchemical losses to the atmosphere by volatilization througha stagnant air boundary layer above the soil surface. The
model is intended to classify and screen organic chemicals
for their relative susceptibility to different loss pathways(volatilization, leaching, and degradation) in soil and air.
SCREEN requires knowledge of the organic carbon partition
coefficient (Koc), Henry’s constant (Kh), and a net first-orderdegradation rate coefficient (μ) or the chemical half-life.
Jury et al. (1983) formulated their transport equations interms of the total concentration cT (= θRTc) as follows:∂cT∂ 2 cT∂c− ve T − μcT
= De
2
∂t
∂x
∂x

ρk aK H
+
θ
θ

(25)

(26)

in which a is the air content, KH is Henry’s law constant accounting for linear equilibrium partitioning between theliquid and air phases, and De and ve are the effective dispersion coefficient and pore-water velocity:
De =

D
RT

ve =

v
RT

(27a)

(27b)
Jury et al. (1983) solved the transport model assumingthat the chemical is uniformly incorporated at concentration
Ci within the upper part (0 ≤ x ≤ L) of a semi-infinite profile (zero concentration below x = L) and subject to a third-

Figure 10. Window showing various input parameters needed for thescreening model (SCREEN) of Jury et al. (1983).

55(4): 1353-1366∂cT= − H e CT
 − De ∂x + ve cT 
 x =0

(28)

where
He =

Da K HdθRT(29)

in which Da is the gaseous diffusion coefficient. Equation
28 accounts for diffusion from the soil surface through astagnant boundary layer of thickness d. The analytical solution of the above problem is (Jury et al., 1983):
Ciexp(−μt )2
 x − L − ve t  x − ve t   ×   erfc  − erfc  2 Dt 
 2 D t 
e
e 



 

cT ( x, t ) =where RT is the total retardation factor given by:
RT = 1 +type boundary condition at the soil surface that accounts forvolatilization into the air:


v x
v 
+ 1 + e  exp  e He De 
 x + L − ve t  x + ve t   × erfc  − erfc 
 





 2 De t  2 De t    ( H + v )( x + H e t )v +  2 + e  exp  e e− μt 
He 
De




 x + (2 H e + ve )t × erfc 
2 De t x + L + (2 H e + ve )t  
H L
− exp  e  erfc 
2 De t  De 

(30)

Here we give one example of the type of calculations
that are typically carried out with SCREEN. The example isa slight modification of the TEST3 project that comes withthe software. Various physical and solute transport parameters are entered in a window, as shown in figure 10, while
figure 11 shows chemical parameter values, which can beselected from a database or provided upon input. The pa-

Figure 11. SCREEN window showing chemical and other parameters
values used for the results shown in figure 11. The catalog of chemicals (upper right corner) can be used to selected chemical parametersfrom a database included in the software.

1361

with the initial and boundary conditions:c( x, y, z ) = Ci(− D

vC |y| ≤ a |z| ≤ b∂c+ v c)|x =0+ =  0∂x|y| > a |z| > b
0Figure 12. SCREEN window showing concentration distributions versus depth of several organic chemicals as obtained with the data
shown in figures 11 and 12.

rameter Koc in figure 11 is used to estimate the distribution
coefficient (k) from the fraction of organic matter (foc) usingk = focKoc.
Results of the transport model using the parameters fromfigures 10 and 11 are shown in figure 12. The example illustrates the interplay of sorption (more downward movement with lower Koc values), degradation (lower maximum
concentrations with increased degradation), and volatilization (lower total amount of mass in the profile and a moredispersed distribution, such as for PCE). Refer to articles
by Jury et al. (1983, 1984a, 1984b) for detailed descriptions
and several applications of their screening model.

THE 3DADE CODE
The 3DADE code (Leij and Bradford, 1994) may be
used to evaluate analytical solutions for three-dimensional
equilibrium solute transport in the subsurface. The solutions pertain to selected cases of three-dimensionaltransport during steady unidirectional water flow in porous
media having uniform flow and transport properties. The
transport equation contains terms to account for solute
movement by advection and dispersion, as well as for solute retardation, first-order decay, and zero-order production. Like CXTFIT, the 3DADE code can be used to solve
direct problems, in which concentrations are calculated as afunction of time and space for specified model parameters,
and for inverse problems, in which the program estimates
selected parameters by fitting one of the analytical solutions to specified experimental data. Analytical solutions
are provided in either Cartesian or cylindrical coordinate
systems.
Figure 13 shows available simulation options that can beconsidered, in this case a rectangular source at the soil surface (or in groundwater perpendicular to the lateral flow direction). The governing transport equation for this example
is given by:
∂c
∂ 2C
∂c
∂ 2C
∂ 2C
− v + Dy
+ Dz
− μc + γ
R = Dx
∂t
∂x
∂x 2
∂y 2
∂z 2

1362

(31)

(32)

(33)

∂c
( ∞, y , z , t ) = 0
∂x

(34)

∂c
( x, ±∞, z, t ) = 0
∂y

(35)

∂c
( x, y, ±∞, t ) = 0
∂z

(36)

where x, y, and z are the three spatial coordinates, with flow
occurring in the x direction; Dx is the longitudinal dispersion coefficient; Dy and Dz are transverse dispersion coefficients; and a and b define the lateral extent of the source at
the inflow boundary. The analytical solution for this example is as follows (Leij et al., 1991):
c( x, y , z, t ) =

μt
− 
γ 
1− e R 

μ



μτ
 ( Rx − vτ )2 
1 −R
−

exp
e
4 DRτ 4 πDx R  τ0 (a + y ) R  ( a − y ) R  
 + erf × erf  2 D y τ 
 2 D y τ  


t

vC0

+


 (b + z ) R  ( b − z ) R  × erf  + erf   dτ

 2 Dz τ 
 2 Dz τ  −t Rx + vτ  vx 
v 2C0
exp 

 erfc 
8 D x Rτ 2 Dx Rτ  Dx 0 (a + y ) R  ( a − y ) R  
 + erf ×  erf  2 D y τ 
 2 D y τ   (b + z ) R  ( b − z ) R  ×  erf  + erf   dτ
 2 Dz τ 
 2 Dz τ  

t

−γ
0

+

-μτ 

eR

 ( Rx − vτ )2 
v2τ


exp

3
 4 Dx Rτ  πDR Rx − vτ  1 
v 2 τ  (37)1vx+
erfc 
 1 +
−2R
 2 Dx Rτ  2 R  Dx Dx R  Rx + vτ  
 vx μτ × exp −  erfc   dt 2 Dx Rτ   Dx R 

TRANSACTIONS OF THE ASABE

Figure 13. Available geometries, initial conditions, and inlet boundary
conditions (BCs) that can be simulated using 3DADE within the
STANMOD software package.

Figure 14 shows graphical output (x-y planar distributions
for z = 0) obtained for this problem (EXAMPLE3 of
3DADE; figure 4 of Leij et al., 1991) assuming the following
values for the geometry and input parameters: v = 10 cm d-1,
R = 1, μ = γ = 0, Dx = 100 cm2 d-1, Dy = Dz = 10 cm2 d-1, Co =1 mg L-1, a = b = 7.5 cm, and t = 1 d. Several other examples, including for circular coordinate systems, are discussed
by Leij et al. (1991) and Leij and Bradford (1994).

THE N3DADE CODEThe N3DADE code of Leij and Toride (1997) is very
similar to 3DADE, except that analytical solutions areevaluated for three-dimensional nonequilibrium solutetransport. Similar to 3DADE, N3DADE evaluates analytical solutions for three-dimensional transport during steadyunidirectional water flow in macroscopically uniform porous media of semi-infinite length in the longitudinal direction, and of infinite length in the transverse directions.
However, unlike 3DADE, only direct (forward) problems
are considered in N3DADE.
As for the CFITIM and CXTFIT codes, nonequilibriumtransport can be simulated in terms of physical (mobileimmobile or dual porosity) formulations, and chemical
nonequilibrium assuming two sites (partial equilibrium,
partial kinetic sorption). The formulations as such are extensions of the CFITIM and CXTFIT one-dimensional
nonequilibrium transport equations to three dimensions.
Refer to the article by Toride et al. (1993) and theN3DADE manual by Leij and Toride (1997) for detailed
discussions of the governing transport equations, the analytical solutions, and several direct applications. Comprehensive sets of specific solutions are presented in thesepublications using Dirac, Heaviside, and exponential functions to describe a variety of initial, boundary, and production profiles. A rectangular or circular inflow area is specified for boundary value problems, while for initial andproduction value problems the respective initial and production profiles are defined for parallelepipedal, cylindrical, or spherical regions of the soil. Because of their complexity, these solutions are not restated here.
The N3DADE code within STANMOD comes with five
preprogrammed examples, all taken from Leij and Toride(1997). They are for calculated concentration distributions,versus time or spatially, for instantaneous solute application

Figure 14. Graphical display of results obtained with STANMOD for the 3DADE example depicted in figure 13.

55(4): 1353-1366

1363

from a disk at the soil surface (a problem with circular geometry), solute application from a rectangular source at thesoil surface (Cartesian geometry), and several initial andproduction value problems where the solute is initially located or produced within some finite circular or rectangularsource within the transport domain. We further note twoother applications of N3DADE to transport of a nonadsorbing dye (rhodamine) and cadmium in an alluvial aquifer by
Pang and Close (1999a, 1999b).

FUTURE DEVELOPMENTS
Notwithstanding their limitations in terms of being restricted to linear problems (and hence homogeneous profiles, unidirectional flow that is constant in time and space,
and the use of simplified initial and boundary conditions),
analytical solutions likely will remain popular for a large
number of applications. As shown also by many of the examples in this STANMOD overview, possible applications
include initial or approximate analyses of subsurface
transport, providing a better understanding of alternativetransport processes and the effects of key transport parameters, predicting transport for large spatial or time scales,and estimating model parameters using inverse methods.
The individual codes discussed in this review and included
in STANMOD were developed for exactly those types ofapplications. The primary contribution of STANMOD is toinclude all of these codes in an easy-to-use, Windowsbased framework for possible use by practitioners not necessarily familiar with the mathematical complexities of theseparate programs. The GUIs of STANMOD are for thispurpose very intuitive. They are used to manage the input
data required to run STANMOD, as well as for editing, parameter allocation, problem execution, and visualization of
results.
All computational programs were written in FORTRAN,
and the graphical interfaces were written in MS Visual
C++. The pre-processing unit includes specification of allnecessary parameters to successfully run the FORTRAN
codes. All post-processing is also carried out in the GUI.
The post-processing unit consists of simple x-y plots forgraphical presentation of the results (and data) and a dialog
window that displays an ASCII output file. The multidimensional 3DADE and N3DADE codes are further supported with output graphics that include 2D contours (isolines or color spectra, as shown in this article) in real or
cross-sectional views for equilibrium, nonequilibrium, andtotal concentrations. Output also includes animation ofgraphic displays at sequential time steps, and line graphs
for selected boundary or internal sections, as well as forplots of variables versus time. Areas of interest can bezoomed into, and the vertical scale can be enlarged forcross-sectional views.
STANMOD is in public domain and can be downloaded
freely (www.pc-progress.com/en/Default.aspx?stanmod).
The software is supported by a discussion forum where users can submit questions or suggestions about the models.
All projects included within STANMOD can also be downloaded separately form the website. While no definite plans

1364

exist at present, new codes geared toward analytical models
may be included in the future. In all, we believe that thedifferent codes within STANMOD are serving an importantrole in the transfer of software technologies to both the scientific community and practitioners.

REFERENCES
Antoniadis, V., and J. D. McKinley. 2003. Measuring heavy metal
migration rates in low-permeable soil. Environ. Chem. Lett.
1(1):103-106.
Bateman, H. 1910. The solution of a system of differential
equations occurring in the theory of radioactive
transformations. Proc. Cambridge Phil. Soc. 15: 423-427.
Bauer P., S. Attinger, and W. Kinzelbach. 2001. Transport of adecay chain in homogenous porous media: Analyticalsolutions. J. Contam. Hydrol. 49(3-4): 217-239.
Cho, C. M. 1971. Convective transport of ammonium withnitrification in soil. Canadian J. Soil Sci. 51(3): 339-350.
Dagan. G. 1993. The Bresler-Dagan model of flow and transport:
Recent theoretical developments. In Water Flow and Solute
Transport in Soils, 13-32. D. Russo and G. Dagan, eds. Berlin,
Germany: Springer-Verlag.
Duke, C. L., R. C. Roback, P. W. Reimus, R. S. Bowman, T. L.
McLing, K. E. Baker, and L. C. Hull. 2007. Elucidation offlow and transport processes in a variably saturated system ofinterlayered sediment and fractured rock using tracer tests.
Vadose Zone J. 6(4): 855-867.
Gamerdinger, A. P., R. J. Wagenet, and M. Th. van Genuchten.
1990. Application of two-site/two-region models for studying
simultaneous nonequilibrium transport and degradation ofpesticides. SSSA J. 54(4): 957-963.
Gaur, A., R. Horton, D. B. Jaynes, J. Lee, and S. A. Al-Jabri.
2003. Using surface time domain reflectometry measurements
to estimate subsurface chemical movement. Vadose Zone J.
2(4): 539- 543.
Higashi, K., and T. H. Pigford. 1980. Analytical models for
migration of radionuclides in geologic sorbing media. J. Nucl.
Sci. Tech. 17(9): 700-709.
Huang, K., N. Toride, and M. Th. van Genuchten. 1995.
Experimental investigation of solute transport in large,
homogeneous and heterogeneous, saturated soil columns.
Transport in Porous Media 18(3): 283-302.
Jacobsen, O. H., F. J. Leij, and M. Th. van Genuchten. 1992a.
Parameter determination for chloride and tritium transport inundisturbed lysimeters during steady flow. Nordic Hydrol.
23(2): 89-104.
Jacobsen, O., H., F. J. Leij, and M. Th. van Genuchten. 1992b.
Lysimeter study of anion transport during steady flow throughlayered coarse-textured soil profiles. Soil Sci. 154(3): 196-205.
Javandel, I., C. Doughty, and Chin-Fu Tsang. 1984. GroundwaterTransport: Handbook of Mathematical Models. Water
Resources Monograph No. 10. Washington, D.C.: American
Geophysical Union.
Jensen, K. H., G. Destouni, and M. Sassner. 1996. Advectiondispersion analysis of solute transport in undisturbed soilmonoliths. Ground Water 34(6): 1090-1097.
Jury, W. A., L. H. Stolzy, and P. Shouse. 1982. A field test of thetransfer function model for predicting solute transport. Water
Resour. Res. 18(2): 369-375.
Jury, W. A., W. F. Spencer, and W. J. Farmer. 1983. Behavior
assessment model for trace organics in soil: I. Description ofmodel. J. Environ. Qual. 12(4):558-564.
Jury, W. A., W. J. Farmer, and W. F. Spencer. 1984a. Behavior
assessment model for trace organics in soil: II. Chemical

TRANSACTIONS OF THE ASABE

classification and parameter sensitivity. J. Environ. Qual.
13(4): 4: 567-572.
Jury, W. A., W. F. Spencer, and W. J. Farmer. 1984b. Behavior
assessment model for trace organics in soil: IV. Review ofexperimental evidence. J. Environ. Qual. 13(4): 580-586.
Kamra, S. K., B. Lennartz, M. Th. van Genuchten, and P.
Widmoser. 2001. Evaluating non-equilibrium solute transport
in small columns. J. Contam. Hydrol. 48(3-4): 189-212.
Kasteel, R., T. Pütz, J. Vanderborght, and H. Vereecken. 2010.
Fate of two herbicides in zero-tension lysimeters and in fieldsoil. J. Environ. Qual. 39(4): 1451-1466.
Kay, J. T., and M. H. Conklin. 2001. Processes of nickel andcobalt uptake by a manganese oxide forming sediment in Pinal
Creek, Globe mining district, Arizona. Environ. Sci. Tech.
35(24): 4719-4725.
Koestel, J. K., J. Moeys, and N. J. Jarvis. 2011. Evaluation of
nonparametric shape measures for solute breakthrough curves.
Vadose Zone J. 10(4): 1261-1275.
Köhne, J. M., S. Schlüter, and H.-J. Vogel. 2011. Predicting solutetransport in structured soil using pore network models. Vadose
Zone J. 10(3): 1082-1096.
Konikow, L. F., and J. D. Bredehoeft. 1992. Ground-water models
cannot be validated. Advances in Water Resources 15(1): 7583.
Lapidus, L., and N. R. Amundson. 1952. Mathematics ofadsorption in beds: VI. The effect of longitudinal diffusion in
ion exchange chromatographic columns. J. Phys. Chem. 56(8):984-988.
Leij, F. J., and S. A. Bradford. 1994. 3DADE: A computer
program for evaluating three-dimensional equilibrium solutetransport in porous media. Research Report No. 134.
Riverside, Cal.: USDA-ARS U.S. Salinity Laboratory.
Leij, F. J., and N. Toride. 1997. N3DADE: A computer programfor evaluating nonequilibrium three-dimensional equilibrium
solute transport in porous media. Research Report No. 143.
Riverside, Cal.: USDA-ARS U.S. Salinity Laboratory.
Leij, F. J., and N. Toride. 1998. Chapter 5: Analytical solutions fornonequilibrium transport models. In Physical Nonequilibriumin Soils: Modeling and Application, 117-156. H. M. Selim and
L. Ma, eds. Chelsea, Mich.: Ann Arbor Press.
Leij, F. J., T. H. Skaggs, and M. Th. van Genuchten. 1991.
Analytical solutions for solute transport in three-dimensionalsemi-infinite porous media. Water Resour. Res. 27(10): 27192733.
Lindstrom, F. T., R. Haque, V. H. Freed, and L. Boersma. 1967.
The movement of some herbicides in soils: Linear diffusionand convection of chemicals in soils. Environ. Sci. Tech. 1(7):561-565.
Marquardt, D. W. 1963. An algorithm for least-squares estimationof nonlinear parameters. SIAM J. Appl. Math. 11(2): 431-441.
Mayer, A., T. Sandman, and M. Breidenbach. 2008. Effect of flow
regime on physical nonequilibrium transport in unsaturated
porous media. Vadose Zone J. 7(3): 981-991.
Oreskes, N., K. Shrader-Frechette, and K. Belitz. 1994.
Verification, validation, and confirmation of numerical models
in the earth sciences. Science 263(5147): 641-646.
Pace, M. N., M. A. Mayes, P. M. Jardine, T. L. Mehlhorn, J. M.
Zachara, and B. N. Bjornstad. 2003. Quantifying the effects ofsmall-scale heterogeneities on flow and transport inundisturbed cores from the Hanford formation. Vadose Zone J.
2(4): 664-676.
Pang, L., and M. Close. 1999a. Field-scale physicalnonequilibrium transport in an alluvia aquifer. J. Contam.
Hydrol. 38(4): 447-464.
Pang, L., and M. Close. 1999b. A field study of nonequilibrium
and facilitated transport of Cd in an alluvial gravel aquifer.

55(4): 1353-1366

Ground Water 37(5): 785-792.
Parker, J. C., and M. Th. van Genuchten. 1984. Determining
transport parameters from laboratory and field tracer
experiments. Bulletin 84-3. Blacksburg, Va.: Virginia
Agricultural Experiment Station.
Pérez Guerrero, J. S., L. C. G. Pimentel, T. H. Skaggs, and M. Th.
van Genuchten. 2009. Analytical solution of the advectiondiffusion transport equation using a change-of-variable andintegral transform technique. Intl. J. Heat and Mass Transfer
52(13-14): 3297-3304.
Pérez-Guerrero, J. S., T. H. Skaggs, and M. Th. van Genuchten.
2010. Analytical solution for multi-species contaminanttransport subject to sequential first-order decay reactions infinite media. Transport in Porous Media 80(2): 373-387.
Pivetz, B. E., and T. S. Steenhuis. 1995. Soil matrix andmacropore biodegradation of 2,4-D. J. Environ. Qual. 24(4):564-570.
Pontedeiro, E. M., P. F. L. Heilbron, and R. M. Cotta. 2007.
Assessment of the mineral industry NORM/TENORM
disposal in hazardous landfills. J. Hazard. Materials 139(3):563-568.
Šejna, M., J. Šimůnek, and M. Th. van Genuchten. 2011. The
HYDRUS software package for simulating two- and threedimensional movement of water, heat, and multiple solutes invariably saturated media. User Manual, Version 2.0. Prague,
Czech Republic: PC Progress.
Seuntjens, P., K. Tirez, J. Šimůnek, M. Th. van Genuchten, C.
Cornelis, and P. Geuzens. 2001. Aging effects on cadmiumtransport in undisturbed contaminated sandy soil columns. J.
Environ. Qual. 30(3): 1040-1050.
Šimůnek, J., and J. W. Hopmans. 2002. Chapter 1.7: Parameter
optimization and nonlinear fitting. In Methods of Soil Analysis:
Part 1. Physical Methods, 139-157. J. H. Dane and G. C. Topp,
eds. Madison, Wisc.: SSSA.
Šimůnek, J., M. Th. van Genuchten, M. Šejna, N. Toride, and F. J.
Leij. 1999. The STANMOD computer software for evaluating
solute transport in porous media using analytical solutions ofconvection-dispersion equation. Versions 1.0 and 2.0.
IGWMC-TPS-71. Golden, Colo.: Colorado School of Mines,
International Ground Water Modeling Center.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna. 2008a.
Development and applications of the HYDRUS andSTANMOD software packages and related codes. Vadose Zone
J. 7(2): 587-600.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M. Th. van
Genuchten. 2008b. The HYDRUS-1D software package forsimulating the movement of water, heat, and multiple solutes
in variably saturated media. Version 4.0. HYDRUS Software
Series 3. Riverside, Cal.: University of California, Departmentof Environmental Sciences.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna. 2012.
HYDRUS: Model use, calibration, and validation. Trans.
ASABE 55(4): 1263-1276.
Skaggs, T. H., and F. J. Leij. 2002. Chapter 6.3: Solute transport:
Theoretical background. In Methods of Soil Analysis: Part 4.
Physical Methods, 1353-1380. J. H. Dane and C.G. Topp, eds.
Madison, Wisc.: SSSA.
Spurlock, F. C., K. Huang, and M. Th. van Genuchten. 1995.
Isotherm nonlinearity and nonequilibrium sorption effects ontransport of fenuron and monuron in soil columns. Environ.
Sci. Tech. 29(4): 1000-1007.
Srinivasan, V., and T. P. Clement. 2008. Analytical solutions forsequentially coupled one-dimensional reactive transport
problems: Part II. Special cases, implementation, and testing.
Adv. Water Resour. 31(2): 219-232.
Srivastava, P., S. M. Sanders, J. H. Dane, Y. Feng, J. Basile, and1365

M. O. Barnett. 2009. Fate and transport of sulfadimethoxine
and ormetoprim in two southeastern United States soils.
Vadose Zone J. 8(1): 32-41.
Tao, Y. Q., Y. R. Bian, X. L. Yang, and F. Wang. 2009. Transportof malathion in homogeneous soil liquid chromatographic
columns: Influence of nonequilibrium sorption. Vadose Zone J.
8(1): 42-51.
Tilahun, K., J. F. Botha, and A. T. P. Bennie. 2005. Transport ofbromide in the Bainsvlei soil: Field experiment anddeterministic/stochastic model simulation. II. Intermittent
water application. Soil Res. 43(1): 81-85.
Toride, N., and F. J. Leij. 1996a. Convective dispersive stream
tube model for field-scale solute transport: I. Moment analysis.
SSSA J. 60(2): 342-352.
Toride, N., and F. J. Leij. 1996b. Convective dispersive stream
tube model for field-scale solute transport: II. Examples andcalibration. SSSA J. 60(2): 352-361.
Toride, N., F. J. Leij, and M. Th. van Genuchten. 1993. A
comprehensive set of analytical solutions for nonequilibrium
solute transport with first-order decay and zero-orderproduction. Water Resour. Res. 29(7): 2167-2182.
Toride, N., F. J. Leij, and M. Th. van Genuchten. 1999. TheCXTFIT code for estimating transport parameters fromlaboratory or field tracer experiments. Version 2.1. Research
Report No. 137. Riverside, Cal.: USDA-ARS U.S. Salinity
Laboratory.
Toride, N., M. Inoue, and F. J. Leij. 2003. Hydrodynamic
dispersion in an unsaturated dune sand. SSSA J. 67(3): 703712.
Vanderborght, J., R. Kasteel, M. Herbst, M. Javaux, D. Thiéry, M.
Vanclooster, C. Mouvet, and H. Vereecken. 2005. A set ofanalytical benchmarks to test numerical models of flow andtransport in soils. Vadose Zone J. 4(1): 206-221.
Vanderborght, J., R. Kasteel, and H. Vereecken. 2006. Stochastic
continuum transport equations for field-scale solute transport:
Overview of theoretical and experimental results. Vadose Zone
J. 5(1):184-203.
van Genuchten, M. Th. 1974. Mass transfer studies in sorbing

1366

porous media. PhD diss. Las Cruces, N.M.: New Mexico State
University.
van Genuchten, M. Th. 1980. Determining transport parametersfrom solute displacement experiments. Research Report No.
118. Riverside, Cal.: USDA-ARS U.S. Salinity Laboratory.
van Genuchten, M. Th. 1981a. Non-equilibrium transportparameters from miscible displacement experiments, Research
Report No. 119. Riverside, Cal.: USDA-ARS U.S. Salinity
Laboratory.
van Genuchten, M. Th. 1981b. Analytical solutions for chemicaltransport with simultaneous adsorption, zero-order production,and first-order decay. J. Hydrol. 49(3-4): 213-233.
van Genuchten, M. Th. 1985. Convective-dispersive transport ofsolutes involved in sequential first-order decay reactions.
Computers and Geosci. 11(2): 129-147.
van Genuchten, M. Th., and W. J. Alves. 1982. Analyticalsolutions of the one-dimensional convective-dispersive solutetransport equation. Tech. Bulletin 1661. Washington, D.C.:
USDA-ARS.
van Genuchten, M. Th., and R. W. Cleary. 1979. Movement ofsolutes in soils: Computer-simulated and laboratory results. In
Soil Chemistry: B. Physico-Chemical Models, 349-386. G. H.
Bolt, ed. Developments in Soil Science 5B. Amsterdam, The
Netherlands: Elsevier Scientific.
van Genuchten, M. Th., and J. C. Parker. 1984. Boundary
conditions for displacement experiments through shortlaboratory soil columns. SSSA J. 48(4): 703-708.
van Genuchten, M. Th., and P. J. Wierenga. 1976. Mass transfer
studies in sorbing porous media: I. Analytical solutions. SSSA
J. 40(4): 473-480.
Wellman, D. M., J. M. Zachara, C. Liu, N. P. Qafoku, S. C. Smith,and S. W. Forrester. 2008. Advective desorption of
uranium(IV) from contaminated Hanford vadose zone
sediments under saturated and unsaturated conditions. Vadose
Zone J. 7(4): 1144-1159.
Wernberg, T. 1998. Multicomponent groundwater transport withchemical equilibrium and kinetics: Model development andevaluation. Hydrol. Sci. J. 43(2): 299-317.

TRANSACTIONS OF THE ASABE