Transp Porous Med (2009) 80:373–387
DOI 10.1007/s11242-009-9368-3

Analytical Solution for Multi-Species Contaminant
Transport Subject to Sequential First-Order Decay
Reactions in Finite Media
Jesús S. Pérez Guerrero · Todd H. Skaggs ·M. Th. van Genuchten

Received: 22 December 2008 / Accepted: 18 February 2009 / Published online: 19 March 2009
© Springer Science+Business Media B.V. 2009

Abstract Transport equations governing the movement of multiple solutes undergoing
sequential first-order decay reactions have relevance in analyzing a variety of subsurface
contaminant transport problems. In this study, a one-dimensional analytical solution for
multi-species transport is obtained for finite porous media and constant boundary conditions.
The solution permits different retardation factors for the various species. The solution procedure involves a classic algebraic substitution that transforms the advection-dispersion partial
differential equation for each species into an equation that is purely diffusive. The new systemof partial differential equations is solved analytically using the Classic Integral Transform
Technique (CITT). Results for a classic test case involving a three-species nitrification chain
are shown to agree with previously reported literature values. Because the new solution was
obtained for a finite domain, it should be especially useful for testing numerical solution
procedures.
Keywords Multi-species transport · Finite domain · Analytical solution ·Integral transform
List of Symbols
b0 j , b1 j Auxiliary coefficients
C0
Dimensional reference solute concentration

J. S. Pérez Guerrero
Radioactive Waste Division, Brazilian Nuclear Energy Commission (DIREJ/DRS/CNEN),
Rio de Janeiro, Brazil
e-mail: jperez@cnen.gov.br
T. H. Skaggs (B)U.S. Salinity Laboratory, USDA-ARS, Riverside, CA, USAe-mail: Todd.Skaggs@ars.usda.gov
M. Th. van Genuchten
Department of Mechanical Engineering, LTTC/COPPE, Federal University of Rio de Janeiro, UFRJ,
Rio de Janeiro, Brazil
e-mail: rvangenuchten@yahoo.com

123

374

Cj
C0 j
c j (x, t)
c0 j
C1 j , C2 j
D
Fj
f¯ji
G j (X )g j (x)
H1 , H2i, j, k
L
Ni
p
Pe
qj
Rj
S̄ ji (t)Tt
T j (x, t)
U
wj
X
x

J. S. Pérez Guerrero et al.

Dimensional solute concentration of the jth species
Dimensional inlet boundary concentration of the jth species
Dimensionless solute concentration of the jth species
Dimensionless inlet boundary concentration of the jth species
Auxiliary coefficients
Dispersion coefficient
Filter function
Integral coefficient
Dimensional initial concentration of the jth species
Dimensionless initial concentration of the jth species
Coefficients
Indices
Domain length
Norm
Constant used in algebraic substitution
Peclet number
Constant used in algebraic substitution
Retardation coefficient for the jth species
Integral coefficient
Dimensional time
Dimensionless time
Unknown function
Constant pore water velocity
Auxiliary coefficient
Dimensional spatial coordinate
Dimensionless spatial coordinate

Greek Symbols
Eigenvalue
βiγjDamkholer number
δik
Kronecker deltaθ j (x, t) Unknown function in purely diffusive equation
Integral transform of the function θ j (x, t)θ̄ ji (t)λjFirst-order decay constant for the jth speciesjCoefficientµi
EigenvalueτEigenvalue
φ j (x)
Auxiliary function
ψi (x)
Eigenfunctionψ̃i (x)
Normalized eigenfunction
Auxiliary functionj (x)1 Introduction
Transport equations governing the movement of multiple solutes undergoing sequential firstorder decay reactions have been studied extensively for several decades because of their123

Analytical Solution for Multi-Species Contaminant Transport

375

