International Journal of Heat and Mass Transfer 52 (2009) 3297–3304

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer
journal homepage: www.elsevier.com/locate/ijhmt

Analytical solution of the advection–diffusion transport equation using
a change-of-variable and integral transform technique
J.S. Pérez Guerrero a, L.C.G. Pimentel b, T.H. Skaggs c,*, M.Th. van Genuchten daRadioactive Waste Division, Brazilian Nuclear Energy Commission, DIREJ/DRS/CNEN, R. General Severiano 90, 22290-901 RJ-Rio de Janeiro, Brazil
Department of Meteorology, Federal University of Rio de Janeiro, Rio de Janeiro, Brazilc
U.S. Salinity Laboratory, USDA-ARS, 450 W. Big Springs Rd, Riverside, CA 92507, USAd
Department of Mechanical Engineering, LTTC/COPPE, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, Brazilb

a r t i c l e

i n f oArticle history:
Received 21 July 2008
Received in revised form 13 January 2009
Available online 11 March 2009
Keywords:
Analytical solution
Transport equation
Integral transforms

a b s t r a c t
This paper presents a formal exact solution of the linear advection–diffusion transport equation with constant coefﬁcients for both transient and steady-state regimes. A classical mathematical substitution
transforms the original advection–diffusion equation into an exclusively diffusive equation. The new diffusive problem is solved analytically using the classic version of Generalized Integral Transform Technique (GITT), resulting in an explicit formal solution. The new solution is shown to converge fasterthan a hybrid analytical–numerical solution previously obtained by applying the GITT directly to theadvection–diffusion transport equation.
Published by Elsevier Ltd.

1. Introduction
Analytical solutions of advective–diffusive transport problems
continue to be of interest in many areas of science and engineering,such as heat and mass transfer and pollutant dispersion in air, soils,
and water. They are useful for a variety of applications [1–5], suchas providing initial or approximate analyses of alternative pollution scenarios, conducting sensitivity analyses to investigate theeffects of various parameters or processes on contaminant transport, extrapolation over large times and distances where numericalsolutions may be impractical, serving as screening models or
benchmark solutions for more complex transport processes thatcannot be solved exactly, and for validating more comprehensive
numerical solutions of the governing transport equations.
The literature presents several methods to analytically solve thepartial differential equations governing transport phenomena[6–10]. For example, the method of separation-of-variables is oneof the oldest and most widely used techniques. Similarly, the classical Green’s function method can be applied to problems withsource terms and inhomogeneous boundary conditions on ﬁnite,semi-inﬁnite, and inﬁnite regions [10,11].
Integral transform techniques, such as the Laplace and Fouriertransform methods, employ a mathematical operator that produces a new function by integrating the product of an existing
function and a kernel function between suitable limits. The kernelof an integral transform, along with the integration limits, distin* Corresponding author. Tel.: +1 951 369 4853; fax: +1 951 342 4964.
E-mail address: Todd.Skaggs@ars.usda.gov (T.H. Skaggs).
0017-9310/$ - see front matter Published by Elsevier Ltd.
doi:10.1016/j.ijheatmasstransfer.2009.02.002

guishes one integral transform from another. Exact solutions of linear diffusion problems by classical integral transform techniques
were reviewed and classiﬁed by Mikhailov and Ozisik [12]. They
identiﬁed and uniﬁed seven classes of problems and demonstrated
many applications in heat and mass diffusion. Cotta [13] generalized and extended the classical integral transform method presented by Mikhailov and Ozisik [12], thereby creating a new
systematic procedure referred to as the Generalized Integral Transform Technique (GITT).
The literature also features the use of mathematical substitutions in obtaining analytical solutions to partial differential equations. Mathematical substitutions can simplify the structure of an
equation, thereby facilitating more ﬂexible applications of certain
solution methods. Compilations of transformations and substitutions are presented by Zwillinger [14] and Polyanin [15]. Most
existing analytical solutions for advection–diffusion transport
problems [3,16–18], including problems with growth and decay
terms, are for semi-inﬁnite or inﬁnite regions, with solutions forﬁnite domains being mostly limited to one-dimensional problems.
The aim of this paper is to present an analytical methodology tosolve advection–diffusion transport problems in a ﬁnite domain forboth transient and steady-state regimes. The proposed methodology uses change-of-variables in combination with the classic version of the Generalized Integral Transform Technique (GITT).
2. Problem formulation
We study transport in a ﬁnite domain and consider a threedimensional linear problem with decay and source terms. The

3298

J.S. Pérez Guerrero et al. / International Journal of Heat and Mass Transfer 52 (2009) 3297–3304

Nomenclatureu, v, wvVx, y, zageneric position on the boundaryCdimensionless solute concentration
Dx, Dy, Dz dispersion–diffusion constantsFﬁlter functionfintegral transform of qiGsource term
G
integral transform of G
coefﬁcients
H1, H2i, jindicesLmathematical operator
domain length
L0
M
function with homogeneous boundary conditions
norm
Ni
Pe
Peclet number
p1, p2, p3 constants for the algebraic substitution
q1, q2, q3 constants for the algebraic substitution
R
coefﬁcientSsource termTscalar quantityttime

