J. Hydrol. Hydromech., 61, 2013, 3, 250–259
DOI: 10.2478/johh-2013-0032

Exact analytical solutions for contaminant transport in rivers
2. Transient storage and decay chain solutions
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: Contaminant transport processes in streams, rivers, and other surface water bodies can be analyzed or predicted using the advection-dispersion equation and related transport models. In part 1 of this two-part series we presented alarge number of one- and multi-dimensional analytical solutions of the standard equilibrium advection-dispersion equation (ADE) with and without terms accounting for zero-order production and first-order decay. The solutions are extended in the current part 2 to advective-dispersive transport with simultaneous first-order mass exchange between the stream
or river and zones with dead water (transient storage models), and to problems involving longitudinal advectivedispersive transport with simultaneous diffusion in fluvial sediments or near-stream subsurface regions comprising ahyporheic zone. Part 2 also provides solutions for one-dimensional advective-dispersive transport of contaminants subject to consecutive decay chain reactions.
Keywords: Contaminant transport; Analytical solutions; Surface water; Transient storage models; Solute decay chains.
INTRODUCTION
Contaminant transport processes in streams and rivers havebeen analyzed and predicted now for several decades using a
range of mathematical models (Bencala, 1983; DeSmedt et al.,2007; Fisher et al., 1979; Runkel et al., 2003; Thomann, 1973).
Because many complex and often nonlinear physical, chemicaland biological processes affect contaminant transport in streams
and rivers, transport model formulations now increasingly
utilize partial differential equations which must be solved using
numerical methods. Exact analytical solutions are generally
available only for more simplified formulations based on idealized representations of the transport domain and transport processes. Yet despite the apparent simplifications, analyticalsolute transport models remain useful for many applications
such as providing initial or approximate analyses of a variety ofcontaminant transport scenarios, especially for longer times,
facilitating transport analyses when insufficient data are available to warrant the use of a comprehensive numerical model,
and for serving as a benchmark to test the accuracy of numerical models.
In this two-part series, we assembled a large number of analytical solutions for modeling solute transport in streams, rivers
and other open surface water bodies. In part 1 (van Genuchten
et al., 2013) we summarized solutions for one- and multidimensional equilibrium transport, with and without zero-orderproduction and first-order decay. In the current part 2 we provide solutions for transport with simultaneous first-order exchange between the river and relatively immobile or stagnant
water zones (transient storage models) and for situations where
exchange with the hyporheic zone is modeled as a diffusion
process (further referred to here as hyporheic zone diffusion
models). Also presented in this paper are several solutions fortransport of solutes subject to consecutive decay chain reactions.

250

Most of the presented solutions were derived from solutions
to mathematically very similar problems in subsurface contaminant transport (Leij and Toride, 1997; Leij et al., 1993; Toride
et al., 1993; van Genuchten, 1985a,b; Weerts et al., 1995).
Except for those pertaining to diffusion in hyporheic zones, allsolutions have been incorporated in the public-domain windows-based STANMOD software package (Šimůnek et al.,2000).
TRANSIENT STORAGE MODELSAll solutions presented in Part 1 (van Genuchten et al.,2013) hold for equilibrium contaminant transport characterized
by relatively symmetrical or sigmoidal concentration distributions versus time or distance, unless modified by production
and degradation processes, or special initial or boundary conditions. This ideal situation generally does not occur in streams
and rivers because of the presence of relatively immobile or
stagnant zones of water connected to the mean stream channel.
Such stagnant zones include pools and eddies along the riverbanks, water isolated behind rocks, gravel or vegetation, or
relatively inaccessible water along an uneven river bottom.
Fluvial sediments and more generally subsurface hyporheic
zones may also contribute to the presence of such relatively
stagnant water. By providing sinks or sources of solute duringtransport, stagnant water zones typically cause tailing in observed concentration distributions, which cannot be predictedwith the conventional equilibrium ADE formulation.
Several conceptual approaches are possible for modeling solute exchange between the river and stagnant water zones. One
popular approach inherent in most transient storage models is to
assume first-order mass transfer between the river and stagnant
water (Bencala and Walters, 1983; De Smedt, 2006; LeGrandMarcq and Laudelout, 1985; Runkel et al., 1996; Thackston andSchnelle, 1970, 2000; Wlostowski et al., 2013, and many others). The resulting model assumes that no advective and disper-Exact analytical solutions for contaminant transport in rivers: 2. Transient storage and decay chain solutions

sive transport occurs in the transient storage zone, and that an
average concentration can be assigned to this zone at any valueof time, t, and longitudinal distance, x, along the stream. A
more refined approach would be to assume that exchange ofcontaminant with stagnant zones, or the hyporheic zone ingeneral, occurs by diffusion (De Smedt, 2007; Grant et al.,2012; Jackman et al., 1984; Wörman, 1998). Simulating thissituation requires coupling the transport equation for the stream
or river with some diffusion model for the near-stream subsurface environment (De Smedt, 2007; Jackman et al., 1984). In
this paper we first consider the relatively traditional quasiempirical approach based on first-order mass transfer (transientstorage models), while subsequently we consider more elaborate models assuming diffusion-based formulations for exchange between the river and the hyporheic zone (referred tohere as hyporheic diffusion models).
Assuming a constant cross-sectional area A of the river and aconstant longitudinal dispersion coefficient Dx, the transientstorage model in its most general form is given by (e.g., Bencala et al., 1983; Runkel and Chapra, 1993; Wlostowski et al.,
2013)

∂C
∂2 C
∂C QL
= Dx
−u
+
(C L − C) − α (C − Cs ) − µC ,
2
∂t
∂x
A
∂x

(1)

∂Cs
A
=α
(C − Cs ) − µ s Cs ,
∂t
As

(2)

tion and first-order degradation are given by Toride et al.
(1993). Two solutions are listed here as an example.
Case E1. We give here first the solution of Eqs. (3) and (4)for the initial condition

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

(5)the third-type inlet condition⎛
Dx ∂C ⎞⎜⎝ C − u ∂x ⎟⎠

