J. Hydrol. Hydromech., 61, 2013, 2, 146–160
DOI: 10.2478/johh-2013-0020

Exact analytical solutions for contaminant transport in rivers
1. The equilibrium advection-dispersion equation
Martinus Th. van Genuchten 1 *, Feike J. Leij 2, Todd H. Skaggs 3, Nobuo Toride 4,
Scott A. Bradford 3, Elizabeth M. Pontedeiro 5
1

Department of Mechanical Engineering, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, RJ, Brazil.
Department of Civil and Construction Engineering Management, California State University, Long Beach, California, USA.
3
U.S. Salinity Laboratory, USDA, ARS, Riverside, California, USA.
4
Faculty of Bioresources, Mie University, Tsu, Japan.
5
Department of Nuclear Engineering, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, RJ, Brazil.
*
Corresponding author. E-mail: rvangenuchten@hotmail.com
2

Abstract: Analytical solutions of the advection-dispersion equation and related models are indispensable for predicting
or analyzing contaminant transport processes in streams and rivers, as well as in other surface water bodies. Many useful
analytical solutions originated in disciplines other than surface-water hydrology, are scattered across the literature, andnot always well known. In this two-part series we provide a discussion of the advection-dispersion equation and related
models for predicting concentration distributions as a function of time and distance, and compile in one place a large
number of analytical solutions. In the current part 1 we present a series of one- and multi-dimensional solutions of thestandard equilibrium advection-dispersion equation with and without terms accounting for zero-order production andfirst-order decay. The solutions may prove useful for simplified analyses of contaminant transport in surface water, andfor mathematical verification of more comprehensive numerical transport models. Part 2 provides solutions for advective-dispersive transport with mass exchange into dead zones, diffusion in hyporheic zones, and consecutive decay chain
reactions.
Keywords: Contaminant transport; Analytical solutions; Surface water; Advection-dispersion equation.
INTRODUCTION
Mathematical models have long proved useful for analyzingor predicting the fate and transport of contaminants in streams
and rivers (Bencala, 1983; Fischer et al., 1979; Runkel, 1998;
Thomann, 1973), including contaminant exchange with fluvial
sediments and the surrounding stream bed (Thackston andSchnelle, 1970; Bencala and Walters, 1983; Wörman, 1998;
Anderson and Phanikumar, 2011), often referred to as the
hyporheic zone (Wörman, 1998; Runkel et al., 2003; Bencala,2005; Gerecht et al., 2011). Such models may be used also forsubsurface streams and karst systems (Field, 1997). Constituents in surface waters may involve a range of specific naturaland anthropogenic chemicals, including toxic trace elements,
radionuclides, industrial solvents, pesticides, nutrients, pathogenic microorganisms, pharmaceuticals, and a variety of water
quality variables such as total salinity or dissolved oxygen.
Because of the many complex and often nonlinear physical,
chemical and biological processes affecting contaminanttransport in streams and rivers, numerical models are now
increasingly used for prediction purposes (e.g., Anderson andPhanikumar, 2011; O’Connor et al., 2009; Runkel, 1998;
Runkel and Chapra, 1993). Still, analytical and quasi-analytical
approaches are useful for simplified analyses of a variety ofcontaminant transport scenarios, especially for relatively longspatial and time scales, when insufficient data are available towarrant the use of a comprehensive numerical model, and fortesting numerical models. The application of analytical approaches is greatly facilitated now also by mathematical software such as Maple, Mathematica and Matlab.
Although the literature contains many analytical solutions
that are applicable to transport in rivers and streams, most of
the solutions originated in other disciplines using mathematically equivalent equations, and hence their existence may not be

146familiar to many in the surface-water hydrology community(e.g. De Smedt et al., 2006; Huang, 2006). In this two-partseries we assembled a large number of one- and multidimensional solutions for contaminant transport in streams andrivers. The solutions are for one-dimensional advectivedispersive transport, longitudinal transport and lateral dispersion, transport with simultaneous first-order exchange withrelative immobile or stagnant water zones (transient storage
models), advective-dispersive transport with simultaneous
diffusion into and out of hyporheic zones, and transport ofsolutes subject to consecutive decay chain reactions. Most of
the solutions were derived from solutions to mathematicallyvery similar problems in subsurface contaminant transport, on
which the authors have worked for many years (e.g., Leij et al.,1991, Toride et al., 1993a, 1999; van Genuchten, 1981). Except
for solutions pertaining to diffusion in hyporheic zones as presented in part 2 (van Genuchten et al., 2013), most or all models have been incorporated in the public-domain windowsbased STANMOD software package (Šimůnek et al., 2000). In
this part 1 we first start with a brief introduction of the advection-dispersion equation governing transport in surface water
bodies, applicable initial and boundary conditions, and a brief
overview of analytical solution techniques.
GOVERNING TRANSPORT EQUATION
Formulation of the Advection-Dispersion equation
The transport of contaminants in surface and subsurface waters is generally described with the Advection-Dispersion Equation (ADE), or appropriate modifications or extensions thereof.
The ADE distinguishes two transport modes: advective
transport as a result of passive movement along with water, anddispersive/diffusive transport to account for diffusion andsmall-scale variations in the flow velocity as well as any otherExact analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equation

processes that contribute to solute spreading. Solute spreadingis generally considered to be a Fickian or Gaussian diffusion/dispersion process. A considerable body of literature on
dispersion processes now exists for both surface and subsurface
contaminant transport (e.g., Bear and Verruijt, 1987; Fischer etal., 1979; Jury and Horton, 2004).
For one-dimensional transport, the solute flux Js (ML-2T-1)
can be written as

J s = uC − Dx

∂C
,
∂x

(1)

where u is the longitudinal fluid flow velocity (LT-1), C is thesolute concentration expressed as mass per unit volume ofwater (ML-3), Dx is the longitudinal dispersion coefficient accounting for the combined effects of ionic or molecular diffusion and hydrodynamic dispersion (L-2T-1), and x is the longitudinal coordinate (L).
The mass balance equation can be formulated in a generalmanner by considering the accumulation or depletion of solutein a control volume over time as a result of the divergence ofthe flux (i.e., net inflow or outflow), possible reactions, and theinjection or extraction of solute along with the fluid phase:∂C= − ∇ ⋅ J s − Rs +RwCe ,∂t(2)

where t is time (T), and Rs represents arbitrary sinks (< 0) or
sources (> 0) of solute (ML-3T-1), while the last term denotes
injection (> 0) or pumping (< 0) of water with constituent concentration Ce at a rate Rw (L3L-3T-1).
The use of analytical methods requires that some of the underlying processes affecting transport are simplified or approximated. Such simplifications include the use of constant values
for u and Dx with respect to time and position, adopting idealized initial and boundary conditions, assuming that production
and degradation are limited to zero- or first-order processes,
and describing solute sorption by fluvial sediments using alinear equilibrium isotherm (or neglecting sorption altogether).
If the last two terms of the mass balance equation are ignored,
substitution of the advective-dispersive flux, Js, yields the
popularly used advection-dispersion equation (ADE) for onedimensional transport as follows

∂C
∂2 C
∂C
= Dx 2 − u
.
∂t
∂x
∂x

(3)

∂S
,
∂t

(5)

in which the distribution coefficient Kd (L3M-1) may be viewed
as the ratio of the sorbed (sediment) concentration and thedissolved (stream) concentration at equilibrium (Bencala et al.,1983). Substituting (1), (4) and (5) into (2) allows the ADE tobe formulated in terms of one dependent variable (i.e., thestream solute concentration) according to

R

∂C
∂2 C
∂C
= Dx
−u
2
∂t
∂x
∂x

(6)

in which the retardation factor R is given by

R = 1+ ρs K d .

(7)

Sorption onto the sediments reduces the apparent advective
and dispersive fluxes by a factor equal to R. This sorption isoften considered negligible, in which case Kd = 0 and the value
for R becomes equal to one. Analytical expressions for theconcentration can generally be obtained only for linear sorption. Because the mathematical problem is trivial with regard to
the value for R (the values of time t in the analytical solutions
are to be divided by R), sorption will not be considered explicitly in the remainder of our study.
Many other processes such as biodegradation or inactivation,
radioactive decay, and production may affect the contaminantconcentration. They can all be included in the sink/source term,
Rs , in Eq. (2). The transport problem can still be solved analytically if this term is described in terms of linear processes. For
example, for relatively simple transport scenario's involving
one or several sets of zero- and first-order rate expressions, thegoverning equation can always be represented in the form:

∂C
∂2 C
∂C
= Dx
−u
− µC + γ ,
2
∂t
∂x
∂x

(8)

where µ is a general first-order decay rate (T-1) and γ a zeroorder production term (ML-3T-1). We assume in this study that µand γ are either zero or always positive. Additional processes
such as nonlinear exchange, precipitation/dissolution, and sorption or chelation of constituents by moving sediment, are notconsidered here since the resulting transport problem generally
will require numerical solution techniques.
Concentration modes and boundary conditions

A variety of solute source or sink terms may need to be implemented in the ADE. One possible process during transport inrivers is sorption of contaminants on sediment along the main
channel. The source/sink term for sorption may then be written
as (e.g., Bencala, 1983)

Rs = ρ s

S = KdC(4)

where S is the sorbed concentration expressed as mass of soluteper mass of sediment readily accessible for sorption (MM-1),
and ρs is the mass of sediment per unit volume of river water(ML-3). For binary exchange in a system with a constant totalsolute concentration, linear equilibrium exchange is given by

The mathematical solution of Eq. (8) and related transport
equations requires auxiliary conditions specifying the initialand boundary condition of the transport problem. It is importantthat the boundary conditions accurately describe the adopted
modes of solute injection and detection. Concentrations areusually expressed in terms of the amount of solute per unit
volume of water. Although the microscopic concentration isbased on an infinitesimally small liquid volume, for practical
purposes a much larger finite averaging scale must be used. In
this way a macroscopic volume-averaged or resident concentration, CV, can be defined as
V

C =

1
∫ ∫ ∫ c dV ,
ΔV

(9)

147

Martinus Th. van Genuchten, Feike J. Leij, Todd H. Skaggs, Nobuo Toride, Scott A. Bradford, Elizabeth M. Pontedeiro

