Transp Porous Med
DOI 10.1007/s11242-010-9553-4

Analytical Solution for Multi-Species Contaminant
Transport in Finite Media with Time-Varying Boundary
Conditions
Jesús S. Pérez Guerrero · Todd H. Skaggs ·M. Th. van Genuchten

Received: 14 August 2009 / Accepted: 17 February 2010
© Springer Science+Business Media B.V. 2010

Abstract Most analytical solutions available for the equations governing the advective–
dispersive transport of multiple solutes undergoing sequential first-order decay reactions havebeen developed for infinite or semi-infinite spatial domains and steady-state boundary conditions. In this study, we present an analytical solution for a finite domain and a time-varying
boundary condition. The solution was found using the Classic Integral Transform Technique(CITT) in combination with a filter function having separable space and time dependencies, implementation of the superposition principle, and using an algebraic transformationthat changes the advection–dispersion equation for each species into a diffusion equation.
The analytical solution was evaluated using a test case from the literature involving a fourradionuclide decay chain. Results show that convergence is slower for advection-dominatedtransport problems. In all cases, the converged results were identical to those obtained withthe previous solution for a semi-infinite domain, except near the exit boundary where differences were expected. Among other applications, the new solution should be useful forbenchmarking numerical solutions because of the adoption of a finite spatial domain.
Keywords Multi-species transport · Finite domain · Analytical solution · Integral transform

Parts of the article were created within the capacity of a US governmental employment and therefore publicdomain. Published by Springer Science+Business Media B.V. 2010.
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

J. S. P. Guerrero et al.

List of Symbols
Ā jm ā jm
Auxiliary coefficients
B1 , B2
Mathematical operators
Auxiliary coefficientsB̄ jm b̄ jm
b jm
Boundary source coefficient
C0
Dimensional reference solute concentration
C j (X, T ) Dimensional solute concentration of the jth speciesc j (x, t)Dimensionless solute concentration of the jth species
D
Dispersion coefficientf¯jmi
Integral coefficient
G j (X )Dimensional initial concentration of the jth species
g j (x)Dimensionless initial concentration of the jth species
H1 , H2
Coefficientsh j (x, t)
Filter function
L
Mathematical operator
L0
Spatial length domain
Ni
Norm
p
Constant used in algebraic substitution
Pe
Peclet number
qj
Constant used in algebraic substitution
Rj
Retardation coefficient for the jth species
Integral coefficientS̄ jmi (t)TDimensional timetDimensionless timeU
Constant pore water velocity
w jm
Auxiliary coefficientr jm
Auxiliary coefficientX
Dimensional spatial coordinatex
Dimensionless spatial coordinate

Greek Symbols
αi
Source leaching coefficient
βi
EigenvalueMathematical operator
γjDamkohler number
δik
Kronecker delta
Auxiliary coefficientδ̄ jm
θ jm (x, t) Unknown function in purely diffusive equation
Integral transform of the function θ j (x, t)θ̄ jmi (t)λjFirst-order decay constant for the jth speciesλmSource decay constantµi
Eigenvalueξ
Auxiliary integration variable
σm
Auxiliary coefficient(T)
Dimensional transient boundary functionjφ j (t)Dimensionless transient boundary function

123Analytical Solution for Multi-Species Contaminant Transport

ϕ jm (x)ψi (x)
ψ̃i (x)ω j (x, t)ω j (x, t)
 jm (x, t)

Spatial function related to the filter
Eigenfunction
Normalized eigenfunction
Unknown function
Unknown function
Unknown function

1 Introduction
Analytical solutions for transport problems involving sequential decay reactions have beendeveloped mostly for steady-state boundary conditions and for infinite or semi-infinite spatial domains. Relatively very little literature is available about analytical solutions for multispecies transport problems for either finite media or time-dependent boundary conditions.
Detailed surveys of the literature were given recently by Srinivasan and Clement (2008a,b)
and Pérez Guerrero et al. (2009a). Hence, we mention here only a relatively few studies thatare most pertinent to our current efforts.
Higashi and Pigford (1980) previously used Laplace transforms to develop an analytic
solution for purely advective multi-species transport in finite media. Again using Laplacetransforms, van Genuchten (1985) later extended the analytical solution to advective–dispersive multi-species transport in a one-dimensional semi-infinite domain. He obtained solutions
for up to four species using either first- (Dirichlet) or third-type (Cauchy) inlet boundary conditions. The solutions for each species additionally featured transient source terms which wererelated through the Bateman equations (Bateman 1910).
Srinivasan and Clement (2008a,b) recently developed a closed-form analytical solution
for the sequential decay problem featuring an arbitrary number of species, spatially varying initial conditions, and a generic, exponentially decaying Bateman-type source boundary.
The coupled, one-dimensional system of partial differential equations was solved for both
Dirichlet and Cauchy boundary conditions using a combination of Laplace transforms andthe matrix diagonalization procedure of Clement (2001) which uncouples the transport equations. The solution was obtained for a semi-infinite domain.
Having a solution for a finite domain and time-varying boundary conditions would bevery useful for benchmarking numerical solutions since numerical solutions generally pertain only to finite domains. We recently used the Classic Integral Transform Technique (CITT)
to obtain an analytical solution for a first-order sequential decay problem with steady-stateboundary conditions and a finite domain (Pérez Guerrero et al. 2009a). The objective of thisstudy is to extend the CITT procedure to obtain an analytical solution for a sequential decay
reaction transport problem with time-varying boundary conditions and a finite domain.
2 Problem Formulation
Consider the transport of first-order sequentially decaying species in a homogeneous finite
porous media of length L 0 , subject to linear equilibrium adsorption processes and a constantadvective velocity. The governing transport equations are given by the following coupledsystem of dimensional partial differential equations:
Rj