advection (convection) components
symbolic coordinates for integration
generic ﬁnite volume
spatial coordinates

Greek symbols
eigenvaluebicdimensionless coefﬁcient
Kronecker deltadijhfunction in purely diffusive equation
hi
integral transform of the function hkgeneric decay constant (Eq. (2))
decay constant for the second test casek*lieigenvaluepnumber Pi
qgeneric initial conditions
auxiliary variabler2
Laplacian operator
wieigenfunction~inormalized eigenfucntion
wboundary conditions can be any combination of the ﬁrst, second, or
third kinds. We begin with formal mathematical deﬁnitions forboth transient and steady-state problems.

as systematized by Cotta [13]. In this technique, termed the Generalized Integral Transform Technique (GITT), the unknown function
is represented in terms of an eigenfunction series expansion. Basically, the GITT has the following steps:

2.1. Transient problem(a) Choose an appropriate auxiliary eigenvalue problem andﬁnd the associated eigenvalues, eigenfunctions, norm, andorthogonalization property.
(b) Develop the integral and inverse transforms.
(c) Transform the partial differential equation into a system ofordinary differential or algebraic equations.
(d) Solve the ordinary differential or algebraic system.
(e) Use the inverse transform to obtain the unknown function.

The transient transport equation for a scalar quantity T(x, y, z, t)
undergoing constant advection (or convection) and dispersion (ordiffusion) is given by:R@Tðx; y; z; tÞþ LTðx; y; z; tÞ þ Sðx; y; z; tÞ ¼ r2 Tðx; y; z; tÞ@tð1Þwhere the operators L and r2 are deﬁned by:

Lu

@
@
@
þv
þ w þ kR@x@y@zr2  D x@2@2@2
þ Dy 2 þ Dz 2
2
@x
@y
@z

ð2Þ
3.2. Transient regime

ð3Þ

in which the parameters u, v, and w are constant velocity coefﬁcients, k is a constant decay term, Dx, Dy, and Dz are constant dispersion (or diffusion) coefﬁcients, and S(x, y, z, t) is a source term. The
coefﬁcient R in Eq. (1) is also a constant parameter (e.g., the retardation factor in many subsurface contaminant transport problems).
A general initial condition can be established as:

Tðx; y; z; 0Þ ¼ qðx; y; zÞConsidering a more general case we assume that the boundary
conditions associated with Eq. (1) are non-homogeneous. To
homogenize the boundary conditions, we rewrite Eq. (1) using
the following ﬁlter strategy:Tðx; y; z; tÞ ¼ Mðx; y; z; tÞ þ Fðx; y; z; tÞwhere F(x, y, z; t) is any function satisfying exactly the original
boundary conditions of T(x, y, z, t). Eq. (1) now becomes:

ð4Þ
R

2.2. Steady-state problem
The transport equation for a steady-state transport regime independent of time, t, follows from Eq. (1) as:

r2 Tðx; y; zÞ ¼ LTðx; y; zÞ þ Sðx; y; zÞð5Þ

3. Mathematical analysis
3.1. Generalized Integral Transform Technique (GITT)
The transient and steady-state problems given by Eqs. (1) and(5), respectively, can be solved using integral transform methods

ð6Þ@Mðx; y; z; tÞþ LMðx; y; z; tÞ ¼ r2 Mðx; y; z; tÞ þ Gðx; y; z; tÞ@tð7Þwhere G(x, y, z, t) is the new source term which includes the original
source term as well as terms containing the ﬁltering function. The
ﬁltered initial condition is given by:Mðx; y; z; 0Þ ¼ qðx; y; zÞ  Fðx; y; z; 0Þð8ÞUsing a change-of-variable, similar to that used by [19,20,10]
we deﬁne M(x, y, z, t) in terms of a new function, h(x, y, z, t):Mðx; y; z; tÞ ¼ hðx; y; z; tÞ exp½p1 x þ p2 y þ p3 z þ tðq1 þ q2 þ q3 Þð9Þ

where p1, p2, p3, q1, q2, and q3 are constants that need to be
determined.
A comparable substitution was used by Brenner [19] to obtain atransient, one-dimensional solution for miscible ﬂuid displace-3299

J.S. Pérez Guerrero et al. / International Journal of Heat and Mass Transfer 52 (2009) 3297–3304

ment in ﬁnite length beds. Selim and Mansell [20] similarly presented an analytical solution for reactive solutes with linear
adsorption, a sink/source term, a ﬁnite domain, and continuousand ﬂux-plug type inlet conditions at the inlet boundary. Recently,
Goltz [21] used an equivalent substitution studying convective–
dispersive solute transport with constant production, ﬁrst-orderdecay, and equilibrium sorption in a porous medium. All of thesestudies used integral transforms and changes-of-variables, butwere limited to unsteady, one-dimensional problems.
The idea of the substitution in Eq. (9) is to transform the transient advection–dispersion (or diffusion) problem into an equivalent heat conduction problem involving a purely diffusive typeequation. Applying the substitution to Eq. (7) leads to:R@hðx; y; z; tÞþ hðx; y; z; tÞ½Rðq1 þ q2 þ q3 þ kÞ@t
þp1 ðDx p1 þ uÞ þ p2 ðDy p2 þ v Þ þ p3 ðDz p3 þ wÞ