relevance in modeling the subsurface transport of a variety of contaminants, including radionuclide decay chains, pesticides and their degradation products, petroleum hydrocarbons, andnitrogen species. Early studies such as those by McLaren (1969) and Cho (1971) focused on
obtaining analytical solutions for the transport equations, whereas more recently numericalsolutions have become increasingly popular using widely available computer software (e.g.,ºŠuimunek
2006). Still, considerable interest remains in analytical solution because they may
be used as relatively simple screening models in regulatory decision making (e.g., Aziz et al.
2000; Skaggs et al. 2007), and because analytical solutions may serve as benchmarks fortesting numerical solutions.
A variety of analytical methods have been used to solve sequential decay transport problems. Cho (1971) used Laplace transforms to obtain an analytical solution for convective-dispersive nitrogen transport involving a three-species nitrification chain. van Genuchten (1985)
similarly used Laplace transforms to obtain the solution for a four-species decay chain, andimplemented the solution in a computer program called CHAIN containing example problems dealing with nitrification–denitrification and radionuclide decay chains. In more recent
work, Lunn et al. (1996) presented analytical solutions obtained with Fourier sine transforms
using a methodology that enabled greater flexibility with respect to initial conditions. Theseauthors were able to produce results identical to those of Cho (1971), while also developing
new solutions for constant and exponential initial contaminant profiles.
A somewhat different approach was followed by Sun et al. (1999a,b) and Sun and Clement (1999) who used an ingenious substitution to transform the multidimensional first-ordersequential decay problem with any number of species into a decoupled set of transport equations for single species. They then solved the decoupled equations using existing analyticalsolutions for single-species transport and first-order decay. Unfortunately, the substitution
can be used only when the retardation coefficients are the same for all species.
Bauer et al. (2001) subsequently presented a method for solving one-, two-, and threedimensional sequentially coupled reactive transport problems with distinct retardation coefficients. Their solution was expressed in terms of a linear combination of fundamental solutionsin the Laplace domain, with coefficients defined by a recursive formula. In some cases it waspossible to transform the solution back to the time domain analytically. However, most transport problems required a numerical approach to obtain the inverse transform. Sun et al. (2001)
avoided numerical inversion of the Laplace transforms by applying the solution scheme of
Bauer et al. (2000) directly in the time domain instead of the Laplace domain. The approachrelied heavily on the use of fundamental analytical solutions for single-species transport
subject to first-order decay.
Clement (2001) wrote the equations for coupled species transport using matrix notationand then used a matrix diagonalization procedure to decouple the equations. Their solution method initially required that different species have identical retardation coefficients.
Quezada et al. (2004) later extended the procedure of Clement (2001) and presented a general
method for solving coupled multidimensional, multi-species reactive transport equations withdifferent retardation factors. The method employs a three-step transformation that involves
Laplace transforms and matrix diagonalization in the Laplace domain. Quezada et al. (2004)
presented an analytical solution for two-species transport; more general cases could be solved
only semi-analytically as they require numerical inversion of the Laplace transform. In related
work, Lu and Sun (2008) developed a model that accounted for equilibrium and kinetic
reactions. They decoupled the transport equations using the matrix diagonalization approach,and then obtained closed-form analytical solutions.
Most recently, Srinivasan and Clement (2008a,b) developed closed-form analytical solutions for sequential decay problems involving an arbitrary number of species and subject

123

376

J. S. Pérez Guerrero et al.

to spatially varying initial conditions and a generic exponentially decaying Bateman-typesource boundary condition as was also used previously by van Genuchten (1985). The coupled one-dimensional system of equations was solved for both Dirichlet and Cauchy boundary
conditions using a combination of Laplace transform and the matrix diagonalization technique of Clement (2001).
Several things are apparent from the above literature survey. First and foremost, the review
shows that for over three decades there has been steady interest in analytical methods forsolving contaminant transport problems with sequential decay reactions, and that interest
continues up to the present. Also, many different analytical methods have been used, with
Laplace transforms, various decoupling strategies, and matrix diagonalization techniques
being the most prevalent. However, one commonality among all methods thus far is that they
have produced analytical solutions only for infinite or semi-infinite spatial domains. Given
that a primary use of analytical solutions is to test numerical models that compute solutions
on finite domains, it would be very useful to also have analytical solutions for finite domains.
Pontedeiro et al. (2007) previously obtained solutions on a finite domain using a semianalytical approach to perform sensitivity analyses and risk assessments of technologicallyenhanced, naturally occurring radioactive material (TENORM) disposed in an industrial
landfill. They solved the decay chain problem using a semi-analytic approach known as
Generalized Integral Transform Technique (GITT), which is a generalization of the Classic
Integral Transform Technique, or CITT (Cotta 1993). Very recently, Cassol et al. (2009)
combined GITT and CITT techniques with Laplace transforms and matrix diagonalization
methods to develop an analytical solution for transient two-dimensional atmospheric pollutant dispersion on a semi-infinite spatial domain.
Using appropriate substitutions, many transport problems can be transformed into a form
that is amenable to analytical solution via the CITT. Pérez Guerrero et al. (2009) recently
demonstrated the latter approach by obtaining the solution of a general linear advectiondispersion equation for single-species transport. The objective of the present study was todevelop the analytical solution for a multi-species advection-dispersion transport equation
comprising sequential first-order decay reactions, distinct retardation coefficients for eachspecies, a finite spatial domain, and steady-state boundary conditions.

2 Problem Formulation
Consider a one-dimensional, homogeneous porous medium having a constant saturated or
unsaturated water content. For multiple reactive solute species that are subject to linear equilibrium sorption and are involved in a sequential first-order decay chain, the governing systemof transport equations is given by:

Rj

∂C j
∂C j
∂ 2C j
−U
=D
∂T
∂ X2
∂X
−λ j R j C j + λ j−1 R j−1 C j−1 ;j = 1, 2, 3, . . . ; λ0 = 0(1a)where C j is the concentration of the jth species, R j is the retardation coefficient for thejth species, D is the dispersion coefficient, U is the constant pore water velocity, λ j is thefirst-order decay constant for the jth species, X is the spatial coordinate, and T is time.
The problem definition is completed by specifying for each species the initial and boundary
conditions on a finite domain of length L as follows:123Analytical Solution for Multi-Species Contaminant Transport377

C j (X, 0) = G j (X )∂C j (0, T )+ U C j (0, T ) = U C0 j
−D
∂X
∂C j (L , T )
=0
∂X

(1b)
(1c)
(1d)

We define the following non-dimensional variables:
x=

Cj
C0 j
X
T
; c0 j =; t=; cj =
L(L/U )
C0
C0(2a–d)