∂C j
∂C j
∂ 2C j−U=D− λ j R j C j + λ j−1 R j−1 C j−1
2
∂T
∂X
∂X
( j = 1, 2, 3, . . . ; λ0 = 0)

(1a)

123

J. S. P. Guerrero et al.

where the index j signifies the particular solute species, R j is the retardation factor, λ j is thefirst-order decay coefficient [T−1 ], D is the dispersion coefficient [LT−2 ], U is the pore-watervelocity [LT−1 ], C j is the solute concentration [ML−3 ], X is position [L], and T is time [T].
The initial condition may be a function of position as follows:
C j (X, 0) = G j (X )(1b)

The boundary conditions are specified as being of the third type for the inlet boundary(X = 0)and the second type for the exit boundary (X = L 0 ):∂C j (0, T )(1c)+ U C j (0, T ) = U j (T )∂X
∂C j (L 0 , T )
=0
(1d)
∂X
The time-varying source functions j (T ) are the same as those used by van Genuchten(1985) for the generalized Bateman equations:−Dj (T ) =jb jm exp(−λm T )

(2a)

m=1

λm = λm + αm(2b)in which the coefficients b jm and λm account for decay reactions (λm ) in the waste source
and the leaching (αm ) of each species from the source.
Equations 1 and 2 can be written in dimensionless form after defining the following
variables:
Cj
X
T
x=
; t=
(3a–c); cj =
L0(L 0 /U )C0
g j (x) =Gj; φ j (t) =C0j (T )
C0(3d,e)

λ j R j L0
U L0
; γj =(3f,g)
D
U
where Pe is the finite domain Peclet number and γ j the Damkohler number for species j.
Equations 1 and 2 then become
Pe =

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 φ j (t)∂x∂c j (1, t)=0∂x
φ j (t) =j

f jm (t)

(4a)
(4b)
(4c)
(4d)
(5a)

m=1



L0
f jm (t) = b jm exp −λm t = b jm exp(−σm t)Uwith σm = λm LU0 .

123

(5b)

Analytical Solution for Multi-Species Contaminant Transport

Define the following operators on c j :
∂
cj
∂t
∂c j
1 ∂ 2c j−
Lc j ≡− γjcjPe ∂ x 2∂x∂c j
B1 c j ≡ −+ Pe c j∂x∂c j
B2 c j ≡∂x
c j ≡ R j

(6)
(7)
(8)
(9)

These operators allows Eq. 4 to be written more compactly asc j = Lc j + γ j−1 c j−1 ( j = 1, 2, 3, . . . ; γ0 = 0)(10a)c j (x, 0) = g j (x)(10b)

B1 c j (0, t) = Pej


f jm (t)

(10c)

m=1

B2 c j (1, t) = 0(10d)

3 Analytical Solution
The solution to Eq. 10 will be obtained using several steps. First, Eq. 10 will be transformed
into a system with homogeneous boundary conditions. Next, according to the procedurereported by Pérez Guerrero et al. (2009a,b), the advection–dispersion system with homogeneous boundary conditions will be transformed into a system of diffusion equations. The
final transformed system is subsequently solved analytically using the CITT (e.g., Ozisik
1980; Mikhailov and Ozisik 1984; Cotta 1993).
3.1 Homogenizing the Boundary Conditions
As noted by Ozisik (1980) and Cotta and Mikhailov (1997), solutions of nonhomogeneous
problems based on eigenfunction expansions may converge slowly or even exhibit anomalous
behavior, especially in the vicinity of the boundaries. The reason is that the source terms are
not accounted for in the auxiliary problems. Therefore, we homogenize the boundary conditions by defining a filter function, h j (x, t), conditioned to reproduce the original boundary
conditions given by Eq. 10c, d for each species:c j (x, t) = ω j (x, t) + h j (x, t) ( j = 1, 2, 3, . . .)(11)

where ω j (x, t)is an unknown function with homogeneous boundary conditions.
For the case of multi-species transport in a finite medium with non-time-dependent boundary conditions, Pérez Guerrero et al. (2009a) found an appropriate filter function h j using
the procedure presented by Sun et al. (1999a,b) and Sun and Clement (1999c). However,using the procedure for time-dependent boundary conditions is not easy because of the needto define the filter over both the time and spatial domains.
For this reason, i.e., for overcoming the above constraint, we implemented a version of themethod suggested by Ozisik (1980) for homogenizing problems with transient boundary conditions. The method involves defining a filter function with temporal and spatial dependencies

123

J. S. P. Guerrero et al.

that are separable. As detailed in Appendix A, the procedure leads to the following form forthe filter function when extended to multi-species transport problems:h j (x, t) =j


f jm (t)ϕ jm (x)