To implement the GITT, we selected an eigenvalue problem
with the same kind of boundary conditions as speciﬁed forh(x, y, z, t).
In that case the problem is given by:

r2 wðx; y; zÞ þ l2 wðx; y; zÞ ¼ 0

The eigenvalue problem given by Eq. (13) has nontrivial solutions
only for certain values of the parameter l  li (i = 1, 2, . . . , 1), called
eigenvalues. The corresponding nontrivial solutions w(x, y, z) wi(x, y, z) are eigenfunctions obeying the following orthogonality
property:Zwi ðx; y; zÞwj ðx; y; zÞdv ¼ Ni dij

where Ni is the normalization integral (or the norm) and dij the Kronecker delta. Using this orthogonality property, the integral
transform pair is readily derived as:hi ðtÞ ¼@hðx; y; z; tÞ@ hðx; y; z; tÞþ ð2Dx p1 þ uÞ
 Dx
@x
@x2
Gðx; y; z; tÞ¼
exp½ðq1 þ q2 þ q3 Þt þ p1 x þ p2 y þ p3 z

v
2Dy

 2

u
q1 ¼ 
þk ;
4Dx R

;

p3 ¼

q2 ¼ 

w
2Dz

v2

ð10Þ

4Dy R

ð11a—cÞ

;q3 ¼ 

w2
4Dz R~ i ðx; y; zÞhi ðtÞ ðInverseÞ
wy; zÞ~ i ðx; y; zÞ ¼ wi ðx;
pﬃﬃﬃﬃﬃ
w
Ni

ð15aÞ
ð15bÞ

ð16Þ

The integral transformation of Eq. (12a) can now be carried out byR~ i ðx; y; zÞdv and using Eq. (15a,b), leading
applying the operator V w
to an inﬁnite system of decoupled ordinary differential equationsof the formdhi ðtÞ
 i ðtÞ;þ l2i hi ðtÞ ¼ G
dt

i ¼ 1; 2; . . .

ð17aÞ

where

 i ðtÞ ¼
GGðx; y; z; tÞ
exp½ðq1 þ q2 þ q3 Þt þ p1 x þ p2 y þ p3 zðTransformÞ

~ i ðx; y; zÞ are the normalized eigenfunctions deﬁned bywhere w

R@hðx; y; z; tÞ¼ r2 hðx; y; z; tÞ
@t

1
X
i¼1

ð11d—fÞ

Note that the coefﬁcient multiplying the term h(x, y, z, t) in Eq. (10)then reduces to zero, which allows us to write the transient transport equation for h(x, y, z, t) as:R~ i ðx; y; zÞhðx; y; z; tÞdvw

hðx; y; z; tÞ ¼

Inspection of Eq. (10) shows that the advection terms can beeliminated by choosing the constants p1, p2, p3, q1, q2 and q3 as
follows:

p2 ¼

Z
V

2

u
;
2Dx

ð14Þ

V

@hðx; y; z; tÞ@ 2 hðx; y; z; tÞþ ð2Dz p3 þ wÞ
 Dz
@z
@z2
2
@hðx; y; z; tÞ@ hðx; y; z; tÞþ ð2Dy p2 þ v Þ
 Dy
@y
@y2

p1 ¼

ð13Þ

Z

~ i ðx; y; zÞwVGðx; y; z; tÞdv
exp½ðq1 þ q2 þ q3 Þt þ p1 x þ p2 y þ p3 z

ð12aÞ

ð17bÞ

The new transport Eq. (12) is a diffusion equation that has a modiﬁed source term which contains all of the advection and decay
information from the original problem. Using Eqs. (8) and (9), themodiﬁed initial condition in terms of h(x, y, z, t) is now:

The initial conditions in Eq. (12b) must also be transformed to give:þhðx; y; z; 0Þ ¼

qðx; y; zÞ  Fðx; y; z; 0Þexpðp1 x þ p2 y þ p3 zÞ

ð12bÞThe boundary conditions must be similarly redeﬁned in terms ofh(x, y, z, t). As a general situation, we consider three kinds of homogeneous boundary conditions in only the x-direction at the genericboundary position at x = a. These conditions are summarized inTable 1, in which g and H are coefﬁcients. Boundary conditions
for the other directions can be obtained by inspection.

Table 1
Boundary conditions at the generic position x = a on the boundary for transient, M M(a, y, z, t), and steady state, M  M(a, y, z), problems.
Boundary conditions for M

Boundary conditions for h

M=0@M@x ¼ 0g @M@x þ HM ¼ 0

h=0@h@x þ p1 h ¼ 0
g @h@x þ ðH þ gp1 Þh ¼ 0

hi ðt ¼ 0Þ ¼ f i ¼

ZV~ i ðx; y; zÞ qðx; y; zÞ  Fðx; y; z; 0Þ dv
wexpðp1 x þ p2 y þ p3 zÞ

ð17cÞThe ordinary differential system Eq. (17a,b) with initial condition
Eq. (17c) has as formal solution: 2  Z t2
 i ðsÞ exp li s ds
hi ðtÞ ¼ exp  li t f i þ 1
G
R 0
R
R

ð18Þ