λj Rj LGj
UL; γj =; g j (x) =D
U
C0(2e–g)

Pe =where C0 is a reference concentration, Pe is termed the Peclet number, and γ j is the ratioof the reaction and advection time scales, often referred to as the Damkholer number. With
these non-dimensional variables, Eq. 1 may be rewritten as:
Rj

∂c j
∂c j
1 ∂ 2c j−=− γ j c j + γ j−1 c j−1 ; j = 1, 2, 3, . . . ; γ0 = 0∂t
Pe ∂ x 2∂xc j (x, 0) = g j (x)∂c j (0, t)−+ Pe c j (0, t) = Pe c0 j∂x
∂c j (1, t)
=0
∂x

(3a)
(3b)
(3c)
(3d)

3 Solution Methodology
In this study, a solution of Eq. 3 was developed using the CITT (Ozisik 1980; Mikhailov
and Ozisik 1984; Cotta 1993) in combination with a classical algebraic transformation thattransforms the governing advection-dispersion equations into purely diffusive equations. The
CITT can be applied to Eq. 3 as follows:(a) define and solve an auxiliary eigenvalue problem related to the original problem;(b) develop an appropriate integral and inverse transform pair;(c) transform each of the partial differential equations such that the space variable is eliminated and an infinite system of decoupled ordinary differential equations is obtained;(d) solve the resulting ordinary differential equation system subject to the transformed
initial condition; and(e) use the inverse transform to obtain the desired solution.
The CITT can be applied to transport problems with homogeneous or non-homogeneous
boundary conditions. When the problem has homogeneous boundary conditions, the eigenfunctions resulting from the auxiliary eigenvalue problem in step (a) automatically include
boundary information. For the case of non-homogeneous boundary conditions, it is necessary
to use the Green’s theorem to incorporate the boundary information, as indicated by Ozisik(1980). However, while the application of the Green’s theorem is formally correct, its use
leads to very slow convergence near the boundary, as illustrated by Cotta (1993). An alternative is to homogenize the boundary conditions using a “filter” function. The idea is to define
a new problem with the non-homogeneous boundary conditions of the original problem in

123

378

J. S. Pérez Guerrero et al.

order to obtain a “filter” function capable of homogenizing the boundary conditions of theoriginal problem. We homogenize the boundary conditions using a filter function obtained
from a steady-state version of the transport equation. Additionally, we utilize an algebraic
substitution to transform the advection-dispersion transport equations into exclusively diffusive equations that are easier to solve analytically.
3.1 The Steady-State Regime and the Filter Function
A filter function capable of homogenizing the boundary condition can be obtained in different
ways. Here, we obtain a filter function from the steady-state version of Eq. 3. We define forthis purpose the following problem:dF j
1 d2 Fj−− γ j F j + γ j−1 F j−1 = 0; j = 1, 2, 3, . . .
Pe dx 2dx
dF j (0)+ Pe F j (0) = Pe c0 j ;−dx
dF j (1)
=0
dx

(4a)
(4b)
(4c)

This is an ordinary differential equation for the function F j which has an analytical solutionthat depends upon F j−1 . That dependence makes obtaining the analytical solution elaborate
for large values of j.
An interesting strategy for solving Eq. 4 is to use the procedure presented by Sun et al.
(1999a,b) and Sun and Clement (1999). As noted in the Introduction, their procedure applies
only when the retardation factors are the same for each species. In this study, we consider the
more general case in which each species has its own retardation factor as indicated by Eqs.
1 and 3. However, we can use the procedure of Sun et al. (1999a,b) and Sun and Clement(1999) to solve the filter problem as Eq. 4 is independent of the retardation factors.
Following the procedure of Sun et al. (1999a,b) and Sun and Clement (1999), we define
the following functions:
F j (x) = φ j (x) − j (x)
⎞
⎛
j−1
j−1


γ
m
⎠
Fk (x) ⎝j (x) =
γm − γ j
k=1

(5)
(6)

m=k

Substitution of Eqs. 5 and 6 into Eq. 4 leads to new equations involving φ j (x) that do nothave recursive source terms.
dφ j
1 d 2φ j− 2 − γjφj = 0
Pe dx 2dx
dφ j (0)d j (0)−+ Pe φ j (0) = b0 j = Pe c0 j −
+ Pe
dx
dx
dφ j (1)d j (1)= b1 j =
dx
dx

123

(7a)
j (0 )(7b)(7c)

Analytical Solution for Multi-Species Contaminant Transport

379

The general analytical general solution for φ j (x) is:√√φ j (x) = C1 j exp[x(Pe − Pe w j )/2] + C2 j exp[x(Pe + Pe w j )/2]

w j = Pe + 4γ j√√
2 exp[(w j Pe − Pe)/2]
C1 j ={b1 j (w j − Pe)j√√+ exp[(w j Pe + Pe)/2)]b0 j (w j + Pe)}√
2 exp(−Pe/2)
C2 j =[exp(Pe/2)b0 j (w j − Pe)j√√+ exp(w j Pe/2)b1 j (w j + Pe)]√√√√ j = Pe{[exp(w j Pe) − 1]Pe + 2[exp(w j Pe) + 1]w j Pe√
+ [exp(w j Pe) − 1]w 2j }

(8a)
(8b)

(8c)