(12)

m=1

The temporally dependent term f jm (t) in this equation is known and given by Eq. 5b, whereas
Appendix A gives details for developing a mathematical formulation for ϕ jm (x). A general
analytical expression for ϕ jm (x), obtained by solving Eqs. A10–A12 with the Mathematica(2007) software, is:




ā jm
b̄ jm
ϕ jm (x) = Ā jm exp x+ B̄ jm exp x
2
2
 

√


x
b̄ jm
ā jm δ̄ jm γ j−1 r jm Peexp ξ+ exp x− Pe ϕ j−1m (ξ )dξ
2
w jm
2


b̄ jm
− exp x2δ̄ jm γ j−1 r jm
w jm

√

1

Pe



x
exp ξ

ā jm
− Pe
2


ϕ j−1m (ξ )dξ

1

(13a)
where w jm , ā jm , b̄ jm are respectively:
w jm =Pe − 4σm r j + 4γ j√ā jm = Pe − w jm Pe√b̄ jm = Pe + w jm Pe

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

The coefficients Ā jm and B̄ jm in Eq. 13a are calculated from the boundary conditions suchas specified in Appendix A. The integrals in Eq. 13a are analytic since the terms ϕ j−1m (ξ ) are
a combination of exponential functions. While it is evident that the above analytical solution
is relatively elaborate for large values of j, strategic grouping of the coefficients and using
the Mathematica software made it possible to generate simplified structures appropriate forcomputational purposes.
Substituting Eq. 11 into Eq. 10 yields a problem reformulated in terms of ω j (x, t):ω j = Lω j + γ j−1 ω j−1 ( j = 1, 2, 3, . . . ; γ0 = 0)
ω j (x, 0) = g j (x) − h(x, 0) = g j (x) −jf jm (0)ϕ jm (x)

(14a)
(14b)

m=1

B1 ω j (0, t) = 0(14c)

B2 ω j (1, t) = 0(14d)

In order to allow the use of a systematized solution procedure as described below, we decouple
the initial condition from the overall solution using the superposition principle. Thus,ω j (x, t) =

j

m=1

123

 jm (x, t)(15)

Analytical Solution for Multi-Species Contaminant Transport

Representing ω j (x, t) in terms of  jm (x, t) allows Eq. 14 to be rewritten as
 jm = L jm + γ j−1  j−1m δ̄ jm( j = 1, 2, 3, . . . ; m = 1, 2, . . . , j; γ0 = 0)(16a)

 jm (x, 0) = g j (x) − f jm (0)ϕ jm (x)(16b)

B1  jm (0, t) = 0(16c)

B2  jm (1, t) = 0(16d)

0 if j = m
1 if j  = m(16e)with the coefficient δ̄ jm defined byδ̄ jm =

In summary, the advection–dispersion problem with time-dependent boundary conditions
given by Eq. 10 was homogenized using a filter function having separable time and spacedependencies. The resulting problem (Eq. 14) with a coupled initial condition (Eq. 14b) was
re-written using the superposition principle so that the condition was decoupled (Eq. 16).
Once  jm (x, t) is found, c j (x, t) can be obtained from Eqs. 11, 12, and 15. The analyticalsolution for  jm (x, t) is obtained below.
3.2 Transforming the Advection–Dispersion Equation into a Diffusion Equation
As detailed in Pérez Guerrero et al. (2009a,b), the following algebraic substitution can beused to transform an advection–dispersion differential system such as Eq. 16 into a purely
diffusive system: jm (x, t) = θ jm (x, t) exp( px + q j t) ( j = 1, 2, 3, . . . ; m = 1, 2, . . . , j)

(17)

In which the coefficients p and q j are given by
Pe
p=
2

Pe
1
qj = −
+ γj
Rj
4

(18a)
(18b)

Using this substitution in Eq. 16 results in the purely diffusive system∂θ jm
1 ∂ 2 θ jm
+ γ j−1 θ j−1m (x, t) exp(− px − q j t)δ̄ jm=∂t
Pe ∂ x 2( j = 1, 2, 3, . . . ; m = 1, 2, . . . , j)

Rj

(19a)

with initial condition and boundary conditions as follows:θ jm (x, 0) = exp(− px)[g j (x) − f jm (0)ϕ j (x)]∂θ jm (0, t)
Pe+θ jm (0, t) = 0−∂x2∂θ jm (1, t)
Pe+θ jm (1, t) = 0
∂x
2

(19b)
(19c)
(19d)

3.3 Integral Transform Solution
Equation 19 can be solved using the CITT. In this procedure, the potential θ jm is expanded
in a series of eigenfunctions. An orthogonality property is established and used to specify

123

J. S. P. Guerrero et al.

a transform and inverse transform rule. The transform is applied to the system of partial
differential equations to remove the spatial variable, thus producing a system of ordinary
differential equations. The above system is solved, and the inverse transformation is applied.
Details of the systematized procedure used here can be found in Ozisik (1980), Mikhailov
and Ozisik (1984), and Cotta (1993). Specific steps in the solution procedure are as follows:(a) The auxiliary eigenvalue problem
A Sturm–Liouville eigenvalue problem is obtained applying a separation-of-variables to thehomogeneous version of Eq. 19b (i.e., without the source term γ j−1 θ j−1m (x, t) exp(− px −q j t)δ̄ jm ),µi2
d2 ψi (x)22(20a,b)
+
β
ψ
(x)
=
0;
β
=
i
i
i
dx 2
(1/Pe)dψi (0)
Pe
−
(20c)
+
ψi (0) = 0dx
2dψi (1)Pe(20d)+
ψi (1) = 0dx
2where ψ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)(21)