The unknowns h(x, y, z, t) and T(x, y, z, t) can now be obtained byusing the inverse formula given by Eq. (15b), and Eqs. (6), (9),respectively, to give:hðx; y; z; tÞ 2  Z t
1
2
X
 i ðsÞ exp li s ds ð19Þ
~ i ðx; y; zÞ exp  li t f i þ 1
¼
w
G
R 0
R
R
i¼1

 2

ux
v y wzuþkTðx; y; z; tÞ ¼ hðx; y; z; tÞ expþþt
2Dx 2Dy 2Dz
4Dx R

v2
w2
þ Fðx; y; z; tÞð20Þþþ
4Dy R 4Dz R

3300

J.S. Pérez Guerrero et al. / International Journal of Heat and Mass Transfer 52 (2009) 3297–3304

3.3. Steady-state regime

4. Test cases

Eq. (5) deﬁned a steady-state transport problem with nonhomogeneous boundary conditions. The boundary conditions canagain be made homogeneous using the ﬁlter strategy:

4.1. First test caseTðx; y; zÞ ¼ Mðx; y; zÞ þ Fðx; y; zÞð21Þ

where F(x, y, z) satisﬁes exactly the original boundary conditions ofT(x, y, z). Therefore, the steady-state problem can be re-written as

As a test case for the general analytical solution, we consider theparticular problem of solving the linearized Burgers equation [13].
The partial differential equation for this test case isR@Tðx; tÞ@Tðx; tÞ@ 2 Tðx; tÞ
þu
¼ Dx
@t
@x
@x2

ð30aÞ

with initial and boundary conditions:LMðx; y; zÞ ¼ r2 Mðx; y; zÞ þ Gðx; y; zÞð22ÞTðx; 0Þ ¼ 1;

0x1

where G(x, y, z) is the new source term containing information about
the original source and the ﬁlter function.
For the steady-state problem we use a change-of-variable similar to Eq. (9) but without the time domain:Tð0; tÞ ¼ 1;

Tð1; tÞ ¼ 0;

Mðx; y; zÞ ¼ hðx; y; zÞ expðp1 x þ p2 y þ p3 zÞu2
v 2 w2þþþ kR hðx; y; zÞ
4Dx 4Dy 4Dz¼ r2 hðx; y; zÞ þ

hi ¼

1
X

~ i ðx; y; zÞhiwðTransformÞ
ðInverseÞ

u2
v 2 w2
 i ; i ¼ 1; 2; . . .
þþþ kR hi þ l2i hi ¼ G
4Dx 4Dy 4Dz
ZGðx; y; zÞi ¼~ i ðx; y; zÞ
dv
w
G
u
w
V
exp 2Dx x þ 2Dv y y þ 2D
z
z

i
G

ð26aÞ

Tðx; y; zÞ ¼ hðx; y; zÞ exp

ð33c; dÞ

u
;
2Dx

u2
q1 ¼ ðÞ
4Dx R

ð34a; bÞbi ¼ ip;pﬃﬃﬃﬃﬃﬃ

li ¼ bi Dx ; i ¼ 1; 2; . . .

ð35a; bÞ

Ni ¼

1
;
2

~i ¼
w

pﬃﬃﬃ
2sinðbi xÞð36; 37Þ

Z 1~ j ðxÞdx ¼ dij~ i ðxÞw
w

ð38Þ

0

i
G

l

t>0

and with the norms and the normalized eigenfunctions as:

Finally, invoking the inverse formula, Eq. (25b), and recalling Eqs.
(21), (23), we obtain:Mð1; tÞ ¼ 0;

ð33bÞ

According to the general solution presented above, we need to
specify the eigenfunction and eigenvalues problem in a form suchas deﬁned by Eq. (13). In this case the eigenvalue problem is a
Sturm–Liouville problem with the set of eigenvalues given by:

ð26bÞ

ð27Þ

v 2 þ w2 þ kR
l2i þ 4Dx þ 4D
4Dz
y

i¼1

Mð0; tÞ ¼ 0;

ð33aÞWe note that these eigenfunctions and eigenvalues satisfy the following orthogonality property:

u2

2
i þ

@Mðx; tÞ@Mðx; tÞ@ 2 Mðx; tÞ
þu
¼ Dx
@t
@x
@x2
Mðx; 0Þ ¼ 1  FðxÞ; 0  x  1

R

p1 ¼

The solution of this equation is:~ i ðx; y; zÞw

ð32Þ

1  exp  Dux

ð25bÞ



1
X

1  exp  Dux ð1  xÞ

The values of p1 and q1 that transform Eq. (33a) into an exclusively diffusive equation are given by:

Note that the transformed potentials 
hi are constants for each valueof i.
R~ i ðx; y; zÞdv to Eq. (24) and using Eq.
Applying the operator V w(25a,b), results in the following transformed equation:

hðx; y; zÞ ¼

ð31b; cÞ

ð25aÞ

i¼1

hi ¼

Fð1Þ ¼ 0

Eq. (30) can be written in terms of M(x, t) by using Eq. (6):Vhðx; y; zÞ ¼

Fð0Þ ¼ 1;

ð24Þ

u
w
x þ 2Dv y y þ 2Dz2Dxz

~ i ðx; y; zÞhðx; y; zÞdv
w