(8d)
(8e)

The filter function F j (x) can now be obtained immediately using Eqs. 5, 6, and 8.
3.2 Homogenization of the Boundary Condition and Use of Algebraic Substitution
Equation 3 can be re-written as a homogeneous boundary condition problem by expressing
the dimensionless concentration c j (x, t) as a linear combination of the filter function and a
new unknown function T j (x, t) as follows:c j (x, t) = T j (x, t) + F j (x)(9)

Substituting Eq. 9 into Eq. 3 gives the following equations for the unknown field T j (x, t):∂Tj∂Tj
1 ∂2Tj=− γ j T j + γ j−1 T j−1 ; j = 1, 2, 3, . . .
−
(10a)
∂t
Pe ∂ x 2∂xNote that in this expression all terms involving F j vanished due to Eq. 4a. The initial andboundary conditions are respectively:Rj

T j (x, 0) = g j (x) − F j (x)(10b)∂ T j (0, t)(10c)−+ Pe T j (0, t) = 0∂x∂ T j (1, t)
=0
(10d)
∂x
The structure of Eq. 10 is the same as that of Eq. 3, with no additional source terms being created by the homogenization procedure. As shown in Appendix A, it is possible to transformEq. 10 from an advection-diffusion equation into a diffusion equation using the following
algebraic substitution (Ozisik 1980; Pérez Guerrero et al. 2009):
T j (x, t) = θ j (x, t) exp( px + q j t)
where:
p=1
Pe; qj = −
2
Rj



(11)

Pe
+ γj
4

(12)

The substitution produces the diffusion equation:Rj∂θ j
1 ∂ 2θ j+ γ j−1 T j−1 (x, t) exp(− px − q j t);=∂t
Pe ∂ x 2

j = 1, 2, 3, . . .

(13a)

123

380

J. S. Pérez Guerrero et al.

The initial and boundary condition for the unknown function θ j (x, t) are, respectively:θ j (x, 0) = exp(− px)[g j (x) − F j (x)]∂θ j (0, t)Pe−+θ j (0, t) = 0∂x2∂θ j (1, t)Pe+θ j (1, t) = 0
∂x
2

(13b)
(13c)
(13d)

Note that both boundary conditions (Eq. 10c, 10d) were modified as a consequence of thealgebraic substitution given by Eq. 11.
3.3 Solution via CITT
The problem defined by Eq. 13b can be solved using the CITT. Following the procedure outlined above, we select a Sturm–Liouville type problem as the auxiliary eigenvalue problem:d 2 ψi (x)+ βi2 ψi (x) = 0;
dx 2
µi2
βi2 =
(1/Pe)

(14a)
(14b)

with boundary conditions:−dψi (0)Pe+ψi (0) = 0dx
2dψi (1)Pe+ψi (1) = 0
dx
2

(14c)
(14d)

where ψi (x) is the eigenfunction and βi the eigenvalue. The analytical solution of this problem is (Ozisik 1980):
ψi = βi cos(βi x) + H1 sin(βi x)(15)The norm and eigenvalues can be calculated from the following equations:(βi2 + H12 ) + H1 + H22βi (H1 + H2 )tan(βi ) = 2
βi − H1 H2

Ni =

(16)
(17)

Pe
respectively, where H1 = Pe
2 and H2 = 2 .
The normalized eigenfunction and the orthogonality property are, respectively:ψ̃i (x) =ψi (x)
1/2

(18)

Ni

1

ψ̃i (x)ψ̃k (x)dx = δik0where δik is the Kronecker delta.

123

(19)

Analytical Solution for Multi-Species Contaminant Transport

381

At this point, we represent the unknown function θ j (x, t) as a series expansion in termsof the eigenfunctions:θ j (x, t) =∞ψ̃i (x)θ̄ ji (t) (I nver se)

(20a)

i=1

where θ̄ ji (t) is the transformed potential. An explicit definition for the transformed potential
1
is obtained by applying the operator 0 ψ̃i (x)dx to both sides of Eq. 20a and then invoking
the orthogonality property given by Eq. 19 (Ozisik 1980; Cotta 1993). This results in:1θ̄ ji (t) =ψ̃i (x)θ j (x, t)dx (T rans f or m)

(20b)

0

Equation 20b is the formula of the integral transformation which we use to transformEq. 13b into a system of ordinary differential equations with dependent variables θ̄ ji (t).
Once that system is solved, θ j (x, t) is obtained using the inverse transform, Eq. 20a.
1
Thus, we now apply the operator 0 ψ̃i (x)dx to Eq. 13b and invoke the integral transform defined by Eq. 20b to obtain the following system of decoupled ordinary differential
equations:Rjd θ̄ ji (t)+ µi2 θ̄ ji (t) = S̄ ji (t); i, j = 1, 2, 3, . . .
dtθ̄ ji (t = 0) = f¯ji

(21a)
(21b)

1

f¯ji =

ψ̃i (x) exp(− px)[g j (x) − F j (x)]dx

(21c)

0

The integral coefficient f¯ji can be evaluated analytically in some cases; in others, it is necessary to calculate this coefficient numerically within a prescribed error tolerance. The termS̄ ji (t) is the transformed source term and is defined by the integral:1S̄ ji (t) =ψ̃i (x)γ j−1 T j−1 (x, t) exp(− px − q j t)dx