The eigenvalues and norm can be calculated from the following equations:
βi (H1 + H2 )βi2 − H1 H2
βi2 + H12 + H1 + H2
Ni =2tan(βi ) =

(22)
(23)

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

(24)

Ni

1
ψ̃i (x)ψ̃k (x)dx = δik

(25)

0

where δi,k is the Kronecker delta.
(b) The integral transform pair
A transform pair is established for the potential θ jm (x, t). The inverse transform is given
by a series expansion of θ jm (x, t) in terms of the eigenfunctions ψi , whereas the forward
transform is found by applying the orthogonality property of Eq. 25. The transform pair isgiven by1
θ̄ jmi (t) =ψ̃i (x)θ jm (x, t)dx (Transform)

(26)

0

θ jm (x, t) =

∞

i=1

123

ψ̃i (x)θ̄ jmi (t) (Inverse)(27)Analytical Solution for Multi-Species Contaminant Transport(c) Integral transform of the system of partial differential equations
1
Application of the operator 0 ψ̃i (x)dx to both sides of Eq. 19 removes the space variable
and leads to the following infinite system of ordinary differential equations:dθ̄ jmi (t)+ µi2 θ̄ jmi (t) = S̄ jmi (t)dtθ̄ jmi (t = 0) = f¯jmi(28b)S̄ jmi (t) = γ j−1 θ̄ j−1mi (t) exp[(q j−1 − q j )t]δ̄ jm

(28c)

Rj

(28a)

where

f¯jmi =

1ψ̃i (x) exp(− px)[g j (x) − f jm (0)ϕ j (x)]dx

(28d)

0

(d) Solution of the transformed system
Solving the system of ordinary differential equations given by Eq. 28, subject to the transformed initial condition, results in the transformed potential:

⎡

 ⎤
1
µi2
µi2
1
θ̄ jmi (t) = exp − t ⎣ f¯jmi +τ dτ ⎦(29)
S̄ jmi (τ ) exp
Rj
Rj
Rj
0

(e) Calculation of the concentration
The inverse formula (27) can be invoked to obtain the desired solution for the potential
θ jm (x, t). Once θ jm (x, t) is obtained, the desired concentration field for species j can becomputed as
c j (x, t) = h(x, t) +j


exp( px + qm t)

m=1

∞


ψ̃i (x)θ̄ jmi (t) ( j = 1, 2, 3, . . .)

(30)

i=1

For computational purposes, the infinite series in Eq. 30 is truncated to a number of terms N
which produces the desired precision.

4 Test Application
The classic four species radionuclide decay chain problem solved by van Genuchten (1985) isused here as a test case for the analytic solution procedure developed above. The test case is an
Table 1 Parameters values for the radionuclide decay chain problem
Parameter

238 Pu( j = 1)234 U
( j = 2)230 Th( j = 3)

226 Ra( j = 4)Retardation coefficient, R j

10000

14000

50000

500

Decay constant, λ j [year−1 ]0.00790.0000028

0.0000087

0.00043

Source decay constant, λ j [year−1 ]0.00890.00100280

0.00100870

0.00143

Pore-water velocity: U = 100 m/year

Dispersion coefficient:
D = 10 m2 /year123

J. S. P. Guerrero et al.
Table 2 Values for the boundary source coefficient b jm
Species, j

b jm



m=1

m=2

m=3

m=4

1.25

****

****

****

−1.25044

1.25044

****

****

0.443684E−3

0.593431−0.593874****−0.516740E−60.120853E−1−0.122637E−1

0.178925E−3



j = 1 238 Pu
j = 2 234 U
j = 3 230 Th
j = 4 226 Ra

Table 3 Convergence of dimensionless concentration for 238 Pu(D = 10 m2 /year)
X (m)Numbers of terms summed (N)N =0

N = 150

N = 200

N = 300

02.812617661E−39

2.812617661E−39

2.812617661E−39

2.812617661E−39

54.660944797E−39

4.660944797E−39

4.660944797E−39

4.660944797E−39

107.723910258E−39

7.723910258E−39

7.723910258E−39

7.723910258E−39

151.279972029E−38

1.279972029E−38

1.279972029E−38

1.279972029E−38

20

2.121112674E−38

2.121112674E−38

2.121112674E−38

2.121112674E−38

25

3.515013512E−38

3.515013512E−38

3.515013512E−38

3.515013512E−3830

5.824923938E−38

5.824923938E−38

5.824923938E−38

5.824923938E−38

35

9.652804682E−38

9.652804682E−38

9.652804682E−38

9.652804682E−38

40

1.599619827E−37

1.599619827E−37

1.599619827E−37

1.599619827E−37452.650818778E−37

2.650818778E−37

2.650818778E−37

2.650818778E−3750

4.392818890E−37

4.392818890E−37