ð31aÞ

with boundary conditionsFðxÞ ¼

Eq. (24) is also a diffusion equation, but with a modiﬁed source termthat contains the advection information of the original problem. The
boundary conditions in terms of h(x, y, z) are summarized in Table 1.
The eigenvalue problem in this case is the same as for the transient problem, and is given by Eq. (13). The integral transform pairis given by:

Z

dFðxÞ
d FðxÞ
¼ Dx
2
dx
dx

This differential equation has the following analytic solution:

Gðx; y; zÞexp

ð30c; dÞ

2

u



t>0

Because Eq. (30c) is not homogeneous, it is necessary to deﬁne a
ﬁlter function to homogenize the boundary condition. A suitable
ﬁlter may be found by solving the steady-state version of Eq. (30a):ð23ÞThis substitution eliminates the advection terms when p1, p2, and p3
are chosen judiciously. In this case the choices are the same as forthe transient case, Eq. (11a–c). Therefore, the new formulation forthe steady transport equation in terms of h(x, y, z) is

ð30bÞ

u2
v 2 þ w2 þ kR
þ 4D
4Dx
4Dz
y

ð28Þ



u
v
w
xþ
yþ
z þ Fðx; y; zÞ
2Dx
2Dz
2Dy

ð29Þ

The next step according to the general analytical solution procedure is to calculate the initial coefﬁcient f i established in Eq. (17c),
whose analytic expression is given by

pﬃﬃﬃ 2iuf ¼  ð1Þ 4 2Dx li expð 2Dx Þi
2 22u þ 4Dx li

ð39Þ

The structure of Eq. (37) implies that for speciﬁed values of u andDx, the absolute value of f i decreases monotonically as the eigen-

3301

J.S. Pérez Guerrero et al. / International Journal of Heat and Mass Transfer 52 (2009) 3297–3304

Tðx; tÞ ¼ exp



X


1
2
u
u2
~ i ðxÞ exp  li t f i þ FðxÞ
x
w
t
2Dx
4Dx R
R
i¼1
ð40Þ

Numerical results were obtained for the following two sets ofparameter values: t = 0.1, R = 1, Dx = 1, u = 1, k = 0 and t = 0.1, R = 1,
Dx = 1, u = 10, k = 0. These cases were chosen to permit a comparison
with previous results presented by Cotta [13]. Cotta obtained a solution for the linearized Burgers equation using a hybrid analytical–
numerical GITT approach in which the transformed inﬁnite systemof ordinary differential equations was truncated and solved numerically using the DIVAPG subroutine from the IMSL Library [22].
Tables 2 and 3 compare the convergence of the present solution
with the hybrid solution of Cotta [13], the latter results shown inparentheses. The parameter N in the tables is the number of terms
summed in the truncated inﬁnite series for both the present analytic solution and the hybrid solution of [13]. Also shown are results for a fully numerical solution obtained by Cotta [13] using
the DMOLCH subroutine from the IMSL Library [22]. In both cases(with u = 1 and u = 10) the new analytical solution required only
N = 5 terms to achieve convergence to six decimal places. In fact,
the solution for only N = 1 term already provides an excellent
approximation to the true solution. The converged analytical solution agreed excellently with the hybrid GITT and the full numericalresults reported by Cotta [13].
The faster convergence demonstrated in Tables 2 and 3 is due tothe change-of-variable used in the present analytic solution thattransformed the original advection–diffusion partial differential
equation into a pure diffusive equation. The hybrid solution of[13] applied the GITT directly to the advection–diffusion equation.
4.2. Second test caseAs a second example, we solve a transport problem that, among
other applications, has been employed to model nutrient and contaminant transport in soils (e.g., van Genuchten [23]). In his paper,van Genuchten [23] used Laplace transforms to derive the analytic

solution for solute transport of up to four species involved insequential decay chain reactions. Here, we compare results fromthe present formulation with earlier results obtained in [23] forthe ﬁrst species, ammonium (NH4+). The transport of ammonium
is of interest both as a plant nutrient and a possible groundwatercontaminant.
The dimensionless transport equation, in terms of the soluteconcentration, is given by:R@Cðx; tÞ @Cðx; tÞ
1 @ 2 Cðx; tÞþþ cCðx; tÞ ¼
@t@x
Pe @x2

ð41aÞ

with initial condition:

Cðx; 0Þ ¼ 0;

t¼0

ð41bÞ

and boundary conditions:@Cð0; tÞ
þ PeCð0; tÞ ¼ Pe;@x@Cð1; tÞ¼0@x

ð41c; dÞThe dimensionless parameters Pe and c in Eq. (39) are deﬁned by:

Pe ¼

L0 u
;
Dx

DIMENSIONLESS CONCENTRATION

value order increases. Such behavior is expected in the GITT
approach.
We can now compose the analytical solution for the linearized
Burgers equation:c¼

k RL0uð42a; bÞ1

Steady-state

0.8

0.6

20 h
0.4

15 h
10 h

0.2

t=5h
0
0

4

8

12

16

20

DISTANCE (cm)Fig. 1. Distribution of the dimensionless concentration at various times(L0 = 20 cm).