⎧C o
⎪
=⎨
x=0+ ⎪ 0⎩

0 < t ≤ to(6)t ≥ toand a semi-infinite profile∂C(∞, t) = 0 .
∂x(7)

The solutions for the stream and the dead storage zone concentrations are, respectively,⎧ 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

(8)

and

where C and Cs are concentrations of the stream and storage
zones, respectively (ML–3), Dx is the longitudinal dispersion
coefficient (L2T–1), u is the longitudinal fluid flow velocity(LT–1), x is the longitudinal coordinate (L), t is time (T), As isthe cross-sectional area of the storage zone (L2), QL is the lateral volumetric inflow rate per unit length (L3T–1L–1), CL is theconcentration of the lateral inflow (ML–3), α is a dead-zone
storage mass transfer coefficient (T–1), and µ and µs are firstorder decay coefficients for the stream and storage zone, respectively (T–1). The terms on the right-hand side denote longitudinal dispersive and advective transport, lateral inflow fromgroundwater, and solute exchange with the storage zone. Several variations of this model have been applied to river transport;
they typically require a numerical solution, particularly when A,
As and QL vary with distance. However, as pointed out by van
Genuchten et al. (1988) and Huang et al. (2006), analyticalsolutions are readily available for mathematically very similar
problems of nonequilibrium transport in porous media (Coatsand Smith, 1964; Leij and Toride, 1997; Toride et al., 1993;
van Genuchten and Wagenet, 1989; van Genuchten and
Wierenga, 1976).
If the lateral inflow and degradation terms in (1) and (2) are
negligible, then the dead zone storage model for constant A, As,
Dx and α reduces to⎧ Ci + (Co − Ci ) B(x,t)
0 < t ≤ to⎪
Cs (x,t) = ⎨⎪⎩ Ci + (Co − Ci ) B(x,t) − Co B(x,t − to ) t > to ,
where

⎧ u2
t
⎡ (x − uτ )2 ⎤⎪A(x,t)= ∫ J (a,b) ⎨exp ⎢ −⎥π Dxτ4Dxτ ⎥⎦
0⎢⎣⎩⎪⎡ x + uτ ⎤ ⎫⎪⎛ ux ⎞u2⎢⎥ ⎬ dτ ,
−
exp ⎜
erfc
⎟
Dx
⎢⎣ 4Dxτ ⎥⎦ ⎪⎝ Dx ⎠⎭⎧ u2t⎡ (x − uτ )2 ⎤⎪B(x,t)= ∫ ⎡⎣1− J (b,a) ⎤⎦ ⎨exp ⎢ −⎥
0⎢⎣ 4Dxτ ⎥⎦⎪⎩ π Dxτ
⎡ x + uτ ⎤ ⎫⎪⎛ ux ⎞u2⎢⎥ ⎬ dτ ,
−
exp ⎜
erfc
⎟
Dx
⎢⎣ 4Dxτ ⎥⎦ ⎪⎝ Dx ⎠
⎭

a =αt

2

∂C
∂ C
∂C
= Dx
−u
− α (C − Cs ) ,
∂t
∂x
∂x 2

(3)

∂C s = α A (C − C )
s .
∂t
As

(4)

(9)

b = α (t – τ )

A
As

.

(10)

(11)

(12)

These expressions contain Goldstein's J-function which isdefined as (Goldstein, 1953):

Many analytical solutions for one-dimensional nonequilibrium transport with or without accounting for zero-order produc-a

J (a,b) = 1− exp(−b) ∫ exp(−ξ ) I0(2 bξ ) dξ

(13)

0

in which I0 is the zero-order modified Bessel function.

251

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

The above solutions are for a third-type inlet boundary condition. For a first-type inlet condition, (10) and (11) must be
replaced by, respectively,⎧t⎡ (x − uτ )2 ⎤ ⎫⎪x2⎪
A(x,t)= ∫ J (a,b) ⎨exp⎢−⎥ ⎬dτ ,
30⎢⎣ 4Dxτ ⎥⎦ ⎪⎭⎪⎩ 4π Dτ(14)⎧t
⎡ (x − uτ )2 ⎤ ⎫⎪x2⎪
B(x,t)= ∫ ⎡⎣1− J (b,a) ⎤⎦ ⎨exp ⎢ −⎥ ⎬dτ .
0⎢⎣ 4Dxτ ⎥⎦ ⎪⎭⎪⎩ 4π Dxτ(15)

v = βuD = β Dxβ=A
A + Asω=αβ L, (16)uwhere L is some characteristic length used to place the transportmodel in dimensionless form (L = 1000 m in the current example). For the calculations of Fig. 1 we hence used in CXTFIT
the parameter values v = 0.83333 m/s, D = 4.16667 m/s2, β =0.83333, and ω = 0, 0.8333 and 8.333 1/s).
We also used the above transient storage model to analyze
the experimental data (exp. I–3) of Brevis et al. (2001) and De
Smedt et al. (2005) that were examined earlier in part 1 using
the equilibrium transport model. The data are for a tracer experiment conducted in the Chillán River in Chile in which 157.1 gof a 20% Rhodamine WT tracer was injected at location x = 0,
with measurements made at L = 4604 m downstream of theinjection point. Fig. 2 shows the data along with fitted curves
obtained with both the transient storage model (solid line) andthe standard ADE equilibrium transport model discussed in part1 (dashed line). Parameters were estimated using the nonlinear
least-squares optimization features of CXTFIT, which provided
estimates of the four parameters given by Eq. (16), as well as
the concentration Co of the applied tracer solution for a given
value of the injection or pulse time to in Eq. (8).
7Transient Storage Model
Equilibrium Model