4.392818890E−37

4.392818890E−3755

7.279583939E−37

7.279583939E−37

7.279583939E−37

7.279583939E−3760

1.206340249E−36

1.206340249E−36

1.206340249E−36

1.206340249E−36

65

1.999093366E−36

1.999093366E−36

1.999093366E−36

1.999093366E−36

70

3.312808546E−36

3.312808545E−36

3.312808545E−36

3.312808545E−3675

5.489838868E−36

5.489838122E−36

5.489838122E−36

5.489838122E−3680

9.097516619E−36

9.097255384E−36

9.097255384E−36

9.097255384E−3685

1.507599961E−35

1.504826035E−35

1.504825228E−35

1.504825228E−35

90

2.498327551E−35−5.724418127E−31

2.405667399E−35

2.405667399E−35

95

4.140117217E−35

3.054489808E−20

3.094978514E−35

3.094978514E−35

100

6.860817975E−35−8.958409266E−10

2.234612976E−35

2.234612976E−35

105

1.136944217E−34−44.76951213

6.612991375E−36

6.612991375E−361101.864865906E−34

1.027415565E13

7.119825525E−37

7.119825525E−37

extended version of the problem studied previously by Higashi and Pigford (1980) in which
the concentration field was obtained for the decays chain 238 Pu → 234 U → 230 Th →226 Ra. The governing system of transport equations is given by Eq. 1 with C = 238 Pu,
1
C2 = 234 U, C3 = 230 Th, and C4 = 226 Ra. Model parameter values for the example application are specified in Tables 1 and 2 and are the same as those employed previously by van
Genuchten (1985).

123

Analytical Solution for Multi-Species Contaminant Transport
Table 4 Convergence of dimensionless concentration for 234 U (D = 10 m2 /year)
X (m)Numbers of terms summed (N)N =0

N = 250

N = 300

N = 400

00.00005599763985

0.00005599763985

0.00005599763985

0.00005599763985

50.0001139083231

0.0001139083231

0.0001139083231

0.0001139083231

100.0002317080881

0.0002317080881

0.0002317080881

0.0002317080881

150.0004713320030

0.0004713320030

0.0004713320030

0.0004713320030

200.0009587660877

0.0009587660877

0.0009587660877

0.0009587660877

25

0.001950286433

0.001950286433

0.001950286433

0.001950286433

30

0.003967200362

0.003967200362

0.003967200362

0.003967200362

35

0.008069931907

0.008069931907

0.008069931907

0.008069931907

400.01641555632

0.01641555632

0.01641555632

0.01641555632

450.03339191609

0.03339191609

0.03339191609

0.03339191609

50

0.06792455098

0.06792456385

0.06792456385

0.06792456385

550.1381008536

0.1381533789

0.1381533789

0.1381533789

60

0.1800304721

0.2788444392

0.2788444392

0.2788444392

65−147.6865682

0.4978363297

0.4978363297

0.4978363297
0.5054199873

70−217564.8353

0.5054199873

0.5054199873

75−3.192736368E8

0.1692657773

0.1692657773

0.1692657773

80−4.685275158E11

0.01366328526

0.01366328526

0.01366328526

85−6.875545178E14

0.0002334133628

0.0002334133628

0.0002334133628

90−1.008972150E18

8.083817021E−7

8.083817021E−7

8.083817021E−7

95−1.480645932E21

5.727475188E−10

5.727475188E−10

5.727475188E−10100−2.172817530E24

8.920290957E−14

8.920290957E−14

8.920290957E−14

105−3.188565151E27−1.753009868E−12

3.736568695E−18

3.736568695E−18

110−3.880320834E30

0.1236300738

7.792826947E−23

7.792826947E−23

The analytical solution of this study was implemented using the symbolic-numerical software package Mathematica 6.0 (Wolfram Research, Inc 2007), with the Mathematica systemvariable $MaxExtraPrecision set to 300. Computations were done on a desktop PC with a
2.0 GHz processor and 1 GB of RAM. The concentration field for each radionuclide wascomputed at time t = 10000 years assuming a domain length L 0 of 110 m. Tables 3, 4, 5,
and 6 show, with 10 significant figures, the convergence of the dimensionless concentrationof the four species for differing numbers of terms N in the truncated infinite summation of
Eq. 30. CPU time duration required to compute the results ranged from 5 s (Table 3) to 75 s
(Table 6). The case N = 0 corresponds the value of the filter function h j (x, t), as given in
Eq. 12.
The results show that the number of terms required to achieve convergence throughout the
spatial domain differed depending upon the species. Convergence of the 238 Pu concentration
required N = 200 terms, 234 U required N = 300 terms, and 230 Th and 226 Ra both required
N = 600 terms. The Peclet number for those calculations was relatively large (1100), reflecting advection-dominated transport. The large Peclet number was responsible for the large
number of terms required for the summation to converge. The converged dimensionless concentrations reported in Tables 3, 4, 5, and 6 for each radionuclide are in exact agreement

123

J. S. P. Guerrero et al.
Table 5 Convergence of dimensionless concentration for 230 Th (D = 10 m2 /year)
X (m)Numbers of terms summed (N)N =0

N = 500

N = 600

N = 7000

1.557998045E−6

1.557998045E−6