Table 2
Convergence comparison of the analytical and hybrid solutions for u = 1, t = 0.1, R = 1, Dx = 1, and k = 0.
x

0.1
0.3
0.5
0.7
0.9

Analytical solution (hybrid solution, Ref. [13])Ref. [13]

N=1

N=5

N = 10N = 15

N = 300.983264
0.925062
0.798233
0.567179
0.216888

0.981048 (0.98145)
0.921078 (0.92105)
0.798211(0.79842)
0.57206 (0.57225)
0.220238 (0.21993)

0.981048 (0.98101)
0.921078 (0.92107)
0.798211(0.79823)
0.57206 (0.57211)
0.220238 (0.22033)

0.981048 (0.98104)
0.921078 (0.92109)
0.798211(0.79821)
0.57206 (0.57207)
0.220238 (0.22026)

0.981048 (0.98105)
0.921078 (0.92108)
0.798211(0.79821)
0.57206 (0.57206)
0.220238 (0.22024)

0.98100.92100.79810.57200.2202

Table 3
Convergence comparison of the analytical and hybrid solutions for u = 10, t = 0.1, R = 1, Dx = 1, and k = 0.
x

0.1
0.3
0.5
0.7
0.9

Analytical solution (hybrid solution, Ref. [13])Ref. [13]

N=1

N=5

N = 10N = 15

N = 300.999941
0.999268
0.99376
0.951251
0.633182

0. 999939 (1.00020)
0.999259 (0.99953)
0.99376 (0.99476)
0.951317 (0.95264)
0.633293 (0.63477)

0.999939 (0.99993)
0.999259 (0.99923)
0.99376 (0.99370)
0.951317 (0.95122)
0.633293 (0.63322)

0. 999939 (0.99994)
0.999259 (0.99927)
0.99376 (0.99377)
0.951317 (0.95135)
0.633293 (0.63322)

0.999939 (0.99994)
0.999259 (0.99926)
0.99376 (0.99376)
0.951317 (0.95132)
0.633293 (0.63329)

0.99990.99930.99380.95130.6333

3302

J.S. Pérez Guerrero et al. / International Journal of Heat and Mass Transfer 52 (2009) 3297–3304

where L0 is the length of the domain, u* is he advective velocity, Dxis the dispersion coefﬁcient, k is the ﬁrst-order decay coefﬁcient,and R is the retardation factor. Parameter values for the transportof ammonium were taken as [23]: u*=1 cm h1, Dx ¼ 0:18 cm2 h1,
R = 2, k ¼ 0:005 h1.
The parameters of the general analytical solution and the variables in the present test problem correspond as:

Dx ¼

1
;
Pe

u ¼ 1;

k¼

cR

ð43a-cÞ

Pe;2

q1 ¼ 



1 Pe
þc
R 4

wi ¼ bi cosðbi xÞ þ H1 sinðbi xÞð44a; bÞThe eigenfunction for this case is obtained from the following
Sturm–Liouville problem:

ð46Þ

The norms and eigenvalues are obtained from the following
equations:

ðb2i þ H21 Þ þ H1 þ H22
bi ðH1 þ H2 Þtanðbi Þ ¼ 2
bi  H1 H2
Ni ¼

Therefore, parameter p1 and q1 are given by:

p1 ¼

It is interesting that, due to the change-of-variables, the boundary
condition at x = 1 is now a third-type boundary condition, similarly
as the condition at x = 0.
The analytic solution of the eigenvalue problem is [10]:

ð47Þ
ð48Þ

respectively, where H1 ¼ Pe
andH2 ¼ Pe.
2
2
The ﬁlter function is obtained by solving the equation:22d wðxÞ
2dx

þ l2 wðxÞ ¼ 0

ð45aÞ

dwð0Þ
þ Pewð0Þ ¼ 0;
dx

ð49aÞ

with the boundary conditions:with boundary conditions:



dFðxÞ
1 d FðxÞþ cFðxÞ ¼
dx
Pe dx2

dwð1Þ
þ Pewð1Þ ¼ 0dx

ð45b; cÞdFð0Þ
þ PeFð0Þ ¼ Pe;
dx

dFð1Þ
¼0
dx

ð49b; cÞThis ordinary differential equation can be solved analytically, leading to the following expression for the ﬁlter function:

Table 4
Solution convergence for L0 = 200 cm and t = 200 h (N = number of terms summed).
X (cm)

Dimensionless concentration

Ref.[23]

N = 250

N = 350

N = 400

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

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.675789878  106
1.045824992  1079.845112917  10100.000000000
1.325484431  10100.0001250537544
336.4850318
3.536666357  108
1.365890562  1014
2.043458682  1020
4.597385308  1026
4.459462906  1032
1.312702502  1038
3.199088984  1044
6.267155179  1050
5.611318166  1056

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.675789878  106
1.045824992  1079.845112917  10100.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000

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.675789878  106
1.045824992  1079.845112917  10100.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000
0.000000000

0.998210.94961
0.903380.859400.817560.777750.73989
0.70387
0.66960
0.63700
0.60599
0.576480.548420.52172
0.49632
0.472150.44901
0.425080.389430.314900.19272
0.07679
0.01794
0.002310.000160.00001
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000

3303