63
Concentration (mg/m )Fig. 1 shows an application of the above solution for a thirdtype boundary condition (Eq. (6)). Like for several examples
presented in part 1 (van Genuchten et al., 2013), the problem
draws upon calculations and parameter values used by De
Smedt et al. (2005). The example involves the injection of 1 kgof a solute in the main channel of a river having a crosssectional area (A) of 10 m2, a connected storage zone (As) of 2
m2, an average flow velocity of 1 m/s in the river, and a dispersion coefficient of 5 m2/s. Fig. 1 shows for three values of thedead-storage zone mass transfer coefficient (α) calculated solute concentrations at a distance 1000 m downstream from theinjection point. The curves were obtained assuming that 1 kgof solute was injected during a time period of only 10 seconds(to = 10 in Eq. (6)). De Smedt et al. (2005) used for this purposea Dirac function, similarly as we selected in part 1 (Fig. 1). The
resulting curves however are essentially identical, i.e., our Fig.
1 here and Fig. 1 of De Smedt et al. (2005). The curve for α = 0
assumes no exchange of solute with the storage zone, and hence
could be calculated also immediately with the equilibriumtransport model.

rium ADE option in CXTFIT. In that case the average porewater velocity v, the dispersion coefficient D, and the dimensionless parameters β and ω in CXTFIT must be defined interms of transient storage model parameters as follows

5
4
3
2
1
0
5000

7500

10000

12500

15000

17500

20000

Time (s)Fig. 2. Observed (solid squares) and fitted (continuous and dashedlines) concentrations for tracer experiment I-3 of Brevis et al.
(2001) and De Smedt et al. (2005).
Fig. 1. Calculated concentration distributions obtained with Eq. (8)for near-instantaneous injection (to = 10 s) of 1 kg of solute (Co =10 g/m3) in the main channel of a stream having a cross-sectional
area of 10 m2, an average flow velocity of 1 m/s, a dispersion
coefficient of 5 m2/d, and a connected transient storage zone of 2
m2, assuming three values of the mass transfer coefficient, α.

Calculations for Fig. 1 were obtained with the CXTFIT codeof Toride et al. (1999) as incorporated in the STANMOD software of Šimůnek et al. (2000). Even though CXTFIT was derived for porous media transport problems, the code is applicable immediately to most or all models listed in this paper, as
well as in part 1, by assuming a volumetric water content of 1.0
in the code. However, some care is needed when the transientstorage model is simulated using the deterministic nonequilib-

252

Assuming a very short injection time of only 10 s, similarly
as in part 1 for the equilibrium analysis, we obtained the following parameter values (with their 95% confidence intervals):v = 0.414 ± 0.002 m/s, D = 2.66 ± 1.40 m2/s, β = 0.8647 ±
0.0304, ω = 2.934 ± 0.936, and Co = 1808 ± 33 mg/m3. The
coefficient of determination (R2) of the fit was 0.997, and theroot mean square error (RMSE) 0.299 mg/m3. Using Eq. (16),the CXTFIT estimates translate to the following parameters in
the transient storage model: u = 0.478 m/s, Dx = 3.07 m2/s, As/A
= 0.156, and α = 3.52 x 01–4 s–1. The total amount of solute
mass (m) injected per m2 cross-sectional can now be calculatedusing m = u Co to, or 8.642 g/m2. Given that a total amount of
157.1 g of tracer was applied to the river, this translates to aneffective cross-sectional area of 157.1 (g)/8.642 (g/m2) or 18.2

Exact analytical solutions for contaminant transport in rivers: 2. Transient storage and decay chain solutions

m2 for the river channel as seen by the transient storage model.
Our estimates for u, A, As are exactly the same as those obtained by De Smedt et al. (2005), while estimates of the dispersion coefficient Dx, and the mass transfer coefficient differedslightly. De Smedt et al. (2005) obtained values of 3.29 ± 1.05
m2/s and 2.87 x 10–4 ± 0.89 x 10-4 s–1 for these two parameters,
respectively.
The estimates of the transport parameters above were obtained assuming a solute injection pulse of 10 s. Similarly as forthe equilibrium analysis in part 1, the same results were obtained when the total mass was assumed to be applied instantaneously as modeled using the Dirac solution of Eqs. (1) and (2)(not further given here) as used also by De Smedt et al. (2005),
and for pulse times to up to about 100 s. We also obtained essentially the same results when assuming a first-type inlet
boundary condition, which again shows that differences between the first- and third-type solutions generally are verysmall at locations having relatively large values of the dimensionless distance variable ux/Dx.
Case E2. We also give the solution here of the general transient storage model given by Eqs. (1) and (2) for degradation inboth the stream and storage zone, as well as lateral inflow (QL
is assumed to be constant). Analytical solutions are given for an
initially solute-free river system (Ci = 0) subject to boundary
conditions (6) and (7). The solutions for the stream and storagezone concentrations are, respectively (Toride et al., 1993),⎧Co A1(x,t) + B1(x,t)
0 < t ≤ to⎪C(x,t) = ⎨⎪⎩Co A1(x,t) + B1(x,t) − Co A1(t − to ) t > to

(18)

(19)

(20)

⎡ αµsτ
⎤Q
i exp ⎢ −− ( µ + L )τ ⎥ E(x,τ )dτ ,A ⎦⎣ α + µs
B2 (x,t)=α QLC L t
∫ ⎡1− J (b,a) ⎤⎦(α + µs ) A 0 ⎣⎡ αµsτ
⎤Q
i exp ⎢ −− ( µ + L )τ ⎥ F(x,τ )dτA ⎦⎣ α + µs

in whicha=

α 2τ
α + µsb=(α + µ s )(t − τ ) A
.
As

(25)

The above solution is for a third-type inlet condition. The solution for a first-type condition is exactly the same, except forEqs. (23) and (24) which must be replaced by, respectively,
2⎡ (x − uτ )2 ⎤x
exp ⎢ −⎥
3
4Dxτ ⎥⎦
4π Dxτ
⎢⎣

(26)

and

⎡ αµsτ
⎤QLC L tQ
− ( µ + L )τ ⎥ F(x,τ ) dτ ,∫ J (a,b) exp ⎢ −A 0
A ⎦⎣ α + µs

α t
∫ ⎡1− J (b,a) ⎤⎦α + µs 0 ⎣