1.557998045E−6

1.557998045E−6

50.00001620869188

0.00001620869188

0.00001620869188

0.00001620869188

100.0001410791354

0.0001410795680

0.0001410795680

0.0001410795680

150.0007747723527

0.0007919552487

0.0007919552487

0.0007919552487

20

93.59969043

0.001298528478

0.001298528478

0.001298528478

25−6.857605987E12

0.001259505239

0.001259505239

0.001259505239

30

5.002484732E23

0.001268644491

0.001268644491

0.001268644491

35−3.634030382E34

0.001271527472

0.001271527472

0.001271527472

40

2.629133751E45

0.001269732017

0.001269732017

0.001269732017

45−1.894414442E56

0.001258502707

0.001258502707

0.001258502707

50

1.359496019E67

0.001228041189

0.001228041189

0.001228041189

55−9.716494897E77

0.001158416815

0.001158416815

0.001158416815

60

6.915816389E88

0.001009486912

0.001009486912

0.001009486912

65−4.901592658E99

0.0007163481928

0.0007163481928

0.0007163481928

70

3.458849405E110

0.0003080549531

0.0003080549531

0.0003080549531

75−2.429675808E121

0.00005314325832

0.00005314325832

0.00005314325832

80

1.698582612E132

2.636721747E−6

2.636721747E−6

2.636721747E−6

85−1.181462475E143

3.120865554E−8

3.120865554E−8

3.120865554E−890

8.173153485E153

7.933489062E−11

8.108033619E−11

8.108033619E−1195−5.620800411E164

0.1448010034

4.563235436E−14

4.563235436E−14
5.930826803E−18100

3.840591902E175−4.795194372E9

5.930826803E−18

105−2.605422444E186−2.882208268E20

2.193954540E−22

2.193954540E−221101.753240837E197

5.211513301E31

4.335249640E−27

4.335249640E−27

Table 6 Convergence of dimensionless concentration for 226 Ra (D = 10 m2 /year)
X (m)Numbers of terms summed (N)N =0

N = 500

N = 600

N = 700

02.397063545E−8

2.397063545E−8

2.397063545E−8

2.397063545E−8

51.691266186E−7

1.691266186E−7

1.691266186E−7

1.691266186E−7

101.511139130E−6

1.511140077E−6

1.511140077E−6

1.511140077E−6

150.00001052414499

0.00001057726603

0.00001057726603

0.00001057726603

20

0.09247382809

0.00003590353701

0.00003590353701

0.00003590353701

25−6.774597698E9

0.00006313026700

0.00006313026700

0.00006313026700

30

4.941931866E20

0.00009010807819

0.00009010807819

0.00009010807819

35−3.590042051E31

0.0001169082851

0.0001169082851

0.0001169082851

40

2.597309249E42

0.0001434364490

0.0001434364490

0.0001434364490

45−1.871483393E53

0.0001695417278

0.0001695417278

0.0001695417278

50

1.343039920E64

0.0001949155053

0.0001949155053

0.0001949155053

123

Analytical Solution for Multi-Species Contaminant Transport
Table 6 continued
X (m)Numbers of terms summed (N)N =0

N = 500

N = 600

N = 700
0.0002189277052

55−9.598881022E74

0.0002189277052

0.0002189277052

60

6.832103489E85

0.0002402964566

0.0002402964566

0.0002402964566

65−4.842261046E96

0.0002565305838

0.0002565305838

0.0002565305838

70

3.416981562E107

0.0002643627198

0.0002643627198

0.0002643627198

75−2.400265656E118

0.0002642455651

0.0002642455651

0.0002642455651

80

1.678022020E129

0.0002610955418

0.0002610955418

0.0002610955418

85−1.167161394E140

0.0002576059772

0.0002576059772

0.0002576059772

90

8.074221075E150

0.0002541482106

0.0002541482106

0.0002541482106

95−5.552763107E161

0.0002785185520

0.0002507349829

0.0002507349829

100

3.794103235E172−919886.1135

0.0002473658487

0.0002473658487

105−2.573884957E183−5.531608337E16

0.0002440402537

0.0002440402537

110

1.732018632E194

1.000039572E28

0.0002408228593

0.0002408228593

Table 7 Comparison of calculated dimensionless concentrations (D = 20m2 /year)
X (m)238 PuPresent solution (N = 100)

234 U

CHAIN

Present solution (N = 150)

CHAIN

0

2.842233909E−39

2.842E−39

5.684118059E−05

5.684E−05

54.735141085E−39

4.735E−39

1.168676484E−04

1.169E−04

107.888710715E−39

7.889E−39

2.402843694E−04

2.403E−04

151.314253486E−38

1.314E−38

4.940338836E−04

4.940E−0420

2.189536778E−38

2.190E−38

1.015752622E−03

1.016E−0325

3.647752397E−38

3.648E−38

2.088426366E−03

2.088E−0330

6.077129045E−38

6.077E−38

4.293884745E−03

4.294E−0335

1.012445293E−37

1.012E−37

8.828391788E−03

8.828E−0340

1.686726517E−37

1.687E−37

1.815150838E−02

1.815E−02452.810074148E−37

2.810E−37

3.731952687E−02

3.732E−0250

4.681563156E−37