J.S. Pérez Guerrero et al. / International Journal of Heat and Mass Transfer 52 (2009) 3297–3304hhpﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ i
pﬃﬃﬃﬃﬃ ipﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
4c Pe exp w Pe þ 12 Pe þ w Pe x þ 2 exp 12 Pe  w Pe x Pe Pe þ w Pe þ 2cFðxÞ ¼
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
Pe þ w Pe þ w Pe þ 2c
 exp w Pe
Pe þ w Pe  w Pe þ 2c þpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
where w ¼ Pe þ 4c.
Finally, we ﬁnd the transformed initial condition coefﬁcient:f ¼ num
ð51aÞ
i
den
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
num ¼ 8 expðw PeÞPeðw  PeÞcbi  4Peðw þ PeÞðPe þ w Pe þ 2cÞbihpﬃﬃﬃﬃﬃ
w pﬃﬃﬃﬃﬃ
þ 2 exp Pe sinðbi Þ Pe5=2 Pe þ w Pe þ 4c
2
pﬃﬃﬃﬃﬃ iþ4wðPe þ w PeÞb2i
ð51bÞ
hpﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃi
den ¼ ðw þ PeÞðPe þ w Pe þ 2c  exp w Pe
pﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃ
 ðPe  w Pe þ 2cÞÞðPew2 þ 4b2i Þ Ni
ð51cÞThe symbolic and numerical computations were made in the Mathematica platform [24]. When we used for this purpose the Mathematica function FindRoot to solve the transcendental Eq. (48) andcompute the eigenvalues, it was necessary to set the Mathematica
parameter WorkingPrecision to 200.
Fig. 1 shows dimensionless concentration proﬁles computed for
different times with a relatively short domain length of L0 = 20 cm.
The ﬁgure shows the concentration distribution progressing toward a linear, steady-state proﬁle.
Table 4 illustrates the convergence of the solution computed fort = 200 h, where N is again the number of terms summed in the truncated series expansion. For comparison purposes, all values less than
1010 were discarded. The results show that convergence is obtained
for the entire spatial domain of L0 = 200 cm with N = 350 terms; theconcentration values did not change when additional (N = 400 ormore) terms were used. The converged values are in complete agree-

ð50Þ

ment with the previous results of van Genuchten [23]. Table 4 also
demonstrates that solution convergence slowly progressed towards
the end of the spatial domain as the number of terms in the series increased. For example, convergence for X 6 135 cm was obtained
Table 6
Solution convergence for L0 = 20 cm and t = 20 h (N = number of terms summed).
X (cm)

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

Dimensionless concentration

Ref. [25]N = 20N = 50N = 100

0.998206
0.988291
0.978469
0.968683
0.958554
0.946242
0.925461
0.881528
0.792957
0.646514
0.458117
0.268772
0.1758860.639369
10.673
164.519
2538.9
39170.
604111.
9.31184  106
1.43399  108

0.998206
0.988291
0.978469
0.968683
0.958554
0.946242
0.925461
0.881528
0.792956
0.646526
0.457931
0.271654
0.131256
0.0506341
0.0153803
0.00364344
0.000668586
0.0000945846
0.0000102798
8.55118  1076.81699  1080.998206
0.988291
0.978469
0.968683
0.958554
0.946242
0.925461
0.881528
0.792956
0.646526
0.457931
0.271654
0.131256
0.0506341
0.0153803
0.00364344
0.000668586
0.0000945846
0.0000102798
8.55118  1076.81699  1080.998210.98829
0.97847
0.96868
0.95855
0.946240.92546
0.88153
0.792960.646530.457930.27165
0.13126
0.050630.01538
0.003640.000670.00009
0.00001
0.00000
0.00000

Table 5
Solution convergence for L0 = 140 cm and t = 200 h (N = number of terms summed).
X (cm)

Dimensionless concentration
N = 100

N = 150

N = 200

N = 250

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

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.4490137609
0.1320360358
103148.4364
3.904198455  1011
1.109849736  1017
3.774497734  1023
3.448237131  1029
2.298888763  1035
5.264401526  1041
4.119631336  1046
5.803900012  1053
3.829622487  1059
4.484729306  1065

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.0001586227102
0.02564681436
2020.240884
2.788998787  1010
2.605574947  1016

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.675789878  106
1.045824992  1079.845112917  101000.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.675789878  106
1.045824992  1079.845112917  1010
0

3304

J.S. Pérez Guerrero et al. / International Journal of Heat and Mass Transfer 52 (2009) 3297–3304