(24)

(17)

where
t
⎡ αµsτ
⎤Q
A1(x,t)= ∫ J (a,b)exp ⎢ −− ( µ + L )τ ⎥ E(x,τ ) dτ ,A ⎦0
⎣ α + µs

(23)

and

F(x,τ ) = 1−

⎧Co A2 (x,t) + B2 (x,t)
0 < t ≤ to⎪
C s(x, t) = ⎨⎪⎩Co A2 (x,t) + B2 (x,t) − Co A2 (t − to ) t > to ,A2 (x,t)=2
2
u exp ⎡⎢ − (x − uτ ) ⎤⎥π Dxτ4Dxτ ⎥⎢⎣⎦
2⎡⎛ ux ⎞
x + uτ ⎤⎥,− u exp ⎜erfc ⎢⎟⎝ Dx ⎠⎢⎣ 4Dxτ ⎥⎦Dx⎡ x − ut ⎤⎡ (x − ut)2 ⎤
1
u 2t
⎥−
F(x,t) = 1− erfc ⎢exp ⎢ −⎥
2
π Dx
4Dx t ⎥⎦⎢⎣ 4Dx t ⎥⎦⎢⎣⎡ x + ut ⎤⎛ ux ⎞
1⎛ux u 2t ⎞⎥+ ⎜ 1++
erfc ⎢⎟ exp ⎜⎟
2⎝
Dx Dx ⎠⎢⎣ 4Dx t ⎥⎦⎝ Dx ⎠E(x,τ ) =and

B1(x,t)=E(x,τ ) =

(21)

(22)

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

Additional one-dimensional analytical solutions for the transient storage model are readily derived from the solutions given
by Toride et al. (1993a) for transport in porous media. They are
all incorporated in the CXTFIT code (Toride et al., 1999) aspart of STANMOD.
Multidimensional versions of the transient storage nonequilibrium model can be obtained by adding terms for the transverse dispersive flux similarly as for the equilibrium case (seeEq. (41) of van Genuchten et al., 2012). Because the nonequilibrium effects generally manifest themselves only in the longitudinal flow direction, solutions for multidimensional transientstorage models can often be deduced from available onedimensional nonequilibrium solutions. We refer to Leij et al.
(1993) for a complete set of analytical solutions for threedimensional nonequilibrium transport, again for porous media
transport. Green’s functions have proved to be particularly
convenient for constructing solutions for nonequilibriumtransport in multiple dimensions and for different mathematicalconditions (cf. Leij and van Genuchten, 2000). Many or most of
the multi-dimensional solutions are included in the 3DADE(Leij and Bradford, 1994) and N3DADE (Leij and Toride,1997) computer programs for equilibrium and nonequilibriumtransport, respectively. Like CXTFIT for one-dimensionaltransport, 3DADE (but not N3DADE) includes parameter estimation capabilities to estimate selected transport parametersfrom observed contaminant concentration distributions versus
distance and/or time.

253

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

HYPORHEIC ZONE DIFFUSION MODELSThe transient storage models in the previous section use afirst-order mass transfer equation to account for solute exchange between the river and relatively stagnant zones along
the river banks or bottom. This conceptual picture can be refined by using Fick's law to describe solute diffusion from thestream into the stagnant water zones, the fluvial sediment ormore generally the entire hyporheic zone (Jackman et al., 1985;
Runkel et al., 2003; Wörman, 1998). Similar problems havebeen described and solved analytically for contaminanttransport in fractured or macroporous media (Sudicky andFrind, 1982; Tang et al., 1981; van Genuchten, 1985a; van
Genuchten et al., 1984). Typically, the porous media solutions
account for advective-dispersive transport through well-defined
fractures or soil macropores with simultaneous diffusion fromthe fractures into the surrounding soil matrix. Many of thesesolutions are again readily applied to river systems. Here we
consider two cases, one for vertical diffusion from the river into
its sediments at the bottom (Case F1), and one for radial diffusion from a semi-circular stream into the surrounding subsurface (Case F2).

centration in the sediment, Cs(x,z,t), is now also a function of
the vertical distance z below the river bottom. Advective
transport as well as longitudinal diffusive transport in the
hyporheic zone are ignored.
We present here solutions of the above hyporheic zonetransport/diffusion model subject to the initial and boundary
conditions

C(x,0) = Cs (x, z,0) = Ci ,⎛Dx ∂C ⎞⎜⎝ C − u ∂x ⎟⎠

⎧ Co
⎪
=⎨
x=0+ ⎪⎩0(33)

0 < t ≤ to(34)

t ≥ to ,∂C(∞, t) = 0 .
∂x

(35)

The analytical solution for the stream concentration is (vanGenuchten, 1985a)⎧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 ,

(36)

where

Fig. 3. Schematic of river system with a rectangular hyporheiczone.

A(x,t) =Case F1. Assuming a rectangular geometry of the river andhyporheic zone (Fig. 3), the transport/diffusion model may be
written as (see also Sudicky and Frind, 1982)∂C∂2 C∂C J z
= Dx
−u
− ,
2
∂t
∂x d
∂x
J z = − θ s Ds

Rs

Cs (x,0,t) = C(x,t),∂Cs (x, zo ,t)
= 0,
∂z

(

⎡⎛ u
⎞
• ⎢⎜
+ η p ⎟ sin γλ 2t − ηm x
2D
⎠
⎢⎣⎝
x

(

)

(37)

)

dλ
−ηm cos γλ 2t − ηm x ⎤⎥
⎦ λ

(28)

(29)

(0 < z ≤ zo ),

(30)

(31)

ηm =

Ωo − Ω12

Ωo =Ω12 + Ω22 ,

(39)

Ds

(40)

γ =

2d 2 Rsηp =Ωo + Ω1,
2,(38a,b)(32)