4.682E−37

7.667966749E−02

7.668E−0255

7.799450274E−37

7.799E−37

1.559649552E−01

1.560E−0160

1.299382752E−36

1.299E−36

2.956393295E−01

2.956E−0165

2.164761203E−36

2.165E−36

4.425247909E−01

4.425E−0170

3.606405149E−36

3.606E−36

4.204066087E−01

4.204E−0175

6.005587956E−36

6.006E−36

2.139721046E−01

2.140E−0180

9.951251375E−36

9.951E−36

5.284934788E−02

5.285E−02851.597732605E−35

1.598E−35

6.010275612E−03

6.011E−03902.293955203E−35

2.294E−35

3.061362252E−04

3.061E−04

95

2.582403850E−35

2.582E−35

6.885135661E−06

6.885E−06

100

1.999021368E−35

1.999E−35

6.795238168E−08

6.795E−081059.680618782E−36

9.681E−36

2.942584009E−10

2.943E−101102.929004446E−36

2.767E−36

6.811823364E−13

5.625E−13

123

J. S. P. Guerrero et al.
Table 7 continued
X (m)230 ThPresent solution (N = 300)

226 RaCHAIN

Present solution (N = 300)

CHAIN

0

1.651367539E−6

1.651E−06

2.486450976E−8

2.486E−08

51.922822867E−5

1.923E−05

2.008439226E−7

2.008E−07

101.714762193E−4

1.715E−04

1.913007451E−6

1.913E−06

157.961635560E−4

7.961E−04

1.208412229E−5

1.208E−05

20

1.268186445E−3

1.268E−03

3.638133476E−5

3.638E−0525

1.260841900E−3

1.261E−03

6.366748116E−5

6.366E−0530

1.268188125E−3

1.268E−03

9.062266913E−5

9.062E−05351.270842955E−3

1.271E−03

1.174019006E−4

1.174E−0440

1.267980669E−3

1.268E−03

1.438914086E−4

1.439E−04451.254346371E−3

1.254E−03

1.699158765E−4

1.699E−0450

1.218541546E−3

1.219E−03

1.951154845E−4

1.951E−04
2.187E−04

55

1.137533012E−3

1.138E−03

2.187512277E−460

9.709188295E−4

9.709E−04

2.393662913E−4

2.394E−04

65

6.855210651E−4

6.855E−04

2.545932531E−4

2.546E−04703.449596549E−4

3.450E−04

2.623426781E−4

2.623E−0475

1.058832537E−4

1.059E−04

2.633518991E−4

2.633E−04801.754331684E−5

1.754E−05

2.609030583E−4

2.609E−04851.446926314E−6

1.446E−06

2.575425035E−4

2.575E−04905.652603367E−8

5.659E−08

2.540953389E−4

2.541E−04

95

1.015680336E−9

1.016E−09

2.506839621E−4

2.507E−04100

8.258661987E−12

8.259E−12

2.473163668E−4

2.473E−04105

3.018191040E−14

3.018E−14

2.439922763E−4

2.440E−04

110

6.031571078E−17

4.966E−17

2.408414731E−4

2.407E−04

with the results reported by van Genuchten (1985) for a semi-infinite domain, except near
the boundary at X = 110 m where, as expected, the results deviated slightly due to the effectof the finite domain. Optimal values of N for achieving convergence to a desired level ofaccuracy can be found using an adaptive procedure described by Cotta (1993). In addition,
it may be possible to improve the rate of convergence by finding better auxiliary eigenvalue
problems and filters (Cotta 1993).
Table 7 gives results for a less advective case. All model parameters are the same as before
except that the dispersion parameter was increased to D = 20 m2 /year. The Peclet numberin this case was Pe = 550, which reduced the number of terms required for convergence byhalf. The required CPU time was 18 s. Also shown in Table 7 are results obtained with the
CHAIN program of van Genuchten (1985) that is contained within the STANMOD software
package (Simunek et al. 1999). CHAIN implements the original analytical solutions of van
Genuchten (1985). The results obtained with the solution of this study were identical to thoseobtained with CHAIN except near the exit boundary, where again slight differences appear
due to the differing spatial domains (finite and infinite). Figure 1 contains plots of the data
presented in Table 7.

123Analytical Solution for Multi-Species Contaminant Transport
3238

Pu (x 1035)

Present Solution
CHAIN
2

1

0
5

U (x 10)

234

3
2
1
0

16

Th (x 104)

12

8

230

Dimensionless Concentration

4

4

0

226

Ra (x 104)

3

2

1

0
0

40

80

120

Distance (m)Fig. 1 Computed results for the decay chain 238 Pu → 234 U → 230 Th → 226 Ra, as given in Table 7(t = 10000 year, D = 20 m2 /year, U = 100 m/year)123

J. S. P. Guerrero et al.