(22a)

0

This expression can be simplified to the form:S̄ ji (t) = γ j−1 θ̄ j−1 i (t) exp[(q j−1 − q j )t]
Equation 21 is an ordinary differential equation whose analytical solution is:
⎡
⎤
t
2
µ
µi2
1
i
τ dτ ⎦
θ̄ ji (t) = exp − t ⎣ f¯ji +S̄ ji (τ ) exp
Rj
Rj
Rj

(22b)

(23)

0

This expression gives a recursive rule [ S̄ ji (t) depends on θ̄ j−1 i (t)] for obtaining the transformed field for each species. Now, invoking Eqs. 9, 11, 20a, and 23, we obtain the explicit
analytical expression for the dimensionless concentration:c j (x, t) = F j (x) + exp( px + q j t)∞ψ̃i (x)θ̄ ji (t);j = 1, 2, 3, . . .

(24)

i=1

123

382
Table 1 Parameter values for thenitrification chain problem

J. S. Pérez Guerrero et al.
Description
Retardation coefficient(R j )
Decay constant, λ j(h−1 )
Pore velocity:
U = 1 cm h−1

NH+
4 ( j = 1)NO−
2 ( j = 2)NO−
3 ( j = 3)

2

1

1

0

0.1

0

Dispersion coefficient: D = 0.18 cm2 h−1

Fig. 1 Steady-state concentration for all the nitrogen species (L = 220 cm)

4 Test Case
The analytical solution presented above for the sequential decay problem was implemented
in a Mathematica (2007) program. Analytical expressions were developed for the transformed fields using Eqs. 22b and 23 for an arbitrary number of species. A generic analyticalexpression for the integral coefficients f¯ji was also developed. All analytical and functional
expressions were implemented for numerical evaluation. To calculate a specific set of eigenvalues βi , the transcendental Eq. 17 was solved numerically using the Mathematica (2007)
built-in function “FindRoot”.
The Mathematica program was first used to solve the three-species nitrification chain−−(NH+
4 → NO2 → NO3 ) transport problem as solved previously by van Genuchten (1985).
The parameter values for this problem are specified in Table 1. Assuming a transport domainof length of L = 220 cm, we calculated the steady-state concentration field for the three nitrogen species (Fig. 1) using Eqs. 5, 6 and 8. The strategy of Sun et al. (1999a,b) and Sun andClement (1999) as implemented here proved to be very efficient for obtaining the steady-statesolution.
Table 2 shows the convergence behavior of the NH+
4 concentration computed for timeT = 200 h, L = 220 cm, and Pe = 1, 222. Table 2 indicates that for this problem, N = 450
terms were needed to reach convergence throughout the spatial domain. Convergence isdefined here as having obtained the correct values up to 10 significant digits. Note thatwith N = 250, convergence was achieved for X ≤ 120 cm, but not beyond that location.
Convergence progressed slowly along the domain as N was increased from N = 250–450.
The region where convergence was slower (X > 120 cm) corresponded to those parts ofthe domain where the concentration of NH+
4 was very close to zero. Table 3 shows the123Analytical Solution for Multi-Species Contaminant Transport

383

Table 2 Solution convergence for the three-species nitrification transport problem assuming L = 220 cm
and T = 200 h (N = Number of terms summed)
X (cm)

0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75
80
85
90
95
100
105
110
115
120
125
130
135
140
145
150
155
160
165
170
175
180
185
190
195
200
205
210
215
220

Dimensionless concentration for NH+
4
N = 250

N = 350

N = 450

N = 500

0.9982064510
0.9496085026
0.9033765583
0.8593954286
0.8175555319
0.7777526219
0.7398875272
0.7038659047
0.6695980046
0.6369984464
0.6059860065
0.5764834154
0.5484171659
0.5217173284
0.4963172806
0.4721485541
0.4490140056
0.4250786668
0.3894312160
0.3149047564
0.1927162768
0.07678511830
0.01794434192
0.002312432594
0.0001586398313
5.675790084E-6
3.007974339E-7
0.001083916543−221818.1891−2.668137442E11−6.668710815E16
2.244041842E23
3.445925381E29
1.558852171E35−2.075877104E41−4.268009016E47−2.727350336E53
1.633178059E59
5.087387277E65
4.212543655E71−8.120405055E76−5.828191022E83−6.048183034E89−5.189939507E94
6.373088273E101

0.9982064510
0.9496085026
0.9033765583
0.8593954286
0.8175555319
0.7777526219
0.7398875272
0.7038659047
0.6695980046
0.6369984464
0.6059860065
0.5764834154
0.5484171659
0.5217173284
0.4963172806
0.4721485541
0.4490140056
0.4250786668
0.3894312160
0.3149047564
0.1927162768
0.07678511830
0.01794434192
0.002312432594
0.0001586398313
5.675789878E-6
1.045824992E-7
9.845112917E-10
4.709534978E-12
1.140537819E-14
1.394604866E-17
8.592942476E-21
2.663974946E-24
4.150608044E-28
3.247053381E-32
1.274515591E-36
2.513402272E-41
1.480383980E-37
2.595176680E-31
3.825694352E-25
5.163387350E-19
6.596017284E-13
8.107330268E-7
0.9676763723
1.127919833E6