where C and Cs are solute concentrations in the stream andsediment, respectively, d is the depth of the river, z is verticaldistance below the river, Jz is the vertical solute diffusive flux
into the hyporheic zone (ML–2T–1), zo is the effective depth of
the sediment, θs is the volumetric water content of the sediments, Ds is the apparent solute diffusion coefficient in thesediment (L2T–1), and Rs is a solute retardation factor accounting for linear sorption/exchange in the sediment (–). Note that,as before, the stream solute concentration, C(x,t), is a functionof longitudinal distance, x, and time, t, but that the solute con-254

1 2u ∞+∫
22 π Dx 0 ⎛
⎞
u
2
+
η
⎜ 2D
p ⎟ + ηm⎝⎠x

in which∂Cs,∂z z=0∂Cs∂ 2 Cs
= Ds
∂t
∂z 2

⎛ ux
⎞
exp ⎜
− η pz⎟
⎝ 2Dx
⎠

Ω1 =

θ Dλ
u2
+ s2 s H1,
4Dx d D

(41a)

x

Ω2 =

Ds λ 22

2d Rs Dx

H1( λ ) =+

θ s Ds λd 2 DxH2 ,sinh( λ zo / d) − sin( λ zo / d),cosh( λ zo / d) + cos( λ zo / d)

(41b)
(42a)

Exact analytical solutions for contaminant transport in rivers: 2. Transient storage and decay chain solutions

H 2 (λ ) =sinh( λ zo / d) + sin( λ zo / d).
cosh( λ zo / d) + cos( λ zo / d)(42b)

The above solution holds for a third-type inlet boundary
condition. The solution for a first-type boundary condition isexactly the same, except that (37) must be replaced byA(x,t) =

(

)

⎛ ux
⎞
1 2 ∞
dλ
+ ∫ exp ⎜− η p z ⎟ sin γλ 2t − ηm x
. (43)
2 π 0
2D
λ
⎝
⎠
x

The same solutions also hold for the slightly simpler situation where diffusion in the sediments occurs over a semiinfinite region, i.e., zo→ ∞. Eqs. (42a,b) then reduce to unity,and Eqs. (41a,b) become

Ω1 =

θ Dλ
u2
+ s2 s ,4Dx d D

Rs

∂Cs Ds ∂ ⎛ ∂Cs ⎞=r∂tr ∂r ⎜⎝ ∂r ⎟⎠(a < r ≤ ro ) ,(49)

Cs (x,a,t) = C(x,t) ,(50)∂Cs (x,ro ,t)
=0,
∂r

(51)

where a and ro are as shown in Fig. 2, and r is the radial coordinate

(44a)

x

Ω2 =

Ds λ 22

2d Rs Dx+

θ s Ds λd 2 Dx.

(44b)

If dispersion in the river system is neglected (Dx→ 0), A(x,t)
for the first- and third type boundary conditions both reduce to(see also Skopp and Warrick, 1974)⎛ θ D xλ ⎞1
A(x,t) = + ∫ exp ⎜ − s s2 H1 ⎟2 π 0
⎝ 2d u⎠2 ∞Fig. 4. Schematic of river system with cylindrical hyporheic zone.

The solutions for this case is the same as for Case F1, except
for the following expressions (van Genuchten et al., 1984)

Ds

(45)

γ =

which holds for t > x/u. Eq. (45) and similar expressions forA(x,t) pertaining to transport problems which neglect dispersion
in the river system are understood to be zero for t < x/u. Eq.
(45) for zo → ∞ may be expressed in a much simpler alternative
form (Grisak and Pickens, 1981; Tang et al., 1981).

Ω1 =⎡ D λ 2 ( ut − x ) H θ D xλ
⎤ dλ
1
s s
⎥
i sin ⎢ s
−
H
2
⎢⎣
⎥⎦ λ
2d 2uRs
2d 2u⎡ θ x DR ⎤s s⎥
A(x,t) = erfc ⎢ s⎢⎣ 2d u 2t − ux ⎥⎦(x < ut) .

(46)

Alternative expressions for the solution of the abovetransport/diffusion problem for a first-type inlet boundary condition were given by Sudicky and Frind (1982) assuming afinite zo, and Tang et al. (1981) for the case where zo → ∞.
These two studies also presented solutions for the sedimentconcentration, Cs.
Case F2. This example is the same as F1, except that thestream has a semi-circular cross-sections as shown in Fig. 4.
Diffusion now takes place in a radial direction. For this case,
Eqs. (28) through (32) must be replaced by

J
∂C
∂2 C
∂C
= Dx
−u
− r2 ,
2
∂t
∂x
∂x
πa
J r = −2π aθ a Da

∂Cs,
∂r r=aa 2 Rs
u2
4Dx2

Ω2 =

,

+

2Ds λ 2ro2 Rs Dx(52)

2θ s Ds λa 2 Dx+

H1 ,2θ s Ds λa 2 Dx

(53a)

H2 ,

(53b)

where

H1 =

H2 =

N 1 ( M1 − M 2 ) + N 2 ( M1 + M 2 )N12 + N 22
N 1 ( M1 + M 2 ) − N 2 ( M1 − M 2 )N12 + N 22

,

(54a)

,

(54b)

M1( λ ) = Ber1(ξoλ )Ker1( λ ) − Bei1(ξ o λ )Kei1( λ )− Ker1(ξ o λ )Ber1( λ ) + Kei1(ξ o λ )Bei1( λ ),M 2 ( λ ) = Ber1(ξoλ )Kei1( λ ) + Bei1(ξ o λ )Ker1( λ )

(47)

(48)

− Ker1(ξ o λ )Bei1( λ ) − Kei1(ξ o λ )Ber1( λ ),N1( λ ) = Bei1(ξoλ )Ker( λ ) + Ber1(ξ o λ )Kei( λ )− Kei1(ξ o λ )Ber( λ ) + Ker1(ξ o λ )Bei( λ ),

(55a)

(55b)

(56a)

255

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