where c is the local concentration (ML-3) in a microscopic
volume element, V, and ΔV is the representative elementaryvolume (Bear and Verruijt, 1987). Spatial or temporal averages
beyond the REV scale may also be employed (Leij and Toride,1995). Unless stated otherwise, we will assume that all concentrations are macroscopic and are of the volume-averaged or
resident type.
Flux-averaged concentrations, CF, may arise when the concentration is determined as the ratio of the solute flux, Js, andthe water flux, u, passing through a given cross section during
an arbitrary time interval (Kreft and Zuber, 1978; Parker andvan Genuchten, 1984a), i.e.,F
C = J s / u.

(10)

A relationship between the volume- and flux-averaged concentrations can be established by substituting the one-dimensional
advective-dispersive flux (Eq. (1)) into (10) to giveV
D x ∂C .
V
C =C −u ∂xFmδ (x − xo ),
A

(12)

(13)

where δ(x–xo) is the Dirac function (L-1) having the property∞∫ δ (x − xo )dx = 1-∞

with δ (x − xo ) = 0 for all x ≠ xo .

(14)

Proper selection of boundary conditions is a somewhat esoteric topic that has received considerable attention in the literature on transport in porous media (e.g., Batu et al., 2013; van
Genuchten and Parker, 1984). Many transport problems involve
the application of a solute with a concentration described by aspatially independent function, g(t), to a surface water bodythat, for modeling purposes, can have a finite, semi-infinite, or
infinite length. The method of application (or the contaminantsource) may also vary widely (pumping or spraying, line orpoint sources, continuous or finite pulse type applications).

148

(15)

Two types of conditions are often used at the inlet of semiinfinite media (0 ≤ x ≤ ∞) or finite media (0 ≤ x ≤ L). These
conditions are based on continuity of either the concentration or
the solute flux across the inlet boundary. We do not discuss
here the simultaneous use of both conditions since this situationis generally too complicated for deriving relatively simple
analytical solutions.
First- or concentration-type inlet boundary conditions, also
referred to as Dirichlet conditions, require the concentration to
be continuous across the interface at all times, i.e.,

C(0,t) = g(t)where f(x) is an arbitrary function versus distance, the simplest
case being the situation where the stream has a constant (zeroor positive) background concentration. An alternative initialcondition, often assumed for one-dimensional stream transport,
is the Dirac condition specifying instantaneous release of aprescribed mass. This relatively hypothetical condition assumes
that a contaminant mass, m, can be distributed instantaneously
over an infinitely small region, for example a thin plane of area
A across the main channel at some longitudinal position xo , i.e.,f (x) =

∂C
( ± ∞, t) = 0.
∂xt > 0.

(16)(11)

Notice that CV and CF become similar for small values of Dx/u,i.e., when dispersion is negligible as compared to advection.
We emphasize that CF is a mathematical entity that depends onthe formulation of the contaminant flux. As will be discussed
shortly, the selection of boundary conditions is related to theconcentration mode.
The initial condition can be formulated as

C(x, 0) = f (x),Boundary conditions are simplest for transport problems defined over infinite domains (–∞ < x < ∞). To ensure that theconcentration remains bounded, the inlet and outlet conditions
can be written in the form

The drawback of this condition is that the (macroscopic) concentration at the interface just inside the river system will, in
reality, not respond instantaneously to changes in the influent
concentration. Condition (16) generally leads to discrepancies
in the mass balance when considered over a finite or semiinfinite system (Batu and van Genuchten, 1990; Parker and van
Genuchten, 1984a). Mass conservation can be ensured by using
a third- or flux-type condition at the inlet (often referred to also
as a Robin or Cauchy type condition):⎛∂C ⎞⎜⎝ uC − D x ∂x ⎟⎠ + = ug(t),
x=0

(17)

where 0+ indicates a position just inside the system being considered. Analytical solutions for first-type conditions are generally somewhat simpler than those for third-type inlet conditions. This probably explains the popularity of first-type conditions despite mass balance errors in the predicted resident concentration distributions.
Some commonly-used functions for the influent concentration, g(t), in (16) or (17) are the Dirac function, a finite squarepulse, and a Heaviside or step function. An instantaneous application of a solute amount m at the inlet, x = 0, at an arbitrary
time, to, across an area A is given byg(t) =mδ (t − to ),
uA

(18)

where δ(t–to) is the Dirac delta function in time (T-1). Instantaneous release into a stream at x = 0 can be described alternatively also with a Dirac type initial condition.
The application of a finite pulse can be formulated as⎧⎪Cg(t) = ⎨ o⎩⎪0

0 < t ≤ tot > to(19)

where to is the time duration of the applied pulse having concentration Co. Note that to approaches infinity for a single stepchange. Because of the linearity of the ADE, solutions forproblems involving an arbitrary number of finite pulses ofdifferent durations (to) and strengths (Co) can be readily inferredExact analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equationfrom the single step solution; the Heaviside function is convenient for this purpose.
The exit boundary condition can be defined in terms of a zero gradient at a finite or infinite distance from the inlet, leading
to the assumption of having a finite or semi-infinite medium,
respectively. The infinite outlet condition is mathematically far
more convenient than the finite condition, and hence has been
far more popular in analytical transport studies, even whenfinite systems are simulated. The infinite outlet condition is
given by

∂C F
∂2 C F
∂C F
= Dx
−
u
∂t
∂x
∂x 2

(20)

The use of this condition for a finite system implies that a semiinfinite fictitious medium with identical properties is placed
beyond the system of interest. No assumptions have to be maderegarding the physical processes at the outlet, x = L, other thanthat the ADE is applicable just like everywhere else at x > 0.
The domain for most transport problems involving rivers isbest described with a semi-infinite or infinite system, rather
than a finite system. However, we shall provide in this paper
also several analytical solutions for finite systems since numerical solutions, by their very nature, can be applied only to finite
domains. Assuming that no dispersion occurs for x > L, theapproximate condition for a finite system is∂C(L,t) = 0.
∂x

(21)

This conditions implies both flux and concentration continuity
at the downstream boundary of the finite transport domain. A
zero-gradient outlet condition is also appealing when notransport occurs across a boundary, such as is the case for solute diffusion or heat transport problems involving insulated
media.
The choice of the boundary conditions is affected by theconcentration mode (Kreft and Zuber, 1978). In practice only
the dependency of the inlet condition is considered (cf. Parker
and van Genuchten, 1984a). A typical problem in terms of theusual resident concentration can be formulated as

∂CV
∂2 CV
∂CV
= Dx
−
u
∂t
∂x
∂x 2

(22)

subject to a uniform initial conditionV

C (x,0) = f (x) = Ci

(23)a third-type inlet boundary condition to avoid mass balance
problems⎛ V∂CV ⎞= ug(t)⎜ uC − D x⎟∂x ⎠ +⎝
x=0(24)

and a semi-infinite or finite outlet condition∂CV(L,t) = 0.
∂x(25a, b)(26)

subject toF∂C(∞,t) = 0.
∂x∂CV(∞,t) = 0,∂x

As shown by Parker and van Genuchten (1984b) this problem
can be readily transformed to the flux-averaged mode using Eq.
(11) to giveC (x,0) = Ci

(27)

C F (0,t) = g(t),

(28)

∂C F
(∞,t) = 0∂xorD ∂2 CV∂C F(L,t) = − x(L,t).
∂xu ∂x 2(29)

Notice that the mathematical problem in terms of CF involves
the simpler first-type inlet condition with mass now being conserved, unlike the use of a first-type condition for transport
problems in terms of CV. This shows that solutions for a semiinfinite system involving a first-type inlet condition represent
mass-conserving solutions if the concentration, inside andoutside the solution domain, is interpreted as being of the fluxaveraged type (Batu et al., 2013; Parker and van Genuchten,
1984a). The transformations leading to Eqs (26) through (29)are less convenient for finite systems as well as for semiinfinite systems having nonuniform initial conditions (Toride etal., 1993b).
Differences between the preferred third-type and first-type
solutions for CV are usually minor except for small values of thedimensionless distance (ux/Dx) and time (u2t/Dx) (e.g., vanGenuchten and Parker, 1984). Since advective transport usually
dominates dispersive transport in rivers and streams, the selection of a first- or third-type condition is of secondary importance for most practical river transport problems, except
near the inlet. However, the issue should not be overlooked forcode verification when small errors in the mass balance areimportant. Hence, analytical solutions will be presented forboth first- and third-type inlet conditions. All concentrations inthis study are understood to be of the resident type; fluxaveraged concentrations can be derived as needed from theresident solutions using Eq. (11).
SOLUTION TECHNIQUES
Problems in one dimension
Before presenting specific analytical solutions, two different
solution techniques will be briefly reviewed. Although a wide
variety of solution procedures can be followed, many linear
transport problems can be solved by either using a transformation of variables, or implementing Laplace transforms.
Most streams can be considered infinite systems in the longitudinal direction. This makes them quite suitable for making atransformation of variables to cast the ADE in a simpler form.
For example, the coordinate transformation

ξ = x − ut
τ =t

(30)

149

Martinus Th. van Genuchten, Feike J. Leij, Todd H. Skaggs, Nobuo Toride, Scott A. Bradford, Elizabeth M. Pontedeirotransforms Eq. (3) subject to⎪⎧ C x < 0
C(x,0) = ⎨ o⎩⎪ 0 x > 0C=∂C(±∞,t) = 0
∂x

(31)

into the heat transfer (or solute diffusion) problem

∂C
∂2 C
= Dx
2
∂τ
∂ξ

(32)

⎧⎪ C ξ < 0C(ξ ,0) = ⎨ o⎩⎪ 0 ξ > 0(33)

Analytical solutions for this and many related problems havebeen widely reported in the literature (e.g. Carslaw and Jaeger,
1959; Crank, 1975).
An alternative transformation

ξ=x − ut

(34)

4Dx t

allows the ADE to be written as
2
d C + 2ξ dC = 0,
2
dξ
dξ

(35)

⎪⎧ C ξ < −∞C(ξ ) = ⎨ o⎩⎪ 0 ξ > ∞(36)

This ordinary differential equation can be solved through areduction of order according todψ+ 2ξψ = 0dξwhere ψ =

dC
.
dξ

(37)

Separation of variables and integration leads to the following
solution of (37):ψ = α1 exp(−ξ 2 ),(38)