with N = 250 terms, while N = 350 terms were required to achieve
convergence over the entire 200 cm domain.
Because in this example the concentration was zero for
X > 140 cm, the simulation could have been performed with a shorter domain (e.g., L0 = 140 cm). Table 5 shows that convergence for
L0 = 140 cm is faster, occurring now with N = 200. The faster convergence is due to the fact that a smaller value of L0 corresponds to asmaller Peclet number, which causes the transport problem to become more diffusive. This faster convergence for smaller L0 can be
exploited when making early time calculations for large domains(i.e., when the solute has moved in only a small fraction of the transport domain). This idea is demonstrated in Table 6, which shows results for the same transport parameters as before at the early timet = 20 h and L0 = 20 cm. Convergence was achieved now with only
N = 50 terms and veriﬁed with N = 100. We note that the converged
results duplicated exactly the earlier analytic solution [23] as implemented within the STANMOD computer software [25]. We note alsothat it was possible to obtain these same results with default values
for the WorkingPrecision parameter.
5. Conclusions
Using the Generalized Integral Transform Technique (GITT), in
its classic formulation, in combination with a simple algebraic substitution, it was possible to obtain a formal exact solution of the linear advection–dispersion (or diffusion) transport equation for bothtransient and steady-state regimes. The mathematical substitution,
which transformed the original advection–diffusion problem, into
an exclusively diffusive problem, proved to be very advantageousby improving convergence of the GITT series solution.
Acknowledgments
The ﬁrst author acknowledges support by the Brazilian National
Research Council (CNPq), the Department of Environmental Sciences of the University of California Riverside (UCR) and the UnitedStates Salinity Laboratory (USSL). This research was supported alsoin part by a fellowship from CAPES to M. Th. van Genuchten.
References[1] M.Th. van Genuchten, Analytical solutions for chemical transport withsimultaneous adsorption, zero-order production and ﬁrst-order decay, J.
Hydrol. 49 (1981) 213–233.

[2] W.A. Jury, W.F. Spencer, W.J. Farmer, Behavior assessment model for traceorganics in soil: I. Model description, J. Environ. Qual. 12 (1983) 558–564.
[3] I. Javandel, C. Doughty, C.F. Tsang, Groundwater Transport: Handbook ofMathematical Models, Water Resources Monograph Series, vol. 10, American
Geophysical Union, Washington, DC.
[4] F.J. Leij, T.H. Skaggs, M.Th. van Genuchten, Analytical solutions for solutetransport in three-dimensional semi-inﬁnite porous media, Water Resour. Res.
27 (10) (1991) 2719–2733.
[5] C.R. Quezada, T.P. Clement, K.K. Lee, Generalized solution to multi-dimensional
multi-species transport equations coupled with a ﬁrst-order reaction network
involving distinct retardation factors, Adv. Water Resour. 27 (5) (2004) 508–
521.
[6] D. Courant, D. Hilbert, Methods of Mathematical Physics, Wiley Interscience
Publications, New York, USA, 1953.
[7] P.M. Morse, H. Feshbach, Methods of Theoretical Physics, McGraw-Hill, New
York, USA, 1953.
[8] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Oxford
University press, Oxford, 1959.
[9] I.N. Sneddon, The Use of Integral Transforms, McGraw-Hill, 1972.
[10] M.N. Ozisik, Heat Conduction, Wiley, New York, 1980.
[11] F.J. Leij, M.Th. van Genuchten, Analytical modeling of nonaqueous phase liquid
dissolution with Green’s functions, Transport in Porous Media 38 (2000) 141–
166.
[12] M.D. Mikhailov, M.N. Ozisik, Uniﬁed Analysis and Solutions of Heat and Mass
Diffusion, John Wiley & Sons, 1984.
[13] R.M. Cotta, Integral Transforms in Computational Heat and Fluid Flow, CRC
Press, Boca Raton, FL, 1993.
[14] D. Zwillinger, Handbook of Differential Equations, Academic Press, 1992.
[15] A.D. Polyanin, Handbook of Linear Partial Differential Equations for Engineers
and Scientists, Chapman & Hall/CRC Press, Boca Raton, 2002.
[16] R.B. Codell, K.T. Key, G.A Whelan, Collection of Mathematical Models for
Dispersion in Surface Water and Groundwater, NUREG 0868, U.S. Nuclear
Regulatory Commission, Washington, DC, 1982.
[17] M. Th. van Genuchten, W.J. Alves, Analytical Solutions of the One-Dimensional
Convective–Dispersive Solute Transport Equation, USDA ARS Technical
Bulletin Number 1661, U.S. Salinity Laboratory, 1982.
[18] E.J. Wexler, Analytical Solutions for One-, Two-, and Three-Dimensional Solute
Transport in Ground-Water Systems with Uniform Flow, U.S. Geological
Survey Techniques of Water-Resources Investigations, 1992. Book 3, Chapter
B7.
[19] H. Brenner, The diffusion model of longitudinal mixing in beds of ﬁnite length.
Numerical values, Chem. Eng. Sci. 17 (1962) 229–243.
[20] H.M. Selim, R.S. Mansell, Analytical solution of the equation of reactive solutes
through soils, Water Resour. Res. 12 (1976) 528–12532.
[21] W.J. Golz, Solute Transport in a Porous Medium: A Mass-Conserving Solution
for the Convection–Dispersion Equation in a Finite Domain. Ph.D. thesis,
Louisiana State University, 2003.
[22] IMSL Library, MATH/LIB, Houston, TX, 1987.
[23] M.Th. van Genuchten, Convective–dispersive transport of solutes involved insequential ﬁrst-order decay reactions, Comput. Geosci. 11 (2) (1985) 129–147.
[24] MATHEMATICA, Wolfram Research, Inc., Champaign, IL, 2004.
[25] J. Simunek, M. Th. van Genuchten, M. Sejna, N. Toride, F.J. Leij, The STANMOD
computer software for 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, Colorado, 1999.