0.9982064510
0.9496085026
0.9033765583
0.8593954286
0.8175555319
0.7777526219
0.7398875272
0.7038659047
0.6695980046
0.6369984464
0.6059860065
0.5764834154
0.5484171659
0.5217173284
0.4963172806
0.4721485541
0.4490140056
0.4250786668
0.3894312160
0.3149047564
0.1927162768
0.07678511830
0.01794434192
0.002312432594
0.0001586398313
5.675789878E-6
1.045824992E-7
9.845112917E-10
4.709534978E-12
1.140537819E-14
1.394604866E-17
8.592942476E-21
2.663974946E-24
4.150608044E-28
3.247053381E-32
1.274515591E-36
2.508537003E-41
2.474588289E-46
1.222969535E-51
3.026982697E-57
3.751109863E-63
2.326790851E-69
7.222868392E-76
1.121853644E-82
1.199389159E-89

0.9982064510
0.9496085026
0.9033765583
0.8593954286
0.8175555319
0.7777526219
0.7398875272
0.7038659047
0.6695980046
0.6369984464
0.6059860065
0.5764834154
0.5484171659
0.5217173284
0.4963172806
0.4721485541
0.4490140056
0.4250786668
0.3894312160
0.3149047564
0.1927162768
0.07678511830
0.01794434192
0.002312432594
0.0001586398313
5.675789878E-6
1.045824992E-7
9.845112917E-10
4.709534978E-12
1.140537819E-14
1.394604866E-17
8.592942476E-21
2.663974946E-24
4.150608044E-28
3.247053381E-32
1.274515591E-36
2.508537003E-41
2.474588289E-46
1.222969535E-51
3.026982697E-57
3.751109863E-63
2.326790851E-69
7.222868392E-76
1.121853644E-82
1.199389159E-89−−
converged values for the three-species nitrification chain (NH+4 → NO2 → NO3 ). With adomain length of L = 220 cm, the computed concentrations in this table for approximately
X < 200 cm were identical to those obtained previously by van Genuchten (1985) for a semiinfinite domain. Closer to the exit boundary, the computed concentrations diverged fromthose computed for the semi-infinite case.

123

384
Table 3 Converged
dimensionless concentration for
the three nitrogen species(L = 220 cm and T = 200 h)

J. S. Pérez Guerrero et al.
X (cm) Dimensionless concentration
NH+
4
0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75
80
85
90
95
100
105
110
115
120
125
130
135
140
145
150
155
160
165
170
175
180
185
190
195
200
205
210
215
220

NO−
2

NO−
3

0.9982064510
0.001731801827
0.00006174718691
0.9496085026
0.03871306202
0.01167843542
0.9033765583
0.05950592502
0.03711751672
0.8593954286
0.07048358153
0.07012098991
0.8175555319
0.07554095946
0.1069035086
0.7777526219
0.07705693685
0.1451904413
0.7398875272
0.07648303679
0.1836294360
0.7038659047
0.07470359556
0.2214304997
0.6695980046
0.07225611283
0.2581458825
0.6369984464
0.06946605852
0.2935354950
0.6059860065
0.06652934238
0.3274846512
0.5764834154
0.06356276291
0.3599538217
0.5484171659
0.06063486556
0.3909479685
0.5217173284
0.05778481534
0.4204978553
0.4963172806
0.05503393446
0.4486486814
0.4721485541
0.05239263717
0.4754531577
0.4490140056
0.04986128376
0.5009661232
0.4250786668
0.04738367837
0.5252204386
0.3894312160
0.04453022133
0.5480350476
0.3149047564
0.03978515558
0.5682353216
0.1927162768
0.03122025618
0.5826020944
0.07678511830
0.01994781008
0.5871394638
0.01794434192
0.01024425185
0.5808288008
0.002312432594 0.004491732945
0.5663503702
0.0001586398313 0.001831048663
0.5470470018
5.675789878E-6 0.0007320388051 0.5249781425
1.045824992E-7 0.0002918152583 0.5010561015
9.845112917E-10 0.0001162802778 0.4756178471
4.709534978E-12 0.00004632002377 0.4487583366
1.140537819E-14 0.00001844287785 0.4204734822
1.394604866E-17 7.337981865E-6
0.3907175422
8.592942476E-21 2.916377722E-6
0.3594259040
2.663974946E-24 1.157092755E-6
0.3265240692
4.150608044E-28 4.578688461E-7
0.2919313874
3.247053381E-32 1.804342084E-7
0.2555653614
1.274515591E-36 7.064376319E-8
0.2173648817
2.508537003E-41 2.737358199E-8
0.1774058936
2.474588289E-46 1.043032216E-8
0.1362795622
1.222969535E-51 3.865819869E-9
0.09580638645
3.026982697E-57 1.368811897E-9
0.05947262630
3.751109863E-63 4.505102185E-10 0.03133947460
2.326790851E-69 1.328924131E-10 0.01349921247
7.222868392E-76 3.370995078E-11 0.004600992689
1.121853644E-82 7.061392726E-12 0.001209079344
1.199389159E-89 1.255589051E-12 0.0002545665546