where α1 is an integration constant. Integration of ψ between ξ(corresponding to the desired concentration, C), and ∞, andmaking use of a zero-concentration as dictated by the infinite
boundary condition, leads to∞

C = ∫ α1 exp(−w2 )dw + α 2 =ξ⎛ x − ut ⎞π
α1erfc ⎜⎟ + α2⎜⎝ 4D t ⎟⎠2(39)

⎡ x − ut ⎤Co⎥.
erfc ⎢
2⎢⎣ 4Dx t ⎥⎦More complicated transformations were previously used by
Brenner (1962) and Selim and Mansell (1976); several can also
be found in Zwillinger (1989).
The complementary error function (erfc) in Eq. (40) frequently occurs in analytical solutions of the ADE. Several
approximations can be found in Gautschi (1964). Its evaluation
is sometimes not straightforward, especially when erfc occurs
in combination with the exponential function and either or bothof these functions have unusually large or small arguments. A
useful routine for evaluating erfc is listed in van Genuchten andAlves (1982). In some cases, especially for situations withexponentially decaying inflow concentrations, the argument of
the erfc function and other terms in the analytical solutions may
involve complex variables with real and imaginary parts. An
accurate subroutine for such cases is given by van Genuchten(1985b).
Linear partial differential equations are often solved with integral transforms. The transformations generally lead to ordinary differential or even algebraic equations which are mucheasier to solve than the original partial differential equations.
Subsequent inversion of the solution to the regular space-time
domain is generally the most complicated part of the solution
procedure. The inversion often can be accomplished in several
ways (cf. Özişik, 1989; Pérez-Guerrero et al., 2009; Sneddon,
1995; Spiegel, 1965). Laplace transforms traditionally havebeen the most popular for solving the ADE (Leij et al., 1993;
Toride et al., 1993a; van Genuchten, 1981).
Problems in several dimensions
Solutions of transport problems in two and three dimensions
are in many cases a direct extension of the one-dimensional
solutions. We consider here three-dimensional problems in
Cartesian coordinates (x,y,z) with uniform flow and dispersion
in the longitudinal direction and dispersion in the transverse
direction. We note that solutions for two-dimensional problems
follow directly from solution of the more general threedimensional examples.
The three-dimensional form of the ADE is given by

∂C
∂2 C
∂C
∂2 C
∂2 C
= Dx 2 − u
+ Dy
+
D
,
z
∂t
∂x
∂x
∂ y2
∂z 2

(41)

where Dy and Dz are dispersion coefficients (L2T-1) in the transverse y and z directions. Although the boundary conditions mayvary depending upon the assumed geometry of the system, theynormally include zero-gradient conditions at the transverse
boundaries. Below are the conditions for a somewhat idealized
three-dimensional system involving a semi-infinite longitudinal
direction and two infinite transverse directions:xin which w is a dummy integration variable, α2 is an additional
integration constant, and erfc represents the complementary
error function (e.g., Gautschi, 1964). The boundary conditions
are used to evaluate α1 and α2. The final solution isC(x, y, z,0) = f (x, y, z),⎛D x ∂C ⎞⎜⎝ C − ω u ∂x ⎟⎠ + = g( y, z,t)
x=0
∂C(∞, y, z,t) = 0,
∂x

150

(40)

(42)

ω = 0,1,

(43)

(44)

Exact analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equation

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

(45)

(46)

where ω equals 0 or 1 for a first- and third-type boundary condition, respectively. Note that f(x,y,z) and g(y,z,t) are arbitrary
initial and influent distributions, respectively. These functions
can be specified in a similar manner as for the one-dimensionalcase except that f and g now depend on multiple coordinates
instead of one.
For three-dimensional infinite systems the use of a moving
longitudinal coordinate (Eq. (30)) readily transforms the ADE
and its auxiliary conditions into a three-dimensional diffusion
problem for which a wealth of analytical solutions exist. Another straightforward method of solving multidimensionaltransport problems is the use of the product rule which statesthat the product of the solutions of (simpler) one-dimensional
problems may directly give the desired three-dimensional solution (e.g., Carslaw and Jaeger, 1959). Consider the three onedimensional problems
2

∂C x
∂ Cx
∂C
= Dx
−u x ,
2
∂t
∂x
∂x

∂C y
∂t

= Dy

∂2 C y
∂ y2

,

∂C z
∂2 Cz
= Dz
.
∂t
∂z 2

(47)

(48)

(49)

The product of the solutions of these three one-dimensional
problems, i.e.,C(x, y, z,t) = C x (x,t) C y ( y,t) C z (z,t)(50)

also satisfies the three-dimensional transport equation. Application of the product rule requires that the boundary conditions behomogeneous and that the initial condition for the threedimensional problem can be written as a product of the three
one-dimensional initial conditions.
Several approaches, in addition to moving coordinate transformation and product rules, have been used in the literature.
These include separation of variables and Fourier series (Bruchand Street, 1967), Green's functions (Leij and van Genuchten,
2000; Sagar, 1982; Yeh and Tsai, 1976), or integral transforms(Cleary, 1973; Leij and Dane, 1990; Leij et al., 1993, Pérez
Guerrero et al., 2013).
For many river systems at least one direction will be finite,which may complicate the mathematical solutions. Several
integral transforms, mostly of the Fourier type, have been usedby Özişik (1989) and Sneddon (1995) for such situations. In
some cases one or two of the finite directions can be neglectedby using a depth-averaged concentration, or by making use of
uniformity in the initial and boundary conditions across thestream. Another approach is the use of fictitious constituent
sources outside the solution domain to ensure a zero-gradient at
the boundaries (e.g., see p. 49 of Fischer et al., 1979). Thisapproach permits one to still use the relatively simple solutions
available for infinite outlet conditions.
SPECIFIC ADE SOLUTIONS
One-dimensional solutions with third-type inlet condition
Consider the ADE for a contaminant subject to a combination of first-order decay and zero-order production to represent
a variety of physical, chemical or biological reactions:

∂C
∂2 C
∂C
= Dx 2 − u− µC + γ .
∂t
∂x
∂x

(51)

Eq. (51) will be solved for the following initial and boundary
conditions:C(x,0) = f (x),

(52)

⎛
∂C ⎞
⎜⎝ uC − ω Dx ∂x ⎟⎠ + = ug(t),
x=0

(53)

∂C
(∞,t) = 0
∂x

or

∂C
(L,t) = 0.
∂x(54a, b)Analytical solutions for a third-type inlet condition (ω = 1)
are given below under Case A, while those for a first-type inlet
condition (ω = 0) are presented under Case B. The first solution
will be for a Dirac delta type initial distribution, from which
solutions for arbitrary initial profiles can be derived. The remaining one-dimensional solutions were selected from van
Genuchten (1981) and van Genuchten and Alves (1982). The
latter publication contains a large number of solutions for situations with and without zero-order production and first-orderdecay. Most of the semi-infinite solutions, as well as others,
were later included in the CXTFIT computer programs (Parker
and van Genuchten, 1984b; Toride et al., 1999), which later
were incorporated into the windows-based STANMOD software package (Simůnek et al., 2000).
Case A0. Semi-infinite domain with zero initial concentration and a Dirac-type input function g(t) = mδ(t)/Au. Neglecting
any production or decay (µ = γ = 0), the solution isC(x,t) =⎡ (x − ut)2 ⎤m ⎧⎪ 1
exp ⎢ −⎥⎨A ⎪ πD t
4Dx t ⎥⎦⎢⎣x⎩
⎡ x + ut ⎤ ⎫⎪⎛ ux ⎞
u
⎢
⎥ ⎬.
−
exp ⎜
erfc
⎟
2Dx
⎢⎣ 4Dx t ⎥⎦ ⎪⎝ Dx ⎠⎭(55)

As mentioned previously, an alternative solution exists for a
Dirac-type initial distribution in an infinite medium.
In case of an arbitrary input function, g(t), the expression forthe concentration can be written as the convolution of the input
signal and the concentration for the Dirac function:

151

Martinus Th. van Genuchten, Feike J. Leij, Todd H. Skaggs, Nobuo Toride, Scott A. Bradford, Elizabeth M. Pontedeiro⎧⎪t⎡ (x − uτ )2 ⎤1
C(x,t) = ∫ ug(t − τ ) ⎨exp ⎢ −⎥
0⎢⎣ 4Dxτ ⎥⎦⎪⎩ π Dxτ
⎡ x + uτ ⎤ ⎫⎪⎛ ux ⎞u⎢⎥ ⎬ dτ .
−
exp ⎜
erfc
⎟
2Dx
⎢⎣ 4Dxτ ⎥⎦ ⎪⎝ Dx ⎠⎭⎡ x − ut ⎤⎡ (x − ut )2 ⎤
1
u 2t
⎥+
⎥
erfc ⎢
exp ⎢ −2
π Dx⎢⎣ 4Dx t ⎥⎦⎢⎣ 4Dx t ⎥⎦⎡ x + ut ⎤⎛ ux ⎞
1⎛ux u 2t ⎞⎥.
− ⎜1 ++
erfc ⎢⎟ exp ⎜⎟
2⎝
Dx Dx ⎠⎢⎣ 4Dx t ⎥⎦⎝ Dx ⎠

A(x,t) =(56)

Fig. 1 shows an example of the Dirac solution given by Eq.
(55) using parameter values selected previously by De Smedt etal. (2005) to test their analytical solution for the same Dirac
inlet problem but assuming a first-type boundary condition (thesolution of which is given by Case B0 below). The example
involves an instantaneous injection (or spill) of 1 kg of a solutein the main channel of a river having a cross-section of 10 m2,
an average flow velocity of 1 m/s and a dispersion coefficientof 5 m2/s. De Smedt et al. (2005) in their analysis also included
the effects of a dead storage zone which we ignore here (thatcase will be discussed in part 2). Fig. 1 shows the solute distributions in the stream after 500, 1000 and 1500 s. The curve fort = 1000 s is essentially the same as the one calculated by De
Smedt et al. (2005) assuming no mass transfer into the deadzone (their Fig. 1 with α = 0).

(59)