N 2 ( λ ) = Bei1(ξoλ )Kei( λ ) − Ber1(ξ o λ )Ker( λ )(56b)− Kei1(ξ o λ )Bei( λ ) − Ker1(ξ o λ )Ber( λ )in which ξo = ro/a, and where Ber, Bei, Ker and Kei represent
Kelvin functions (Olver, 1970).
For the case of a semi-infinite radial diffusion system (ro→∞), H1 and H2 reduce to2 ⎡⎣ Ker( λ ) Ker '( λ ) + Kei( λ ) Kei'( λ ) ⎤⎦H1( λ ) = −H 2 (λ ) =Ker 2 ( λ ) + Kei 2 ( λ )2 ⎡⎣ Kei( λ ) Ker '( λ ) − Ker( λ ) Kei'( λ ) ⎤⎦Ker 2 ( λ ) + Kei 2 ( λ ),.

where Ci represents the concentration of the i-th species, and µi
is the i-th first-order degradation coefficient. Eqs. (62a,b) aresolved for the initial and boundary conditions

Ci (x, 0) = 0(i = 2,3,4) , (62b)(i = 1,4) ,(63)⎧ g i(t)
⎪
= ⎨
⎪⎩0
x=0+

0 < t ≤ to(64)t ≥ to ,∂Ci(∞, t) = 0
∂x

(65)

in which gi(t) are the prescribed input concentration functions.
We give here the solution for the very general situation where
the input concentrations are given by

g1(t) = B1 e

− λ1t

(x < ut) , (59)g3 (t) = B4 e

,

(66a)

− λ1t

+ B3 e− λ1t+ B5 e− λ1t+ B8 eg2 (t) = B2 ewhich was found to give accurate results for the condition

− λ2t

,

(66b)

− λ2t

+ B6 e

− λ3t

− λ2t

+ B9 e

,

(66c)

2

ua Rs
Ds

.

(60)

The above general solutions for rectangular and cylindricalstreambeds involve integrals of expressions that are the productof decaying exponential functions and rapidly oscillating sinusoidal functions. Because of the oscillatory properties of theintegrands as a function of the integration parameter λ, direct
numerical integration of the integrals using Gaussian quadrature or related techniques often leads to very poor results, evenwith an excessive number of quadrature terms. We refer to
Rasmuson and Neretnieks (1981) and van Genuchten et al.
(1984) for efficient procedures to evaluate the integrals.
CONSECUTIVE DECAY CHAINS
The final set of analytical solutions in this paper concerns
the movement of solutes involved in sequential first-orderdecay reactions. Typical examples are the transport of interacting nitrogen species, pesticides and their degradation products,
chlorinated hydrocarbons, organic phosphates, pharmaceuticals,
radionuclides, and interacting toxic trace elements (e.g., fromacid mine drainage or other anthropogenic or natural sources).
We consider here the transport of four species (Ei, i = 1,..,4)
involved in the consecutive decay chain

256

(62a)

(57b)

which also hold for ro →∞, provided H1 and H2 are given by(57a,b). A good approximation of (58) for ro →∞ is (vanGenuchten et al., 1984)x ≤ ut < x +

∂C1
∂2 C1
∂C
= Dx
− u 1 − µ1C1 ,
∂t
∂x
∂x 2
∂Ci
∂ 2 Ci∂C= Dx− u i + µ i−1Ci−1 − µ iCi
2
∂t
∂x
∂x

(58)

⎡ D λ 2 ( ut − x )⎤ dλ
2 θ s Ds xλi sin ⎢ s 2
−
H2 ⎥
2
⎢⎣
⎥⎦ λ
a uRs
a u

⎡θ x D R ⎤⎛ θ D x⎞s s ⎥A(x,t) = exp ⎜ − s 2 s ⎟ erfc ⎢ s⎝⎢⎣ a u 2t − ux ⎥⎦a u ⎠The transport equations for this problem are⎛Dx ∂Ci ⎞⎜⎝ Ci − u ∂x ⎟⎠⎛2 θ s Ds xλ ⎞1 2 ∞+ ∫ exp ⎜ −
H1 ⎟
2 π 0
a 2u
⎝
⎠

(61)

(57a)

The solutions again simplify when longitudinal dispersion in
the river is assumed to be negligible (Dx = 0). A(x,t) in Eq. (37)becomes thenA(x,t) =E1 → E2 → E3 → E4 .

g4 (t) = B7 e

− λ3t

+ B10 e

− λ4t

,

(66d)

where Bi and λi are constants. The multiple terms in Eqs. (66)
may describe possible decay processes in the contaminantsource, and/or account for finite release rates from the source
into the river system. For one particular release mechanism, theconstants Bi are related to each other through the Bateman
equations (Bateman, 1910). A detailed description of this situation is given by van Genuchten (1985b).
The analytical solution of the above transport problem is(van Genuchten, 1985b)⎧C * (x,t)
0 < t ≤ to
⎪ i
Ci (x,t) = ⎨⎪C * (x,t) − e− λito C * (x,t − t ) t > t ,
i
o
o
⎩ i

(67)

where

C1* = B1F110 ,(68)

Exact analytical solutions for contaminant transport in rivers: 2. Transient storage and decay chain solutions

C2* = B2 F210 + B3 F220 +µ1B1(S − F210 + F110 ) ,µ2 − µ1 12(69)

C3* = B4 F310 + B5 F320 + B6 F330 ++µ2 B2(S − F + F210 )µ3 − µ2 23 310

+

µ2 B3
(S − F + F220 ),µ3 − µ2 23 320(70)

Here we give one illustrative hypothetical application of theconsecutive decay chain solution given by Eqs. (67)–(76). Fig.
5 shows calculated distributions versus distance for a threemember decay chain (E1→E2→E3) assuming a finite pulse typeinjection (to = 500 s) of solute E1 with concentration Co=1 g/m3,
using mostly the same parameters as before for Fig. 1 (i.e., u =1 m/s and Dx = 5 m2/s), and with values of 0.004, 0.001 and
0.002 s–1 for the degradation coefficients of solutes E1, E2 andE3, respectively.

C4* = B7 F410 + B8 F420 + B9 F430 + B10 F440 ++µ3 B4(S − F410 + F310 ) +µ4 − µ3 34