Table 4 illustrates for different domain lengths the divergence of the concentration at theexit boundary from the semi-infinite results. Table 4 shows converged dimensionless concentrations for X = L, T = 200 h, and domain lengths L = 100, L = 160, and L = 200 cm, as
well as results reported by van Genuchten (1985) for the same time T and position X but witha semi-infinite domain. For all three of these domain lengths, the computed concentrations
for positions X < (L − 5) cm were identical to the semi-infinite case (and results not shown).

123Analytical Solution for Multi-Species Contaminant Transport

385

Table 4 For time T = 200 h and various finite domain lengths L, comparison of dimensionless concentrations
at X = L with results obtained by van Genuchten (1985) for position X and a semi-infinite domain
X (cm)

100
160
200

L (cm)

100
∞
160
∞
200
∞

Table 5 Calculated
dimensionless concentration for
T = 50 h and L = 110 cm

Dimensionless concentration
NH+
4

NO−
2

NO−
3

0.1974489849
0.19272
3.284869976E-24
0.00000
5.005954795E-63
0.00000

0.03160207912
0.031221.194407808E-6
0.00001
4.686175169E-10
0.00000

0.5822557877
0.58260
0.3277406677
0.32652
0.03215551034
0.03134

X (cm) Dimensionless concentration
NH+
4
0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75
80
85
90
95
100
105
110

NO−
2

NO−
3

0.9982064510
0.001731801827
0.00006174718691
0.9496085026
0.03871306202
0.01167843542
0.9033763767
0.05950592248
0.03711751630
0.8590855526
0.07047834522
0.07011995330
0.7813793474
0.07471180447
0.1066758904
0.3980357655
0.06321867920
0.1389760794
0.03721508150
0.03218574711
0.1463644870
0.0003259844933 0.01199491648
0.1251952903
2.116560254E-7
0.004038071439
0.08898433433
9.362204927E-12 0.001163571008
0.04786979317
2.720387050E-17 0.0002378086225 0.01540368887
5.098177485E-24 0.00002569193818 0.002211598187
6.098329323E-32 1.106360548E-6
0.0001125983943
4.626179563E-41 1.571458652E-8
1.768548405E-6
2.216160379E-51 6.608912388E-11 7.926994338E-9
6.684354841E-63 7.749196672E-14 9.700019589E-12
1.266667903E-75 2.446642598E-17 3.155992504E-15
1.505598773E-89 2.036299547E-21 2.685192944E-19
1.121133343E-104 4.406821766E-26 5.908976694E-24
5.224910888E-121 2.457022876E-31 3.337684639E-29
1.522752711E-138 3.506493253E-37 4.812926628E-35
2.773488562E-157 1.274872527E-43 1.764637898E-41
5.1438380E-177
1.622019705E-50 2.260630087E-48

At the exit boundary, the computed concentrations diverged from the semi-infinite results(Table 4). Convergence for L = 160 cm was achieved with N = 250 terms in the series while
L = 220 cm required N = 450 terms. With all other parameters constant, the number theterms required for convergence was found to depend on the domain length. This is because
the domain length is contained inside of the Peclet number. Increasing L makes the effective
diffusion coefficient (1/Pe) in Eq. 13a is smaller, which necessitates the use of more terms
in the series to achieve convergence. Similar finding were obtained by van Genuchten andAlves (1982) for single-species transport in a finite domain.
For benchmark purposes, Table 5 represents the concentration fields for the three nitrogenspecies at a much smaller time (T = 50 h) and a domain L of 110 cm.

123

386

J. S. Pérez Guerrero et al.

5 Conclusion
Using the algebraic substitution given by Eq. 11 within the CITT, it was possible to obtainan analytical solution of a multi-species advection-dispersion transport equation featuring
sequential first-order decay reactions, distinct linear sorption reactions for each species, afinite domain, and steady-state boundary conditions. The convergence properties of the solution were examined using a classic test case involving a three-species nitrification chain.
Results showed that the solution for each species converged to previously reported literature
values. The presented solution procedure can be extended to decay chain problems subject
to transient boundary conditions. Because of the scarcity of available solutions for finite
domains, the solution procedure and the results presented in this study should be especially
useful for testing numerical solutions.
Acknowledgements The first author acknowledges support given by the Brazilian National Research Council, the University of California at Riverside (UCR) and the United States Salinity Laboratory of the USDA
Agricultural Research Service.

Appendix A
Consider the following advection-diffusion equation with constant coefficients R, A1 , A2 ,B1 , and source term S(x, t):
∂2 M
∂M
∂M
= A2 2 − A1
− A0 M + S(x, t)∂t∂x∂xUsing the algebraic substitution:R

M(x, t) = M(x, t) exp( px + qt)

(25)

(26)

Equation 25 can be reformulated as a diffusion equation in terms of M(x, t):∂2 M∂M= A2 2 + S(x, t) exp( px + qt)∂t∂x
where the coefficients p and q are given by:
R

p=

A1
1
; q=−
2 A2
R

A21
+ A0
4 A2

(27)

(28)