Here we illustrate one possible application of the above analytical solution by fitting the transport parameters in Eqs. (58)
and (59) to observed data obtained by Brevis et al. (2001) as
discussed by De Smedt et al. (2005). The data pertain to a tracer
experiment in the Chillán River in Chile using a 20% Rhodamine WT solution injected simultaneously at multiple points
along a lateral cross-section of the river at x = 0. Fig. 2 shows
the observed data of one of their experiments (Exp. I-3), in thiscase for the application of 157.1 g of tracer at x = 0, with measurements made 4604 m downstream of the injection point.
Because of the relative short duration of the experiments, photodegradation and sorption of Rhodamine WT were presumed
to be negligible (De Smedt et al., 2005).
7

3Concentration (mg/m )

6
5
4
3
2
1
0
5000

7500

10000

12500

15000

17500

20000

Time (s)Fig. 1. Calculated concentration distributions obtained with Eq.
(55) for instantaneous injection of 1 kg of solute in a stream having
a cross-sectional area of 10 m2 assuming a third-type Dirac input
boundary condition (u = 1 m s-1, Dx = 5 m2 s-1).

Case A1. Semi-infinite domain with uniform initial concentration, f(x) = Ci and no production or decay. The inlet concentration function, g(t), is of the pulse type (a Heaviside type function) with constant concentration Co , i.e.,⎧⎪Cg(t) = ⎨ o⎪⎩0

0 < t ≤ tot > to .

(57)

The solution of this example is⎧Ci + (Co − Ci ) A(x,t)
0 < t ≤ to⎪C(x,t) = ⎨⎪ C + (C − C ) A(x,t) − C A(x,t − t ) t > t ,
o
i
o
o
o
⎩ i
where

152

(58)

Fig. 2. Observed (solid squares) and fitted (continuous line) concentrations for tracer experiment I-3 of Brevis et al. (2001) and De
Smedt et al. (2005).

Fig. 2 shows the data along with the fitted curve based on
Eqs (58) and (59). Parameters were estimated using the nonlinear least-squares optimization features of the CXTFIT code
(Toride et al., 1999) within STANMOD. The code allows one
to estimate not only the longitudinal flow velocity u and thelongitudinal dispersion coefficient Dx, but also the concentration Co of the applied tracer solution for a given value of theinjection or pulse time to in Eq. (58). Assuming a very short
injection period of only 10 s, we obtained the following parameter values (with their 95% confidence intervals): u = 0.426 ±0.003 m/s, Dx = 11.4 ± 1.4 m2/s, and Co = 1729 ± 86 mg/m3.
The coefficient of determination (R2) of the fit was 0.983, andthe root mean square error (RMSE) 0.382 mg/m3). One may
verify that the total amount of solute mass (m) injected per m2
cross-sectional area hence equals m = u Co to, or 7.37 g/m2.
Given that a total amount of 157.1 g of tracer was applied to theriver, this translates to an effective cross-sectional area of 157.1
g/7.37 (g/m2) or 21.3 m2 of the channel as calculated with theequilibrium ADE transport model.
The above calculations assume applicability of Eqs. (58) and(59) using a solute injection pulse of 10 s. Exactly the same
results were obtained when the total mass was assumed to be
applied instantaneously (as estimated with CXTFIT using theExact analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equation

Dirac solution given by Eq. (55)), or spread over a time period
to of 100 s. We also ran the problem assuming a first-type inlet
boundary condition (Cases B0 and B1 given later in this paper).
Results were again the same or very close to those obtainedwith the third-type boundary conditions. For example, for the
first-type boundary condition (Case B0) we obtained u = 0.424± 0.003 m/s, Dx = 11.3 ± 1.3 m2/s, and Co = 1729 ± 860 mg/m3(R2 = 0.970; RMSE = 0.382), which essentially duplicates theresults for the third-type solutions. This shows that differences
between the first- and third-type boundary conditions fortransport problems with relatively large values of the system
Peclet number, uL/Dx (van Genuchten and Parker, 1984) generally are very small. The calculated curve in Fig. 2 assuming
equilibrium transport provided a reasonable fit of the data,albeit not perfect. In part 2 (van Genuchten et al., 2013) we will
analyze the same problem in terms of the transient storage
model, similarly as was considered by De Smedt et al. (2005).
Case A2. The same as case A1 (i.e., for a semi-infinite domain with uniform initial concentration Ci and a pulse-type
input function with concentration Co), but now with zero-orderproduction and first-order decay. The solution is⎧⎛γ⎞
0 < t ≤ to⎪ ⎜ Co − ⎟ A(x,t) + B(x,t)µ⎠⎪⎝⎪
C(x,t) = ⎨⎪⎪⎛ C − γ ⎞ A(x,t) + B(x,t)⎪⎩⎜⎝ o µ ⎟⎠− Co A(x,t − to ) t > to ,assuming that the solute was injected during a finite (pulse)
period to of 1000 s, and that the solute was subject to first-orderdecay (µ ≥ 0), but without any first-order production (γ = 0).
Case A3. Steady-state solution for a semi-infinite domain
with a third-type inlet concentration Co, and zero- and firstorder decay. The analytical solution can be written as

C(x) =⎡ (u − ξ )x ⎤γ ⎛γ ⎞ ⎛ 2u ⎞
+ C −
exp ⎢
⎥.
µ ⎜⎝ o µ ⎟⎠ ⎜⎝ u + ξ ⎟⎠⎣ 2Dx ⎦

(64)

Case A4. Semi-infinite domain with a nonuniform initialconcentration, f(x)=E(x), and a pulse-type inlet condition like
Case A1 with g(t) given by Eq. (57). As initial condition we use
the steady-state profile resulting from continuous application ofa solute at background concentration, Cb. This initial profile
follows directly from Case A3 as:C(x,0) ≡ E(x) =⎡ (u − ξ )x ⎤γ ⎛γ ⎞ ⎛ 2u ⎞
+ C −
exp ⎢
⎥.
µ ⎜⎝ b µ ⎟⎠ ⎜⎝ u + ξ ⎟⎠⎣ 2Dx ⎦

(65)

(60)

where

⎡ x − ξt ⎤⎡ (u − ξ )x ⎤
u
⎥
exp ⎢
⎥ erfc ⎢u+ξ
⎢⎣ 4Dx t ⎥⎦⎣ 2Dx ⎦
⎡ x + ξt ⎤⎡ (u + ξ )x ⎤
u
⎥
+
exp ⎢
⎥ erfc ⎢u−ξ⎢⎣ 4Dx t ⎥⎦⎣ 2Dx ⎦
⎡ x + ut ⎤
⎛ ux
⎞
u2
⎥.
+
exp ⎜
− µt ⎟ erfc ⎢2 µ Dx⎢⎣ 4Dx t ⎥⎦⎝ Dx⎠

A(x,t) =B(x,t) =(61)⎧ 1⎡ x − ut ⎤
γ ⎛
γ⎞
⎪
⎥
+ ⎜ Ci − ⎟ exp ( − µt ) ⎨1− erfc ⎢µ ⎝µ⎠⎢⎣ 4Dx t ⎥⎦
⎪⎩ 2
−
+

The analytical solution for this case is⎧ (Co − Cb ) A(x,t) + E(x)0 < t ≤ to⎪C(x,t) = ⎨⎪ (C − C ) A(x,t) + E(x) − C A(x,t − t ) t > t , (66)
b
o
o
o
⎩ o

⎡ (x − ut)2 ⎤u 2t
exp ⎢ −⎥π Dx
4Dx t ⎥⎦⎢⎣

⎡ x + ut ⎤ ⎫⎪⎛ ux ⎞
1⎛ux u 2t ⎞⎥ ⎬.
+erfc ⎢⎜ 1+⎟ exp ⎜⎟
2⎝
Dx Dx ⎠⎢⎣ 4Dx t ⎥⎦ ⎪⎝ Dx ⎠⎭(62)

in whichξ = u 1+ 4 µ Dx / u 2 .

Fig. 3. Effect of first-order degradation (µ) on calculated concentration distributions in a stream after injection of a solute having a
concentration of 10 g/m3 for 1000 s (u = 1 m s-1, Dx = 5 m2 s-1, γ =0 g m-3 s-1, to = 1000 s).

(63)

Fig. 3 shows an application of this solution using the sameproblem and parameter values as for Case A0 (Fig. 1), but now

where A(x,t) is given by Eq. (61) and E(x) by (65).
Case A5. Semi-infinite domain with uniform initial concentration, f(x)=Ci, no production and decay, and an influent concentration that decreases exponentially with time:g(t) = Cb + Co exp(− λ t)(67)

where λ is a constant (T-1). The solution for this case isC(x,t) = Ci + (Cb − Ci ) A(x,t) + Co B(x,t),

(68)

153

Martinus Th. van Genuchten, Feike J. Leij, Todd H. Skaggs, Nobuo Toride, Scott A. Bradford, Elizabeth M. Pontedeiro

where

the following approximate solution may be used (Brenner,
1962; van Genuchten and Alves, 1982):⎡ x − ut ⎤⎡ (x − ut)2 ⎤
1
u 2t
⎥+
erfc ⎢
exp ⎢ −⎥
2
π Dx
4Dx t ⎥⎦⎢⎣ 4Dx t ⎥⎦⎢⎣⎡ x + ut ⎤⎛ ux ⎞
1⎛ux u 2t ⎞⎥.
− ⎜ 1++
erfc ⎢⎟ exp ⎜⎟
2⎝
Dx Dx ⎠⎢⎣ 4Dx t ⎥⎦⎝ Dx ⎠

A(x,t) =⎡ x − ut ⎤⎡ (x − ut)2 ⎤C(x,t) 1
u 2t⎥+⎥= erfc ⎢exp ⎢ −
Co2
π Dx⎢⎣ 4Dx t ⎥⎦⎢⎣ 4Dx t ⎥⎦⎡ x + ut ⎤⎛ ux ⎞
1⎛ux u 2t ⎞⎥− ⎜ 1++
erfc ⎢⎟ exp ⎜⎟
2⎝
Dx Dx ⎠⎢⎣ 4Dx t ⎥⎦⎝ Dx ⎠(69)⎧ u⎡ x −κ t ⎤⎛ (u − κ )x ⎞⎪⎥B(x,t) = exp(− λ t) ⎨exp ⎜erfc ⎢⎟u +κ
⎢⎣ 4Dx t ⎥⎦⎝ 2Dx ⎠⎩⎪+⎡ uL (2L − x + ut)2 ⎤⎤4u 2t ⎡u(2L − x + ut) ⎥ exp ⎢−⎥⎢1+
π Dx ⎣ 4Dx
4Dx t
⎢⎣ Dx
⎥⎦
⎦