+

µ3 B5
(S − F420 + F320 ) +µ4 − µ3 34

+

(71)

µ3 B6
(S − F430 + F330 ),µ4 − µ3 34

Sij = F jji − Fiji ,

(72)

⎧ u
⎡ x − ξt ⎤⎡ (u − ξ )x ⎤⎪⎥Fijk = exp(−aijk t) ⎨exp ⎢⎥ erfc ⎢⎢⎣ 4Dx t ⎥⎦⎣ 2Dx ⎦⎪⎩ u + ξ⎡ x + ξt ⎤⎡ (u + ξ )x ⎤
u
⎥
+
exp ⎢
(73)
⎥ erfc ⎢u−ξ⎢⎣ 4Dx t ⎥⎦⎣ 2Dx ⎦
⎡ x + ut ⎤ ⎫⎪
⎛ ux
⎞
2u 2
⎥⎬
+ 2
exp ⎜
− µit ⎟ erfc ⎢( µi ≠ aijk ),2⎢⎣ 4Dx t ⎥⎦ ⎪
⎝ Dx
⎠
ξ −u
⎭

⎧1
⎡ x − ut ⎤⎡ (x − ut)2 ⎤
u 2t
⎪
⎥+
Fijk = exp(−aijk t) ⎨ erfc ⎢exp ⎢ −⎥π Dx
4Dx t ⎥⎦⎢⎣ 4Dx t ⎥⎦⎢⎣⎪⎩ 2
⎡ x + ut ⎤ ⎫⎪⎛ ux ⎞
1⎛ux u 2t ⎞⎢⎥ ⎬ ( µi = aijk ),− ⎜1 +
+
exp
erfc
⎟
⎜D ⎟
2⎝
Dx Dx ⎠⎢⎣ 4Dx t ⎥⎦ ⎪
⎝ x⎠
⎭
(74)

⎧⎪ λ jaijk = ⎨µ − µj
⎩⎪ i

k=0
k >0

ξ = u 2 + 4Dx ( µi − aijk ) .

(75)

(76)

The solution for a first-type inlet boundary condition is thesame as the above third-type inlet solution, except for the term
Fijk (Eqs. (73) and (74)), which should be replaced by⎧1⎡ x − ξt ⎤⎡ (u − ξ )x ⎤⎪⎥Fijk = exp(−aijk t) ⎨ exp ⎢⎥ erfc ⎢⎢⎣ 4Dx t ⎥⎦⎣ 2Dx ⎦⎪⎩ 2
⎡ x + ξ t ⎤ ⎫⎪⎡ (u + ξ )x ⎤1
⎥ ⎬.
+ exp ⎢
⎥ erfc ⎢2⎢⎣ 4Dx t ⎥⎦ ⎪⎣ 2Dx ⎦⎭

(77)

Fig. 5. Calculated distributions versus distance for a three-member
decay chain (E1→E2→E3) upon injection of a finite pulse (to = 500s) of solute E1 having a concentration Co of 1 g/m3 (u = 1 m/s, Dx =5 m2/s, µ1 = 0.004 s–1, µ2 = 0.001 s–1, and µ3 = 0.002 s–1).

A useful computer program for evaluating the above consecutive chain transport solutions is the CHAIN code of van
Genuchten (1985b), also incorporated into STANMOD. We
note that the decay chain solutions in this section assume nosorption of all species involved. The CHAIN code is for the
more realistic case where sorption can occur, and where theretardation factors (Ri) of the individual species can be different. Finally, we note similar decay chain solutions for multidimensional transport with unequal retardation factors, but for
first-type boundary conditions and without the Bateman source
equations, are discussed in a recent paper by Quezada et al.
(2004).
CONCLUDING REMARKSIn this two-part paper we collected a large number of analytical solutions for contaminant transport in rivers and surface
water bodies. The solutions in part 1 pertain to one-dimensional
longitudinal transport in streams and rivers, and longitudinaltransport and lateral dispersion in rivers and larger surface
water bodies. The current part 2 focused on nonequilibriumtransport caused by the presence of stagnant water zones (transient zone models), and simultaneous longitudinal advectivedispersive transport in a river and diffusion into and out of thehyporheic zone. We also provided several solutions for thetransport of solutes involved in consecutive decays chains.
Most of the solutions were derived from solutions to mathematically very similar problems in subsurface contaminanttransport. Except for solutions pertaining to diffusion in
hyporheic zones, all solutions have been incorporated in the
public-domain windows-based STANMOD software package(Šimůnek et al., 2000). This software package also includes
parameter estimation capabilities, and hence may be a convenient tool for analyzing observed contaminant concentration

257

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

distributions versus distance and/or time. While inherently lessflexible than more comprehensive numerical models for contaminant transport in streams and rivers, we believe that theanalytical solution can be very useful for simplified analyses ofalternative contaminant transport scenarios, as well as for testing of numerical models.
REFERENCES