References
Aziz, C.E., Newell, C.J., Gonzales, J.R., Haas, P., Clement, T.P., Sun, Y.: BIOCHLOR—Natural attenuation
decision support system v1.0, user’s manual. U.S. Environmental Protection Agency Report EPA/600/R00/008. EPA Center for Subsurface Modeling Support (CSMOS), Ada, Oklahoma (2000)
Bauer, P., Attinger, S., Kinzelbach, W.: Transport of a decay chain in homogenous porous media: Analyticalsolutions. J. Contam. Hydrol. 49, 217–239 (2001). doi:10.1016/S0169-7722(00)00195-9
Cassol, M., Wortmann, S., Rizza, U.: Short communication: Analytic modeling of two-dimensional transient
atmospheric pollutant dispersion by double GITT and Laplace Transform techniques. Environ. Model.
Softw. 24(1), 144–151 (2009). doi:10.1016/j.envsoft.2008.06.001
Cho, C.M.: Convective transport of ammonium with nitrification in soil. Can. J. Soil Sci. 51(3), 339–350 (1971)
Clement, T.P.: Generalized solution to multispecies transport equations coupled with a first-order reaction
network. Water Resour. Res. 37, 157–163 (2001). doi:10.1029/2000WR900239
Cotta, R.M.: Integral Transforms in Computational Heat and Fluid Flow. CRC Press, Boca Raton, FL (1993)

123

Analytical Solution for Multi-Species Contaminant Transport

387

Lu, X., Sun, Y.: A solution of transport in porous media with equilibrium and kinetic reactions. Transp. PorousMedia 72(2), 199–206 (2008). doi:10.1007/s11242-007-9144-1
Lunn, M., Lunn, R.J., Mackay, R.: Determining analytic solutions of multiple species contaminant transport,with sorption and decay. J. Hydrol. (Amst.) 180, 195–210 (1996). doi:10.1016/0022-1694(95)02891-9
Mathematica: Wolfram Research, Inc., Champaign, IL (2007)
McLaren, A.D.: Nitrification in soil-systems approaching a steady state. Soil Sci Soc. Am. Proc. 33, 551–
556 (1969)
Mikhailov, M.D., Ozisik, M.N.: Unified Analysis and Solutions of Heat and Mass Diffusion. Wiley, New York(1984)
Ozisik, M.N.: Heat Conduction. Wiley, New York (1980)
Pérez Guerrero, J.S., Pimentel, L.C.G., Skaggs, T.H., van Genuchten, M.Th.: Analytical solution of theadvection-diffusion transport equation using a change-of-variable and integral transform technique. Int.
J. Heat Mass Transf. (2009). doi:10.1016/j.ijheatmasstransfer.2009.02.002
Pontedeiro, E.M., Simunek, J., van Genuchten, M.T., Cotta, R.M.: Performance Assessment Modeling of aRadioactive Mining Waste Disposal Site, Soil Science Society of America Annual Meeting. November
4–8, 2007, New Orleans, Louisiana, CDROM (2007)
Quezada, C.R., Clement, T.P., Lee, K.K.: Generalized solution to multi-dimensional multi-species transport
equations coupled with a first-order reaction network involving distinct retardation factors. Adv. Water
Resour. 27, 508–521 (2004)ºŠuimunek,J., Šejna, M., van Genuchten, M.T.: The HYDRUS Software Package for Simulating Two- andThree-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media. Technical Manual, Version 1.0, PC Progress, Prague, Czech Republic (2006)
Skaggs, T.H., Jarvis, N.J., Pontedeiro, E.M., van Genuchten, M.T., Cotta, R.M.: Analytical advection-dispersion model for transport and plant uptake of contaminants in the root zone. Vadose Zone J. 6, 890–898(2007). doi:10.2136/vzj2007.0124
Srinivasan, V., Clement, T.P.: Analytical solutions for sequentially coupled one-dimensional reactive transport
problems – Part I: Mathematical derivations. Adv. Water Resour. 31, 203–218 (2008a)
Srinivasan, V., Clement, T.P.: Analytical solutions for sequentially coupled one-dimensional reactive transport
problems – Part II: Special cases, implementation and testing. Adv. Water Resour. 31, 219–232 (2008b)
Sun, Y., Clement, T.P.: A decomposition method for solving coupled multi-species reactive transport problems. Transp. Porous Media 37, 327–346 (1999)
Sun, Y., Petersen, J.N., Clement, T.P.: Analytical solutions for multiple species reactive transport in multiple
dimensions. J. Contam. Hydrol. 35, 429–440 (1999a)
Sun, Y., Petersen, J.N., Clement, T.P., Skeen, R.S.: Development of analytical solutions for multi-species
transport with serial and parallel reactions. Water Resour. Res. 35(1), 185–190 (1999b)
Sun, Y., Buscheck, T.A., Mansoor, K., Lu, X.: Analytical Solutions for Sequentially Reactive Transport withDifferent Retardation Factors. UCRL-ID-145211, LLNL. http://www.doe.gov/bridge (2001)
van Genuchten, M.T.: Convective-dispersive transport of solutes involved in sequential first-order decay reactions. Comput. Geosci. 11(2), 129–147 (1985). doi:10.1016/0098-3004(85)90003-2
van Genuchten, M.T., Alves, W.J.: Analytical solutions of the one-dimensional convective-dispersive solutetransport equation. Techn. Bull. 1661, Agricultural Research Service, U.S. Department of Agriculture,
151 pp (1982)

123