⎡ x + κ t ⎤ ⎫⎪⎛ (u + κ )x ⎞
u
⎥⎬
+
exp ⎜
erfc ⎢⎟u −κ⎢⎣ 4Dx t ⎥⎦ ⎪⎝ 2Dx ⎠⎭
⎡ x + ut ⎤⎛ ux ⎞
u2⎥−
exp ⎜
erfc ⎢⎟
2λ Dx⎢⎣ 4Dx t ⎥⎦⎝ Dx ⎠•(70)

κ = u 1− 4λ Dx / u 2 .

(71)

Note that this example reduces to Case A1 (only for t < to)when Co = 0.
Case A6. Finite domain with zero initial concentration, f(x)= 0, a third-type inlet condition, and no zero- and first-order
rate processes (i.e., µ = γ = 0). The inlet condition is again thesame as for Case A1, i.e., a pulse type concentration with constant Co. The outlet condition is now a zero gradient at x = L as
given by Eq. (21). The solution for this case is (Brenner, 1962):2⎡ ux
u 2t β m Dx t ⎤G(x)exp ⎢
−
+
⎥
L2 ⎥⎦
⎢⎣ 2Dx 4DxC(x,t)
2uL ∞= 1−, (72a)∑22⎤
Co
Dx m−1 ⎡⎤⎡⎛⎞⎛⎞uL ⎥ ⎢ 2uL ⎥⎢ β 2 + uL+
β +⎢ m ⎜⎝ 2Dx ⎟⎠Dx ⎥ ⎢ m ⎜⎝ 2Dx ⎟⎠ ⎥
⎣
⎦⎣
⎦

(75)

One-dimensional solutions with first-type inlet condition
Case B0. Semi-infinite domain with zero initial concentration, a first-type (ω = 0) inlet condition, and a Dirac-type input
concentration. The solution for the case of no production or
decay isC(x,t) =m
Aut

C(x,t) = ∫0⎡⎛ β x ⎞ uL
⎛ β x⎞ ⎤G(x) = β m ⎢ β m cos ⎜ m ⎟ +sin ⎜ m ⎟ ⎥⎝ L ⎠ 2Dx ⎝ L ⎠ ⎥⎦
⎢⎣

(72b)

in which the eigenvalues βm are the positive roots of2

β mDx uL
+
= 0.
uL
4Dx

(73)

or

uL
> 100
Dx

(76)

⎡ (x − uτ )2 ⎤exp ⎢⎥ dτ .
⎢⎣ 4Dxτ ⎥⎦4π Dxτ 3

xg(t − τ )(77)Case B1. Semi-infinite domain with uniform initial concentration, f(x) = Ci, and a pulse-type solute application (i.e., Eq.
(57) as for Case A1). The solution is⎧Ci + (Co − Ci ) A(x,t)
0 < t ≤ to⎪C(x,t) = ⎨⎪⎩Ci + (Co − Ci ) A(x,t) − Co A(x,t − to ) t > to

(78)

where

The above series solution converges very slowly when advection dominates dispersion as is often the case for stream
transport problems. For eitheruLut> 5 + 40Dx
L

⎡ (x − ut)2 ⎤exp ⎢ −⎥.
⎢⎣ 4Dx t ⎥⎦4π Dx t 3x

For arbitrary input profiles the expression for the concentration
can again be written as the convolution of the input signal andthe concentration for the Dirac input:

where

154

⎡ (2L − x) + ut ⎤ ⎫⎪⎛ uL ⎞⎥ ⎬.
exp ⎜erfc ⎢⎟⎢⎣⎝ Dx ⎠
4Dx t ⎥⎦ ⎪⎭Although in most cases considerably more complicated, alarge number of analytical solutions are also available for finite
systems with zero-order production and/or first order decay(Pérez Guerrero et al., 2009; van Genuchten and Alves, 1982).

andβ m cot( β m) −u ⎪⎧ ⎡
3utu
2⎤+(2L − x + ut ) ⎥⎨ ⎢ 2L − x +Dx ⎩⎪ ⎣
2 4Dx
⎦

−

(74)

A(x,t) =⎡ x + ut ⎤⎛ ux ⎞⎡ x − ut ⎤ 1
1⎢⎥.
erfc ⎢
+
exp
erfc
⎥
⎜D ⎟
2
⎢⎣ 4Dx t ⎥⎦⎝ x ⎠⎣ 4Dxt ⎦ 2(79)Case B2. Semi-infinite domain with uniform initial concentration, f(x) = Ci, and a pulse-type solute application. The solution for a first-type inlet condition isExact analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equation⎧⎛γ⎞
0 < t ≤ to⎪ ⎜ Co − ⎟ A(x,t) + B(x,t)µ⎠⎪⎝⎪
C(x,t) = ⎨⎪⎪⎛ C − γ ⎞ A(x,t) + B(x,t)⎪⎩⎜⎝ o µ ⎟⎠− Co A(x,t − to ) t > to ,and decay, and an influent concentration that decreases exponentially with time as defined for Case A5 (Eq. (67)). The
solution for this case is(80)

A(x,t) =⎡ x − ξt ⎤⎡ (u − ξ )x ⎤
1
⎥
exp ⎢
⎥ erfc ⎢2⎢⎣ 4Dx t ⎥⎦⎣ 2Dx ⎦
⎡ x + ξt ⎤⎡ (u + ξ )x ⎤1⎥
+ exp ⎢
⎥ erfc ⎢2⎢⎣ 4Dx t ⎥⎦⎣ 2Dx ⎦

A(x,t) =

(81)

and

B(x,t) =γ ⎛γ⎞+ ⎜ Ci − ⎟ exp(− µt)
µ ⎝
µ⎠

⎧ 1
⎡ x − ut ⎤ 1⎡ x + ut ⎤ ⎫⎪⎛ ux ⎞
⎪
⎢
⎥
⎢
⎥⎬
1−
erfc
−
exp
erfc
⎨
⎜D ⎟
⎢⎣ 4Dx t ⎥⎦ 2⎢⎣ 4Dx t ⎥⎦ ⎪
⎝ x⎠
⎪⎩ 2
⎭

(82)

in which ξ is as defined previously by Eq. (63).
Case B3. Steady-state solution for a semi-infinite domain
with uniform initial concentration, f(x) = Ci. The analyticalsolution for a step-type inlet condition with constant concentration Co is⎡ (u − ξ )x ⎤γ ⎛γ⎞
+ ⎜ Co − ⎟ exp ⎢
⎥.
µ ⎝
µ⎠
⎣ 2Dx ⎦

(83)

Case B4. Semi-infinite domain with nonuniform initial concentration, f(x) = E(x), and again a pulse-type inlet condition ofCase A1 with constant Co. As initial condition we use thesteady-state profile resulting from a continuous application of
the solute at background concentration, Cb. This initial profile
follows from Case B2 as

C(x,0) ≡ E(x) =⎡ (u − ξ )x ⎤
γ ⎛
γ⎞
+ C −
exp ⎢
⎥,
µ ⎜⎝ b µ ⎟⎠⎣ 4Dx ⎦

(84)

where Cb, as before, represents the background concentration
before the pulse of concentration Co is applied. The solution forthis case is

(

)

(

)

⎧ Co − Cb A(x,t) + E(x)0 < t ≤ to⎪⎪C(x,t) = ⎨⎪⎪⎩ Co − Cb A(x,t) + E(x)− Co A(x,t − to ) t > to ,

(86)

where

where

C(x) =

C(x,t) = Ci + (Cb − Ci ) A(x,t) + Co B(x,t),(85)

where A(x,t) is given by Eq. (81) and E(x) by (84).
Case B5. Semi-infinite domain with uniform initial concentration, f(x) = Ci, a first-type inlet condition with no production⎡ x − ut ⎤ 1⎡ x + ut ⎤⎛ ux ⎞
1⎥ + exp ⎜⎥,erfc ⎢
erfc ⎢⎟
2⎢⎣ 4Dx t ⎥⎦ 2⎢⎣ 4Dx t ⎥⎦⎝ Dx ⎠⎧1⎡ x −κ t ⎤⎛ (u − κ )x ⎞⎪⎥B(x,t) = exp(− λ t) ⎨ exp ⎜erfc ⎢⎟⎢⎣ 4Dx t ⎥⎦⎝ 2Dx ⎠⎪⎩ 2
⎡ x + κ t ⎤ ⎫⎪⎛ (u + κ )x ⎞1⎢⎥⎬
+ exp ⎜
erfc⎟
2⎢⎣ 4Dx t ⎥⎦ ⎪⎝ 2Dx ⎠
⎭

(87)

(88)

and κ as given by Eq. (71).
Case B6. Finite domain with uniform initial concentration,f(x) = 0, a first-type inlet condition, and no zero- or first-orderproduction or degradation processes. The inlet condition is thesame as for Case A1 (Eq. (57)), i.e., a pulse type application
with constant Co. The solution is (Cleary and Adrian, 1973):2⎡ ux⎛ β x⎞
u 2t β m Dx t ⎤2β m sin ⎜ m ⎟ exp ⎢
−
−
⎥
⎝ L ⎠
∞
L2 ⎥⎦
⎢⎣ 2Dx 4DxC(x,t)=1 − ∑
C0
m=1

2⎛ uL ⎞uLβ m2 + ⎜⎟ + 2D2D⎝x⎠
x, (89)

where the eigenvalues βm are the positive roots ofβ m cot( β m) +

uL
= 0.
2Dx

(90)