Bateman, H., 1910. The solution of a system of differential
equations occurring in the theory of radioactive transformations. Proc. Cambridge Philos. Soc., 15(V), 423–427.
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., Walters, R., 1983. Simulation of transport in a
mountain pool-and-riffle stream: A transient storage model.
Water Resour. Res., 19, 718–724.
Brevis, W., Debels, P., Vargas, J., Link, O., 2001. Comparisonof methods for estimating the longitudinal dispersion coefficient of the Chillán River, Chile. (Comparación de Métodos
para estimar el coeficiente de dispersión longitudinal en el
Río Chillán, Chile.) Anales de XV Congreso Chileno deIngenieriá Hidráulica, Sociedad Chilena de Ingenieriá
Hidráulica, Concepción, Chile, 155–164. (In Spanish.)Coats, K.H., Smith, B.D., 1964. Dead-end pore volume anddispersion in porous media. Soc. Petrol. Eng. J., 4, 73–84.
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., 2006. Analytical solution for transport of decaying solutes in rivers with transient storage. J. Hydrol., 330,
672–680.
De Smedt, F., 2007. Analytical solution and analysis of solutetransport in rivers affected by diffusive transfer in the
hyporheic zone. J. Hydrol., 339, 29–38.
Field, M.S., 1997. Risk assessment methodology for karst aquifers: (2) Solute transport modeling. Env. Modeling and Assessment, 47, 23–37.
Goldstein, S., 1953. On the mathematics of exchange processes
in fixed columns. I. Mathematical solutions and asymptotic
expansions. Proc. R. Soc. London Ser., A. 219, 151–185.
Grant, S.B., Stewardson, M.J., Marusic, I., 2012. Effective
diffusivity and mass flux across the sediment-water interface
in streams. Water Resour. Res., 48, W05548, doi:
10.1029/2011WR011148.
Grisak, G.E., Pickens, J.F., 1981. Analytical solution for solutetransport through fractured media with matrix diffusion. J.
Hydrol., 52, 47–57.
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.
Jackman, A.P., Walters, R.A., Kennedy, V.C., 1984. Transportand concentration controls for chloride, strontium, potassium
and lead in Uvas Creek, a small cobble-bed stream in Santa
Clara County, California, USA., 2. Mathematical modeling.
J. Hydrol., 75, 111–141.
LeGrand-Marcq, C., Laudelout, H., 1985. Longitudinal dispersion in a forest stream. J. Hydrol., 78, 317–324.
Leij, F.J., Toride, N., 1995. Discrete time- and length-averaged
solutions of the advection-dispersion equation. Water Resour. Res., 31, 1713–1724.

258

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., Toride, N., van Genuchten, M.Th., 1993. Analyticalsolutions for non-equilibrium solute transport in threedimensional porous media. J. Hydrol., 151, 193–228.
Leij, F. J., van Genuchten M.Th., 2000. Analytical modeling of
nonaqueous phase liquid dissolution with Green’s functions.
Transport in Porous Media, 38, 141–166.
Olver, F.W.J., 1964. Bessel functions of integer order. In:
Abramowitz, M. Stegun, I.A., Handbook of Mathematical
Functions, 355–433. Appl. Math. Ser. 55. National Bureauof Standards, Washington, D.C.
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.
Quezada, C.R., Clement, T.P., Lee, K.-K., 2004. Generalized
solutions to multi-dimensional multi-species transport equations coupled with a first-order reaction network involving
distinct retardation factors. Adv. Water Resour., 27, 507–
520.
Rasmuson, A., Neretnieks, I., 1981. Migration of radionuclides
in fissured rock: The influence of micropore diffusion andlongitudinal dispersion. J. Geophys. Res., 86(B5), 3749–
3758.
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, 211–215.
Runkel, R.R., McKnight, D., Rajaram, H. (eds.). 2003. Modeling hyporheic processes. Special issue. Adv. Water Resour.,26(9), 901–1040.
Š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.hydrus2d.com.
Skopp, J., Warrick, A.W., 1974. A two-phase model for the
miscible displacement of reactive solutes in soil. Soil Sci.
Soc. Am. J., 38(4), 545–550.
Sudicky, E.A., Frind, E.O., 1982. Contaminant transport infractured porous media: Analytical solution for a system ofparallel fractures. Water Resour. Res., 18, 1634–1642.
Tang, D.H., Frind, E.O., Sudicky, E.A., 1981. Contaminanttransport in fractured porous media. Water Resour. Res., 17,555–564.
Thackston, E.L., Schnelle Jr., K.B., 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.
Toride, N., Leij, F.J., van Genuchten, M.Th., 1993. A comprehensive set of analytical solutions for nonequilibrium soluteExact analytical solutions for contaminant transport in rivers: 2. Transient storage and decay chain solutionstransport with first-order decay and zero-order production.
Water Resour. Res., 29(7), 2167–2182.
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. 1985a. A general approach for modeling solute transport in structured soils. Proc. 17th Int. Congress, Hydrogeology of Rocks of Low Permeability, Jan. 7–
11, Tucson, Arizona. Memoires Int. Assoc. Hydrogeol.,
17(2), 513–526.
Van Genuchten, M.Th., 1985b. Convective-dispersive transportof solutes involved in sequential first-order decay reactions.
Computers & Geosciences, 11(2), 129–147.
Van Genuchten, M.Th., Wagenet, R.J., 1989. Two-site/tworegion models for pesticide transport and degradation: Theoretical development and analytical solutions. Soil Sci. Soc.
Am. J., 53(5), 1303–1310.
Van Genuchten, M.Th., Wierenga, P.J., 1976. Mass transfer
studies in sorbing porous media. I. Analytical solutions. Soil
Sci. Soc. Am. J., 40(4), 473–480.
Van Genuchten, M.Th., Gorelick, S.M., Yeh, W.-G., 1988.
Application of parameter estimation techniques to solutetransport studies. In: Donn G. DeCoursey (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., Tang, D.H., Guennelon, R., 1984.
Some exact solutions for solute transport through soils containing large cylindrical macropores. Water Resour. Res.,
20(3), 335–346.
Van Genuchten, M.Th., Leij, F.J., Skaggs, T.H., Toride, N.,
Bradford, S.A., Pontedeiro, E.M., 2013. Exaxt analytical solutions for contaminant transport in rivers. 1. The equilibrium advection-dispersion equation. J. Hydrol. Hydromech.,61(2), 146–160.
Wörman, A., 1998. Analytical solution and timescale fortransport of reacting solutes in rivers and streams. Water Resour. Res., 34(10), 2703–2716.
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.
Wlostowski, A.N., Gooseff, M. Wagener, T., 2013. Influence ofconstant rate versus slug injection experiment type on parameter identificability in a 1-D transient storage model for
stream solute transport. Water Resour. Res., 49(1-5), doi:
10.1002/wrcr.20103.
Zwillinger, D., 1989. Handbook of Differential Equations.
Academic Press, San Diego, CA.
Received 8 August 2012
Accepted 23 February 2013259