5 Conclusions
Using the CITT in combination with a filter function with separable space and time-dependencies, the superposition principle, and a classic algebraic substitution, we were able to obtainan analytical solution for the transport of multiple solutes undergoing sequential first-orderdecay reactions for a finite domain and time-varying inlet condition. The analytical solution
is general and permits different values for the retardation coefficients of each species.
The convergence of the solution was investigated by means of a test case in which transport was advection-dominated. Those conditions typically require a large number of terms
for the eigen-expansion series to converge. Another test case with more dispersion required
fewer terms. In all the cases, the converged results were exactly the same as those computed
using analytical solutions published previously by van Genuchten (1985) for a semi-infinite
domain, except near the exit boundary where differences were expected. Among other applications, the new solution should be useful for benchmarking numerical solutions involving
the finite boundary conditions.
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. The research was supported also in part by a fellowship from CAPES to M.
Th. van Genuchten.

Appendix A
Obtaining a solution of Eq. 10 requires a filter function, h j (x, t), capable of homogenizing
the boundary conditions (Eqs. 10c, 10d). Additionally, it is desirable that the filter function
be as similar as possible to c j (x, t). In order to accomplish these objectives, we consider asystem for h j with the same structure as Eq. 10,
h j = Lh j + γ j−1 h j−1
B1 h j (0, t) = Pe

j


(A1)

f jm (t)

(A2)

m=1

B2 h j (1, t) = 0(A3)

but impose the following functional form on the solution:h j (x, t) =j


f jm (t)ϕ jm (x)

(A4)

m=1

This representation of h j , in which the time and spatial dependences are separated into
different terms in the series, is an extension of the split procedure presented by Ozisik (1980).
In Eq. A4, the function associated with the time domain ( f jm ) is known and given by Eq. 5b,
while ϕ jm (x) is unknown and needs to be found. Substituting Eq. A4 into Eq. A1 gives:
j

m=1

123

 f jm (t)ϕ jm (x) =

j

m=1

f jm (t)Lϕ jm (x) + γ j−1

j−1

m=1

f j−1m (t)ϕ j−1m (x)(A5)

Analytical Solution for Multi-Species Contaminant Transport

Rearranging this result gives:jR j σm f jm (t)ϕ jm (x) + f jm (t)Lϕ jm (x) + δ̄ jm γ j−1 f j−1m (t)ϕ j−1m (x) = 0 (A6)m=1

For the boundary condition we similarly obtain:j

f jm (t)B1 ϕ jm (0) − Pe f jm (t) = 0

(A7)

m=1
j


f jm (t)B2 ϕ jm (1) = 0(A8)

f j−1m (t)b j−1m=f jm (t)
b jm

(A9)

m=1

r jm =

Therefore, the equation and boundary conditions which determinate ϕ jm (x) are given by
R j σm ϕ jm (x) + Lϕ jm (x) + δ̄ jm γ j−1 r jm ϕ j−1m (x) = 0( j = 1, 2, . . . ; m = 1, 2, . . . , j)

(A10)

B1 ϕ jm (0) − Pe = 0(A11)

B2 ϕ jm (1) = 0(A12)A solution for ϕ jm (x) is given in the main text.

References
Bateman, H.: The solution of a system of differential equations occurring in the theory of radioactive transformations. Proc. Camb. Philos. Soc. 15, 423–427 (1910)
Clement, T.P.: Generalized solution to multispecies transport equations coupled with a first-order reaction
network. Water Resour. Res. 37, 157–163 (2001)
Cotta, R.M.: Integral Transforms in Computational Heat and Fluid Flow. CRC Press, Boca Raton (1993)
Cotta, R.M., Mikhailov, M.D.: Heat Conduction: Lumped Analysis, Integral Transforms, Symbolic Computations. Wiley-Interscience, New York (1997)
Higashi, K., Pigford, T.: Analytical models for migration of radionuclides in geologic sorbing media. J. Nucl.
Sci. Technol. 17, 700–709 (1980)
Mikhailov, M.D., Ozisik, M.N.: Unified Analysis and Solutions of Heat and Mass Diffusion. Wiley, NewYork (1984)
Ozisik, M.N.: Heat Conduction. Wiley, New York (1980)
Pérez Guerrero, J.S, Skaggs, T.H., van Genuchten, M.Th.: Analytical solution for multi-species contaminanttransport subject to sequential first-order decay reactions in finite media. Transp. Porous Media (2009a).
doi:10.1007/s11242-009-9368-3
Pérez Guerrero, J.S., Pimentel, L.C.G., Skaggs, T.H., van Genuchten, M.Th.: Analytical solution of advection–diffusion transport equation using change-of-variable and integral transform. Int. J. Heat Mass
Transf. 52, 3297–3304 (2009b)
Simunek, J., van Genuchten, M. Th., Sejna, M., Toride, N., Leij, F.J. The STANMOD Computer Softwarefor Evaluating Solute Transport in Porous Media Using Analytical Solutions of Convection–Dispersion Equation, Versions 1.0 and 2.0, IGWMC-TPS-71. International Ground Water Modeling Center,
Colorado School of Mines, Golden, 32 pp (1999)
Srinivasan, V., Clement, T.P.: Analytical solutions for sequentially coupled one-dimensional reactive transport
problems—Part I: Mathematical derivations. Adv. Water Res. 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 Res. 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)

123

J. S. P. Guerrero et al.
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)
van Genuchten, M.Th.: Convective–dispersive transport of solutes involved in sequential first-order decay
reactions. Comput. Geosci. 11(2), 129–147 (1985)
Wolfram Research, Inc. Mathematica, Version 6.0, Champaign (2007)123