As for Case A6, this series solution converges very slowly
when advective transport dominates dispersion. For conditions
given by (74), the solution is accurately approximated by themuch more convenient (van Genuchten and Alves, 1982):⎡ x − ut ⎤ 1⎡ x + ut ⎤⎛ ux ⎞
C(x,t) 1⎥ + exp ⎜⎥= erfc ⎢erfc ⎢⎟
Co
2⎢⎣ 4Dx t ⎥⎦ 2⎢⎣ 4Dx t ⎥⎦⎝ Dx ⎠⎡ (2L − x) + ut ⎤⎛ uL ⎞1⎡u(2L − x) u 2t ⎤⎥+ ⎢2 ++
erfc ⎢⎥ exp ⎜
⎟
2 ⎢⎣
Dx
Dx ⎥⎦
⎢⎣
⎝ Dx ⎠
4Dx t ⎥⎦⎡ uL (2L − x + ut )2 ⎤
u 2t
⎥.
−
exp ⎢
−
π Dx
4Dx t
⎢⎣ Dx
⎥⎦

(91)

4.3 Two-dimensional solutions
Case C1. Infinite domain with Dirac-type initial concentration and no zero- and first-order rate processes. The longitudinal flow direction (x) and the transverse direction (y) are assumed to be infinite. An amount m is released instantaneously
over the depth the river (0 < z < d) at x = 0 and y = 0. The initialcondition for this line source problem can be written as155

Martinus Th. van Genuchten, Feike J. Leij, Todd H. Skaggs, Nobuo Toride, Scott A. Bradford, Elizabeth M. Pontedeiromδ (x, y)dg(x, y) =(92)

in which the Dirac delta function (δ) has the dimension of L-2.
The analytical solution for the concentration can be found withthe product rule asC(x, y,t) =⎡ (x − ut)2y2 ⎤exp ⎢ −−⎥.
4Dx t
4D y t ⎥⎦4π dt Dx D y
⎢⎣
m

(93)

Arbitrary initial profiles can be viewed as a superposition of
instantaneous releases and can therefore be obtained by integrating the analytical solution over time (see also Cases A0 and
B0).
Case C2. Infinite domain with a continuous line source at x
= 0 and y = 0 from which a nonreactive constituent is released
at a rate of Wo (MT-1). The steady-state solution is given by(Sayre, 1973)
C(x, y) =⎡ u⎛ ux ⎞
Wo
exp ⎜
K0 ⎢
⎟
⎢ 4Dx
⎝ 2Dx ⎠
2π d Dx D y⎣

x 2 y 2 ⎤⎥+,Dx D y ⎥⎦(94)

where K0[ ] is the zero-order modified Bessel function of thesecond kind. If longitudinal dispersion is negligible, the following approximate solution may be usedC(x, y) =if

y 2 Dx2

x Dy⎛ uy 2 ⎞exp ⎜ −⎟
2d π uD y x⎝ 4D y x ⎠g(x, y, z) = m δ (x, y, z)<< 1 and x >>C(x, y, z;t) =

∞

(

)

(

)

where B is the channel width (L), with transverse solution domain –B/2 ≤ y ≤ B/2, and n is the number of reflection cycles.
Usually, a value of 4 or 5 for n yields sufficiently accurate
results.
Three-dimensional solutions
The governing equation for three-dimensional transport issimilar as for two-dimensional transport except that a second
transverse dispersion term with coefficient Dz is included. The(98)

Solutions for arbitrary initial profiles can again be obtained byintegrating the point source solution over time.
Case D2. Semi-infinite domain with a constant water depth,d, for which a nonreactive material is continuously released at a
rate Wo (MT-1) from a point source at x = 0, y = yo, and z = zo. It
is also assumed that longitudinal dispersion is negligible (i.e.,
Dx = 0) as compared to advective transport. The steady-statesolution for this case is (Onishi et al., 1981).

(95)⎡C * x, nB − y + (−1) n y⎤o⎢⎥,C(x, y) = C (x, y − yo ) + ∑ ⎢⎥ (96)
n = 1 +C * x, − nB − y + (−1) n y
⎢⎣
⎥⎦
o
*

m
8π t π Dx D y Dz t

2 ⎤⎡ (x − ut)2
y2
exp ⎢ −−− z ⎥.
4Dx t
4D y t 4Dz t ⎥⎦⎢⎣C(x, y, z) =

2Dx
.
u

(97)

in which the Dirac delta function (δ) now has the dimension of
L-3. The expression for the concentration follows again from theproduct rule as (Onishi et al., 1981; Sayre, 1973)

Wo

Case C3. Semi-infinite longitudinal and finite transverse
domain having a continuous line source at x = 0 and y = yo (notnecessarily at y = 0) from which a contaminant is released at a
rate Wo (MT-1). The steady-state solution is very similar to thesolution of Case C2; however, part of the contaminant eventually will reach the banks of the river because of its finite width.
The mirror-image technique (e.g., Fischer et al., 1979) may be
used to ensure that the constituent beyond the river bank, as
predicted with the solution of Case C2 (denoted here as C*) for
the infinite transverse condition, is reflected back. The solution
can be written in terms of C* with modified y variable (Sayre,1973):

156

three-dimensional problem can be readily solved for an infinite
z direction.
Case D1. Infinite domain with Dirac-type initial concentration and a first-type inlet condition with no zero- and first orderprocesses. The x, y and z directions are infinite for all practical
purposes. An amount m is released instantaneously at x = y = z= 0. The initial condition for this point source can be written as⎧⎪⎛ ( y − y )2 u ⎞Woo⎟⎨exp ⎜ −4D
4π x D y Dz ⎪⎩⎝yx ⎠⎛ ( y + y )2 u ⎞ ⎫⎪o
+ exp ⎜ −⎟⎬4D y x ⎠ ⎪
⎝
⎭
⎧⎪

∞

•

⎛ (z − z − 2md)2 u ⎞
o
⎟
4Dz x
⎝
⎠

∑ ⎨exp ⎜ −m=−∞ ⎪⎩⎛ (z + z − 2md)2 u ⎞ ⎫⎪o
+ exp ⎜ −
⎟ ⎬.
4Dz x
⎝
⎠ ⎪⎭

(99)

Case D3. Semi-infinite domain with the solute initially uniformly distributed in a parallelepipedal region with no inflow ofmass at the upstream boundary. The y and z directions are assumed to be infinite. The problem can be solved with the product rule. The longitudinal transport equation, as given by Eq.
(56), is solved subject to⎧ C1/3⎪ i
C x (x,0) = ⎨⎪ 0⎩

x1 < x < x2

(100)

elsewhere

and⎛Dx ∂C x ⎞⎜⎝ C x − ω u ∂x ⎟⎠

x=0+

=0

∂C x
(∞, t) = 0,
∂x

(101)

where x1 and x2 are two arbitrary locations between which thesolute is originally present. The solution can be deduced direct-

Exact analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equation

ly from examples A5 and A6 in van Genuchten and Alves(1982). The solution for a first-type inlet condition (ω = 0) is
given byC x (x,t) =

Ci1/3
2

C(x, y, z, t) = C x (x,t) C y ( y,t) C z (z,t).

⎛ x − x − ut ⎞⎛ x − x − ut ⎞2
1
erfc ⎜⎟ − erfc ⎜⎟⎜⎝⎜⎝ 4D t ⎟⎠4Dx t ⎟⎠x

⎛ x + x + ut ⎞⎛ x + x + ut ⎞ ⎤⎛ ux ⎞ ⎡2
1⎢+ exp ⎜
erfc
−
erfc
⎜
⎟
⎜
⎟⎥
⎟
⎜⎝
⎟⎠
⎜⎝ 4D t ⎟⎠ ⎥⎝ Dx ⎠ ⎢
4D
t
x
x
⎣
⎦

(102)

whereas the solution for a third-type inlet condition (ω = 1) is

C x (x,t) =and for homogeneous boundary conditions can now be written
as the product of individual solutions as follows⎛ ux ⎞ ⎡⎛⎞Ci1/3 ⎧⎪u⎨ exp ⎜⎟ ⎢⎜ 1+ D (x + x1 + ut)⎟
2 ⎪D
⎝ x ⎠ ⎢⎣⎝⎠x⎩
⎛ x + x + ut ⎞1
• erfc ⎜⎟⎜⎝4Dx t ⎟⎠

Case D4. Semi-infinite domain with zero initial concentration, f(x) = Ci, and a third-type inlet condition. The inlet condition is of the pulse type with constant concentration Co, i.e.,⎧Co | y | < a, | z | < b, 0 < t ≤ to⎪g( y, z, t) = ⎨⎪0
otherwise.
⎩C(x, y, z,t) =⎛ ux ⎞ ⎡⎛ (x + x2 + ut) ⎞
4u 2t
+
exp ⎜
⎢exp ⎜ −
⎟
⎟
π Dx
4Dx t
⎝ Dx ⎠ ⎢⎣
⎝
⎠

•

∂C y
∂y

(104)

(105)

The same solution can be used to express Cz in the z-direction
where, for example, the solute is initially between –b and b.
The solution of the three-dimensional transport equation
given by (41) for solution domain {0 < x < ∞, –∞ < y < ∞, –∞ <z < ∞, t > 0}, subject to the initial condition

x1 < x < x2 | y | < a, | z | < belsewhere⎛ uτ − x ⎞γ t⎟∫ exp(− µτ )erfc ⎜⎜⎝ 4D τ ⎟⎠20

⎛ x + uτ ⎞⎛ ux ⎞⎡u ⎤+ ⎢1+ (x + uτ )erfc ⎜⎟⎥ exp ⎜⎟⎜⎝ 4D τ ⎟⎠Dx ⎦⎝ Dx ⎠⎣x−⎛ (x − uτ )2 ⎞ ⎤4u 2τ
exp ⎜ −⎟ ⎥ dτ ,π Dx
4Dxτ ⎠ ⎥
⎝
⎦

(109)

where

(±∞,0) = 0

⎡
⎛
⎞
⎛
⎞⎤
Ci1/3 ⎢y−ay+a ⎥⎟ − erfc ⎜⎟ .
erfc ⎜
2 ⎢⎜⎝ 4D y t ⎟⎠⎜⎝ 4D y t ⎟⎠ ⎥
⎢⎣
⎥⎦

⎧C i
⎪
C(x, y, z,0) = ⎨
⎪0
⎩

⎡
⎛
⎞
⎛
⎞⎤
⎢erfc ⎜ z − b ⎟ −erfc ⎜ z + b ⎟ ⎥ dτ⎜⎝ 4D τ ⎟⎠⎜⎝ 4D τ ⎟⎠ ⎥
⎢
z
z
⎣
⎦

x

(103)

can be readily adopted from the literature (e.g., Eq. [2.15] of
Crank, 1975):

C y ( y,t) =⎧⎪ 1⎡ (x − uτ )2 ⎤uCo t
exp ⎢ −⎥
∫ exp(− µτ ) ⎨4 P(t )⎢⎣ 4Dxτ ⎥⎦⎪⎩ π Dxτ
⎛ x + uτ ⎞ ⎫⎪⎛ ux ⎞
u−exp ⎜erfc ⎜⎟⎬⎟⎜⎝ 4D τ ⎟⎠ ⎪2Dx
⎝ Dx ⎠
x
⎭
⎡
⎛
⎞
⎛
⎞⎤
y−a
y+a ⎥
⎢
⎟ −erfc ⎜⎟• erfc ⎜⎢⎜⎝ 4D yτ ⎟⎠⎜⎝ 4D yτ ⎟⎠ ⎥
⎣⎢
⎦⎥

+

The solution of the transverse one-dimensional transport
equation given by Eq. (48) subject to the condition⎧ C1/3 − a < y < a⎪ i
C y ( y,0) = ⎨and⎪ 0elsewhere⎩(108)

The solution for this case is (Leij et al., 1991)⎛ x + x + ut ⎞ ⎤⎛⎞u2− ⎜1 +(x + x2 + ut)⎟ erfc ⎜
⎟⎥
⎜⎝
Dx
⎝
⎠
4Dx t ⎟⎠ ⎥⎦⎛ x − x − ut ⎞⎛ x − x − ut ⎞2
1
+ erfc ⎜⎟ − erfc ⎜⎟⎜⎝⎜⎝4Dx t ⎟⎠
4Dx t ⎟⎠

⎛ ⎡ R(x + x ) + ut ⎤ 2 ⎞ ⎤ ⎫⎦ ⎟ ⎥ ⎪.
1
− exp ⎜ − ⎣
⎬
4Dx t
⎜
⎟ ⎥⎪
⎝
⎠ ⎥⎦ ⎭

(107)

(106)

⎪⎧ 0 0 < t ≤ toP(t) = ⎨⎪⎩ t − to t > to .

(110)

Here we illustrate a two-dimensional application of Eq.
(109) involving the injection of a solute over a finite horizontal
cross-section (–10 ≤ y < 10 m), 0 ≤ z < ∞) of an initially solute
free river channel or other open surface water body. As before
we use a mean transport velocity of 1 m s-1, and longitudinal(Dx) and lateral (Dy) dispersion coefficients of 5 and 0.1 m2 s-1,
respectively. Production or decay processes are neglected (µ = γ= 0). The two-dimensional (x, y) concentration distribution attime 1500 s is shown in Fig. 4. Calculations were carried out
using the 3DADE code (Leij and Bradford, 1994) within
STANMOD assuming a = 10 m and b = 1000 m (the latterlarge value ensuring a two-dimensional transport geometry).
Fig. 5 further shows for this example the concentration distribution along the longitudinal (x) coordinate in the middle ofthe channel (i.e., for y = z = 0). The width (b) of the solutesource along the z-coordinate perpendicular to the river surface(Eq. (108)) was taken to be very large in this example (1000m), thus creating a two-dimensional (x, y) transport geometry.
Smaller values of the source width, b (e.g., 5 or 2 m) would

157

Martinus Th. van Genuchten, Feike J. Leij, Todd H. Skaggs, Nobuo Toride, Scott A. Bradford, Elizabeth M. Pontedeiro

create a three-dimensional geometry and significantly lower theintermediate concentrations in Fig. 4 (between approximately
200 and 1300 m) because of lateral dispersion away from thecenter line. The same effect would occur if higher values of thelateral dispersion coefficient (Dy) were used.

Fig. 4. Two-dimensional (x, y) concentration distributions at thesurface of a stream upon continuous injection of a solute having a
concentration of 10 g/m3 into a cross-section –10 ≤ y ≤ 10 (m), –1000 ≤ z < 1000 m as obtained with Eq. (109) using the 3DADE
module within STANMOD (a = 10 m, b = 1000 m, u = 1 m/s, Dx =5 m2/s, Dy = Dz = 0.1 cm2/s, µ = γ = 0).

Fig. 5. Concentration distribution along the longitudinal (x) coordinate in the middle of the stream (y = z = 0) for the problem depicted in Fig. 4.

The above example is just one possible scenario for predicting two- and three-dimensional transport in streams and larger
surface water bodies. We refer to Leij and Dane (1990), Leij etal. (1991), Leij and Bradford (1994) and Weerts et al. (1995)
for many other solutions, including for nonequilibriumtransport as discussed in Part 2. Many of the solutions are included in the 3DADE (Leij and Bradford, 1994) and N3DADE(Leij and Toride, 1997) modules within the STANMOD software package, which we used for all calculations in this study.
Fig. 6 gives an overview of the various geometries that can beconsidered with this software for multidimensional transport
problems.

158

CONCLUDING REMARKSThe solutions in this paper pertain to one-dimensional longitudinal equilibrium transport in streams and rivers, and longitudinal transport and lateral dispersion in rivers and larger surface
water bodies. This ideal situation generally does not occur instreams and rivers because of the presence of relatively immobile or stagnant (dead) zones of water connected to the mean
stream. Such stagnant zones include pools and eddies along theriver banks, water isolated behind rocks, gravel or vegetation,or relatively inaccessible water along an uneven river bottom.
Fluvial sediments and more generally the subsurface hyporheic
zone may also contribute to the presence of such relatively
stagnant water zones. By providing sinks or sources of soluteduring transport, stagnant water zones typically cause tailing inobserved concentration distributions, which cannot be predictedwith the conventional equilibrium ADE formulations. Part 2 ofthis two-part series (van Genuchten et al., 2013) will focus on
nonequilibrium transport caused by the presence of such stagnant water zones (transient storage models), and also includes
models that consider simultaneous longitudinal advectivedispersive transport in a river and lateral diffusion into and outof the hyporheic zone. Part 2 additionally provides several
solutions for the transport of solutes involved in consecutive
decay chain reactions.

Fig. 6. Available geometries, initial conditions, and inlet boundary
conditions (BC’s) that can be simulated using the 3DADE module
within the STANMOD software package. Both transient (timedependent) and steady-state (time-independent) problems can beconsidered.

Most of the analytical solutions presented in this part 1 were
derived from solutions to mathematically very similar problems
in subsurface contaminant transport. All part 1 solutions havebeen incorporated in the public-domain windows-based
STANMOD software package (Šimůnek et al., 2000). This
software package also includes parameter estimation capabilities (e.g., Fig. 2), and hence provides a convenient tool for
analyzing observed contaminant concentration distributions
versus distance and/or time.
While inherently less flexible than more comprehensive numerical models for contaminant transport in streams and rivers,
we believe that the analytical solutions assembled in this study
can be very useful for a variety of applications, such as forproviding initial or approximate analyses of alternative pollution scenarios, conducting sensitivity analyses to investigate theeffects of various parameters or processes on contaminantExact analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equationtransport in rivers and streams, extrapolating results over large
times and spatial scales where numerical solutions become
impractical, serving as screening models, providing benchmark
solutions for more complex transport processes that cannot be
solved analytically, and for validating more comprehensive
numerical solutions of the governing transport equations.
REFERENCES

Anderson, E.J., Phanikumar, M.S., 2011. Surface storage dynamics in large rivers: Comparing three-dimensional particle
transport, one-dimensional fractional derivative, and multirate transient storage models. Water Resour. Res., 47,
W09511, doi: 10.1029/2010WR010228.
Batu, V., van Genuchten, M.Th., 1990. First- and third-typeboundary conditions in two-dimensional solute transportmodeling. Water Resour. Res., 26(2), 339–350.
Batu, V., van Genuchten, M.Th., Parker, J.C., 2013. Authors
reply. Ground Water, 51, 1, 1–9.
Bear, J., Verruijt, A., 1987. Modeling Groundwater Flow and
Pollution. Kluwer Acad. Publ., Norwell, MA.
Bencala, K.E., 1983. Simulation of solute transport in a mountain pool-and-riffle stream with a kinetic mass transfer model for sorption. Water Resour. Res., 19, 732–738.
Bencala, K.E., 2005. Hyporheic exchange flows. In: Anderson,
M.P. (Ed.): Encyclopedia of Hydrology Science, pp. 1733––1740, John Wiley, New York.
Bencala, K., Walters, R., 1983. Simulation of transport in a
mountain pool-and-riffle stream: A transient storage model.
Water Resour. Res., 19, 718–724.
Bencala, K.E., Jackman, A.P., Kennedy, V.C., Avanzino, R.J.,
Zellweger, G.W., 1983. Kinetic analysis of strontium andpotassium sorption onto sands and gravels in a natural channel. Water Resour. Res., 19, 725–731.
Brenner, H., 1962. The diffusion model of longitudinal mixing
in beds of finite length. Numerical values. Chem. Eng. Sci.,
17, 229–243.
Bruch, J.C., Street. R.L., 1967. Two-dimensional dispersion. J.
Sanit. Eng. Div. Proc. Am. Soc. Civ. Eng., 93(A6), 17–39.
Carslaw, H.S., Jaeger, J.C., 1959. Conduction of Heat in Solids.
Clarendon Press, Oxford.
Cleary, R.W., 1973. Unsteady-state, multi-dimensional modeling of water quality in rivers. Water Resour. Program, Tech.
Rep. 1-73-1, Dept. of Civil Engineering, Princeton Univ.,
Princeton, NJ.
Cleary, R.W., Adrian, D.D., 1973. Analytical solution of theconvective-dispersive equation for cation adsorption in soils.
Soil Sci. Soc. Am. Proc., 37, 197–199.
Crank, J., 1975. The Mathematics of Diffusion. Clarendon
Press, Oxford.
De Smedt, F., Brevis, W., Debels, P., 2005. Analytical solution
for solute transport resulting from instantaneous injection in
streams with transient storage. J. Hydrol., 315, 25–39.
De Smedt, F., Brevis, W., Debels, P., 2006. Reply to commenton ‘Analytical solution for solute transport resulting from instantaneous injection in streams with transient storage’. J.
Hydrol., 330, 761–762.
Field, M.S., 1997. Risk assessment methodology for karst aquifers: (2) Solute transport modeling. Env. Modeling and Assessment, 47, 23–37.
Fischer, H.B., List, E., Koh, R.C.Y., Imberger, J., Brooks,
N.H., 1979. Mixing in Inland and Coastal Waters. Academic
Press, Orlando, FL.
Gautschi, W. 1964. Error function and Fresnel integrals. p.
295–329. In: Abramowitz M., Stegun, I.A.: Handbook ofMathematical Functions. Appl. Math. Ser. 55. National Bureau of Standards, Washington, D.C.
Gerecht, K.E., Cardenas, M.B., Guswa, A.J., Sawyer, A.H.,
Nowinski, J.D., Swanson, T.E., 2011. Dynamics of hyporheic flow and heat transport across a bed-to-bank continuum ina large regulated river. Water Resour. Res., 47, W03524,
doi: 10.1029/2010WR009794.
Huang, J., Goltz M.N., Roberts, P.V., 2006. Comment on ‘Analytical solution for solute transport resulting from instantaneous injection in streams with transient storage’ by F. De
Smedt, W. Brevis and D. Debels. J. Hydrol., 330, 759–760.
Jury, W.A., Horton, R.H., 2004. Soil Physics. 6th Ed. John
Wiley & Sons, Inc., Hoboken, NJ, 370 pp.
Kreft, A., Zuber, A., 1978. On the physical meaning of thedispersion equation and its solutions for different initial andboundary conditions. Chem. Eng. Sci., 33, 1471–1480.
LeGrand-Marcq, C., Laudelout, H., 1985. Longitudinal dispersion in a forest stream. J. Hydrol., 78, 317–324.
Leij, F.J., Bradford, S.A., 1994. 3DADE: A computer programfor evaluating three-dimensional equilibrium solute transportin porous media. Research Report 134, U.S. Salinity Laboratory, USDA, ARS, Riverside CA, 81 pp.
Leij, F.J., Dane, J.H., 1990. Analytical solution of the onedimensional advection equation and two- or threedimensional dispersion equation. Water Resour. Res., 26,1475–1482.
Leij, F.J., Toride, N., 1995. Discrete time- and length-averaged
solutions of the advection-dispersion equation. Water Resour. Res., 31, 1713–1724.
Leij, F.J., Toride, N., 1997. N3DADE: A computer program forevaluating
three-dimensional
nonequilibrium
solutetransport in porous media. Research Report 143, U.S. Salinity Laboratory, USDA, ARS, Riverside CA. 116 pp.
Leij, F.J., van Genuchten, M.Th., 2000. Analytical modeling of
nonaqueous phase liquid dissolution with Green’s functions.
Transport in Porous Media, 38, 41–166.
Leij, F.J., Skaggs, T.H., van Genuchten, M.Th., 1991. Analytical solutions for solute transport in three-dimensional semiinfinite porous media. Water Resour. Res., 27, 2719–2733.
Leij, F.J., Toride, N., van Genuchten, M.Th., 1993. Analyticalsolutions for non-equilibrium solute transport in threedimensional porous media. J. Hydrol., 151, 193–228.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear parameters. SIAM J., 11, 431–441.
O’Connor, B.L., Hondzo, M., Harvey, J.W., 2009. Predictive
Modeling of Transient Storage and Nutrient Uptake: Implications for Stream Restoration. Journal of Hydraulic Engineering, 136, 1018–1032.
Onishi, Y., Serne, R.J., Arnold, E.M., Cowan, C.E., Thompson,
F.L., 1981. Critical review: Radionuclide transport, sedimenttransport, and water quality mathematical modeling; and radionuclide adsorption/desorption mechanisms. NUREG/CR1322, PNL-2901, Pacific Northwest Laboratory, Richland,WA.
Özişik, M.N., 1989. Boundary Value Problems of Heat Conduction. Dover, New York.
Parker, J.C., van Genuchten, M.Th., 1984a. Flux-averaged
concentrations in continuum approaches to solute transport.
Water Resour. Res., 20, 866–872.
Parker, J.C., van Genuchten, M.Th., 1984b. Determining
transport parameters from laboratory and field tracer experiments. Virginia Agric. Exp. Stat., Bull., 84–3, Blacksburg,
VA.
Pérez Guerrero, J.S., Pimentel, L.C.G., Skaggs, T.H., vanGenuchten, M.Th., 2009. Analytical solution of the advec-

159

Martinus Th. van Genuchten, Feike J. Leij, Todd H. Skaggs, Nobuo Toride, Scott A. Bradford, Elizabeth M. Pontedeiro

tion-diffusion transport equation using a change-of-variable
and integral transform technique. Int. J. Heat Mass Transf.,
52, 3297–3304.
Pérez Guerrero, J.S., Pimentel, L.C.G., Skaggs, T.H., vanGenuchten M.Th., 2009. Analytical solution of the advection-diffusion transport equation using a change-of-variable
and integral transform technique. Int. J. Heat Mass Transf.,
52, 3297–3304.
Pérez Guerrero, J.S., Pontedeiro, E.M., van Genuchten, M.Th.,
Skaggs, T.H. 2013. Analytical solutions of the onedimensional advection-dispersion solute transport equation
subject to time-dependent boundary conditions. Chem. Eng.
J., doi: http://dx.doi.org/10.1016/j.cej.2013.01.095.
Runkel, R.L. 1998. One-dimensional transport with inflow andstorage (OTIS): A solute transport model for streams andrivers. USGS Water Resources Investigations Report 984018, U.S. Geological Survey, Denver, CO. Available fromhttp://co.water.usgs.gov/otis/.
Runkel, R.L., Chapra, S.C., 1993. An efficient numerical solution of the transient storage equations for solute transport insmall streams. Water Resour. Res., 29, 11–215.
Runkel, R.R., McKnight, D., Rajaram, H. (Eds.), 2003. Modeling Hyporheic Processes. Special issue, Adv. Water Resour.,26, 9, 901–1040.
Sagar, B., 1982. Dispersion in three dimensions: approximate
analytic solutions. J. Hydraul. Div. Proc. Am. Soc. Civ.
Eng., 108(HY1), 47–62.
Sayre, W.W., 1973. Natural mixing processes in rivers. In:
H.W. Shen (ed.): Environmental Impact on Rivers, University of Colorado, Fort Collins, CO.
Selim, H.M., Mansell, R.S., 1976. Analytical solution of theequation for transport of reactive solutes through soils. Water Resour. Res., 12, 528–532.
Šimůnek, J., van Genuchten, M.Th., Sejna M., Toride, N., Leij,
F.J., 2000. The STANMOD computer software for evaluating solute transport in porous media using analytical solutions of convection-dispersion equation. Version 2.0, IGWMC-TPS-71, Int. Ground Water Modeling Center (IGWMC), Colorado School of Mines, Golden, Colorado, 32
pp. Available from http://www.pc-progress.com/en/Default.
aspx?stanmod.
Sneddon, I.N., 1995. Fourier Transforms. Dover Publications,
New York.
Spiegel, M.R., 1965. Laplace Transforms. McGraw-Hill, NewYork.
Thackston, E.L., Schnelle, K.B. Jr., 1970. Predicting the effects
of dead zones on stream mixing. J. Sanit. Eng. Div., Am.
Soc. Civil Eng., 96(SA2), 319–331.
Thomann, R.V., 1973. Effect of longitudinal dispersion on
dynamic water quality response of streams and rivers. Water
Resour. Res., 9(2), 355–366.

160

Toride, N., Leij, F.J., van Genuchten, M.Th., 1993a. A comprehensive set of analytical solutions for nonequilibrium solutetransport with first-order decay and zero-order production.
Water Resour. Res., 29(7), 2167–2182.
Toride, N., Leij, F.J., van Genuchten, M.Th., 1993b. Fluxaveraged concentrations for transport in soils having nonuniform initial solute concentrations. Soil Sci. Soc. Am., 57(6),1406–1409.
Toride, N., Leij, F.J., van Genuchten, M.Th., 1999. The CXTFIT code for estimating transport parameters from laboratory
or field tracer experiments. Version 2.1. Research Report
No. 137, U.S. Salinity Laboratory, USDA, ARS, Riverside,
CA. 121 pp.
Van Genuchten, M.Th., 1981. Analytical solutions for chemicaltransport with simultaneous adsorption, zero-order production and first-order decay. J. Hydrol., 49, 213–233.
Van Genuchten, M.Th., 1985. Convective-dispersive transportof solutes involved in sequential first-order decay reactions.
Computers & Geosciences, 11(2), 129–147.
Van Genuchten, M. Th., Alves, W.J., 1982. Analytical Solutions of the One-Dimensional Convective-Dispersive Solute
Transport Equation. Techn. Bull., 1661, U.S. Department ofAgriculture, Washington, DC.
Van Genuchten, M.Th., Parker, J.C., 1984. Boundary conditions for displacement experiments through short laboratory
columns. Soil Sci. Soc. Am. J., 48, 703–708.
Van Genuchten, M.Th., Gorelick, S.M., Yeh, W.-G. 1988.
Application of parameter estimation techniques to solutetransport studies. In: DeCoursey, D.G. (Ed.): Proc. Int.
Symp. Water Quality Modeling of Agricultural Non-Point
Sources, Part 2, June 19–23, Utah State University, Logan,UT.
Van Genuchten, M.Th., Leij, F.J., Skaggs, T.H., Toride, N.,
Bradford, S.A., Pontedeiro, E.M., 2013. Exact analytical solutions for contaminant transport in rivers. 2. Transient storage and decay chain solutions. J. Hydrol. Hydromech.,61(3), (in press).
Yeh, G.T., Tsai, Y.-J., 1976. Analytical three-dimensionaltransient modeling of effluent discharges. Water Resour.
Res., 12(3), 533–540.
Weerts, A.H., Leij, F.J., Toride, N., van Genuchten, M.Th.,
1995. Analytical models for contaminant transport in soils
and groundwater. Unpublished Report, U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, 337 pp.
Wörman, A. 1998. Analytical solution and timescale fortransport of reacting solutes in rivers and streams. Water Resour. Res., 34(10), 2703–2716.
Zwillinger, D., 1989. Handbook of Differential Equations.
Academic Press, San Diego, CA.
Received 8 August 2012
Accepted 23 February 2013