Transport in Porous Media 38: 141–166, 2000.
c 2000 Kluwer Academic Publishers. Printed in the Netherlands.

141

Analytical Modeling of Nonaqueous Phase
Liquid Dissolution with Green’s Functions
FEIKE J. LEIJ and MARTINUS TH. VAN GENUCHTEN
U.S. Salinity Laboratory, U.S. Department of Agriculture, Agricultural Research Service, 450 W.
Big Springs Road, Riverside, CA 92507-4617, U.S.A., e-mail: fleij@ussl.ars.usda.gov(Received: 28 July 1998)
Abstract. Equilibrium and bicontinuum nonequilibrium formulations of the advection–dispersion
equation (ADE) have been widely used to describe subsurface solute transport. The Green’s Function Method (GFM) is particularly attractive to solve the ADE because of its flexibility to deal
with arbitrary initial and boundary conditions, and its relative simplicity to formulate solutions
for multi-dimensional problems. The Green’s functions that are presented can be used for a wide
range of problems involving equilibrium and nonequilibrium transport in semi-infinite and infinite
media. The GFM is applied to analytically model multi-dimensional transport from persistent solutesources typical of nonaqueous phase liquids (NAPLs). Specific solutions are derived for transport
from a rectangular source (parallel to the flow direction) of persistent contamination using first-,
second-, or third-type boundary or source input conditions. Away from the source, the first- andthird-type condition cannot be expected to represent the exact surface condition. The second-typecondition has the disadvantage that the diffusive flux from the source needs to be specified a priori.
Near the source, the third-type condition appears most suitable to model NAPL dissolution into
the medium. The solute flux from the pool, and hence the concentration in the medium, depends
strongly on the mass transfer coefficient. For all conditions, the concentration profiles indicate that
nonequilibrium conditions tend to reduce the maximum solute concentration and the total amountof solute that enters the porous medium from the source. On the other hand, during nonequilibrium transport the solute may spread over a larger area of the medium compared to equilibriumtransport.
Key words: Green’s functions, NAPL dissolution, nonequilibrium transport, initial value problem,
boundary value problem.
Abbreviations: ADE – advection–dispersion equation; GFM – Green’s function method.

Nomenclature
a
b
c
C
D
f
G
g

auxiliary variable.
auxiliary variable.
dimensional concentration [ML−3 ].
dimensionless solute concentration (=c/c0 ).
diffusion coefficient [L2 T−1 ].
initial concentration.
Green’s function.
input concentration for NAPL pool.

142

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTENh
I0 , I1
k
k∗
L
P
R
s
s∗
T
u
W
X
Y
Z

function in convolution integral.
modified Bessel function of orders zero and one.
dimensionless mass-transfer coefficient.
mass-transfer coefficient [LT−1 ].
characteristic length.
Peclet number.
retardation factor.
Laplace variable.
modified Laplace variable.
dimensionless time.
arbitrary function.
dimensionless vertical velocity.
dimensionless longitudinal position.
dimensionless transversal position.
dimensionless transversal position.

Greek
β
0
δ
ζ
η
µ
ξ
τ
ω
3

dimensionless variable for partitioning in nonequilibrium transport models.
prescribed gradient for second-type condition.
Dirac delta function.
dummy variable.
dummy variable.
first-order degradation term.
dummy variable.
dummy variable.
dimensionless mass transfer coefficient.
Auxiliary term in Green’s function to quantify nonequilibrium effects.

Subscriptseeffective value.
Llongitudinal.
0reference or characteristic value.
ssolubility or upper value.
Ttransversal.
uequilibrium solution for unit input.
X
X-direction.
YY-direction.
Z
Z-direction.
1
mobile region in two-region nonequilibrium model.
2
immobile region in two-region nonequilibrium model.
Superscripts
ffirst-type condition.
ssecond-type condition.
tthird-type condition.
Special Symbols
L
differential operator for advection–dispersion equation.
adjoint operator for advection–dispersion equation.
L∗L or− Laplace transform.

ANALYTICAL MODELING OF NAPL DISSOLUTION

143

1. Introduction
The deterministic advection–dispersion model has been widely used to describe
the movement of dissolved chemicals, heat and other substances in porous media.
Mathematical solutions of this model allow us to quantify the spatial and temporal
distribution of the solute concentration in porous media, thus providing important
information for optimally monitoring or managing the subsurface environment.
Whereas, numerical solutions are often used for general problems, analytical solutions may be employed for idealized conditions. Analytical solutions are valuable
to evaluate numerical solution techniques, gain a better understanding of the importance of different transport parameters and mathematical conditions, estimate
model parameters with inverse methods, and study transport for conditions where
numerical methods may not yield accurate or reliable results. Techniques to derive
analytical solutions of the advection–dispersion equation (ADE) typically involve
Fourier analysis and integral transformation procedures (cf. Cleary and Adrian,
1973; Güven et al., 1984; Leij et al., 1993). Relatively few solutions have beenobtained with Green’s functions (Yeh and Tsai, 1976; Galya, 1987; Ellsworth andButters, 1993). Yeh (1981) presented several Green’s functions for a source termin the ADE.
The introduction and movement of nonaqueous phase liquids (NAPLs) in thesubsurface presents an especially serious threat to the long-term quality of soil
and groundwater systems. Spills of hydrocarbon solvents can create hazardous andlong-lasting sources of groundwater contamination (Bradford et al., 1998). Due toits density and viscosity, a dense nonaqueous phase liquid (DNAPL) tends to sinkto the bottom of aquifers and often remains there as a contamination pool becauseof its low solubility in water. Similarly, a light nonaqueous phase liquid (LNAPL)
may form a relatively immobile pool on top of the aquifer. Pollution from oil tanks,
landfills, and other buried sources are further examples of gradual and continuous
contamination from a fixed source. In this study we use the advection–dispersion
model to describe NAPL movement in the aqueous phase after dissolution from a
simple rectangular pool at the top or bottom of the aquifer.
Analytical solutions for such persistent contamination problems can be derived
by integrating instantaneous sources with respect to time or space (e.g., Prakash,
1984; Yeh, 1981). This approach does not allow the specification of different typesof boundary conditions at the interface of the source and the porous medium.
Persistent NAPL contamination has been described as a boundary value problem
by, among others, Chrysikopoulos et al. (1994), Holman and Javandel (1996), andShan and Javandel (1997). Solutions in these studies were obtained with a variety
of integral transforms. However, the solutions may be obtained far more conveniently with the Green’s function method (GFM). This method offers flexibility interms of handling changes in initial and boundary conditions. Furthermore, very
general and concise expressions for the solute concentration are possible with the
GFM, particularly for multidimensional problems. Several pertinent texts exist on

144

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

the method (Greenberg, 1971; Stakgold, 1979; Roach, 1982; Beck et al., 1992).
Expressions for Green’s functions have been derived for the standard ADE (Leijet al., 1999). However, bicontinuum transport models may predict solute transport
more accurately by distinguishing a ‘mobile’ or flowing region, with advective
and dispersive solute transport, and an ‘immobile’ or stagnant region (possibly a
different fluid phase). Solute exchange between the regions is described by a firstorder diffusion equation (Coats and Smith, 1964; Van Genuchten and Wierenga,1976).
The objectives of this study are to (1) formulate Green’s functions for the equilibrium and nonequilibrium ADE, (2) use the GFM to derive analytical solutions
for solute movement from a persistent contaminant source, and (3) illustrate theeffects of alternative source boundary conditions, mass transfer parameters, andnonequilibrium parameters on the solute distribution.

2. Mathematical Formulation
2.1.

EQUILIBRIUM TRANSPORT

The dimensionless ADE for a solute subject to unidirectional flow, three-dimensional dispersion, linear sorption and first-order degradation in a porous medium
whose transport properties are constant in time and space, may be given
by
R

∂C
∂C
1 ∂ 2C1 ∂ 2C1 ∂ 2C
−
+
− µC,
=
+
∂T
PX ∂X 2 ∂X PY ∂Y 2
PZ ∂Z 2(1)

where R is a solute retardation factor, T is time, X, Y and Z denote the position
along the longitudinal and two transversal coordinates with corresponding Peclet
numbers PX , PY and PZ , and µ is a degradation factor. The model may be extended
to include a zero-order source/sink term without complicating the solution procedure according to the GFM. The transformation from dimensional to dimensionlessparameters has been shown elsewhere (e.g., Toride et al., 1993).
2.1.1. Green’s Function Method
The differential operator, L, for the previously defined ADE is
L=R

∂
∂
1 ∂2
1 ∂2
1 ∂2
+
−
+ µ.
−
−
∂T
PX ∂X 2 ∂X PY ∂Y 2 PZ ∂Z 2

(2)

The solution to L(C) = 0 yields the solute concentration C. An alternative
adjoint problem, which may be easier to solve, is obtained through multiplication of L(C) by a Green’s function G and subsequent integration over the solution
region (X, Y, Z, T). This procedure can be schematically represented as

ANALYTICAL MODELING OF NAPL DISSOLUTION

145(cf. Greenberg, 1971)
ZZZZ
GL(C) dζ dη dξ dτZZZZ= boundary terms +CL∗ (G) dζ dη dξ dτ.

(3)

Integration by parts yields the boundary terms and the formally adjoint differential
operator, which is defined as
L∗ = −R

∂
∂
1 ∂2
1 ∂2
1 ∂2
−
−
+ µ.
−
−
∂τ
PX ∂ξ 2 ∂ξ
PY ∂η2 PZ ∂ζ 2

(4)

The dummy variables (ξ , η, ζ , τ ) correspond to the regular independent variables(X, Y, Z, T). The adjoint operator is similar to the ADE operator except for thesign of the storage and advection terms. A solution for C may be obtained fromthe equivalent problem in terms of the adjoint operator by requiring that the adjoint
differential operator on G is equal to the following four-dimensional Dirac delta
function:L∗ (G) = δ(X − ξ, Y − η, Z − ζ, T − τ ).

(5)

The Green’s function G(X, Y, Z, T; ξ , η, ζ , τ ) denotes the concentration at (X, Y, Z,
T) as the result of instant solute release at τ for a unit source located at (ξ , η, ζ ).
The concentration for a spatial domain with longitudinal boundaries at X1 , X2 andtransversal boundaries at Y1 , Y2 , and Z1 , Z2 , is expressed according to the GFM as ξ →X2Z T Z Y2 Z Z2 
G ∂C
C ∂G
C(X, Y, Z, T ) =dζ dη dτ+
−
− GC
PX ∂ξ
PX ∂ξ
Y1
Z1
0
ξ →X1
 η→Y2Z T Z X2 Z Z 2 
G ∂C
C ∂G
+
dζ dξ dτ +−PY ∂η
PY ∂η η→Y1
X1
Z1
0
 ζ →Z2Z T Z X2 Z Y2 
G ∂C
C ∂G
+
dη dξ dτ +−PZ ∂ζ
PZ ∂ζ ζ →Z1
X1
Y1
0
Z X2 Z Y2 Z Z2+
RGCdζ dη dξ.
(6)
X1

Y1

Z1

τ =0

Equation (6) follows from Equation (3) when L(C) = 0 and L∗ (G) is given by(5). The first three terms on the right-hand side arise for boundary value problems in the X, Y and Z directions, while the fourth term is associated with theinitial value problem; this expression for C can be further simplified once theboundary and initial conditions are known. The appropriate Green’s function isdetermined by solving Equation (5). This equation is homogeneous except for a
discontinuity when (ξ = X, η = Y , ζ = Z, τ = T ). Because the boundary conditions for G are homogeneous, the Green’s function may be written as

146

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

the product of Green’s functions for individual directions (cf. Beck et al., 1992):
G(X, Y, Z, T ; ξ, η, ζ, τ ) = GX (X, T ; ξ, τ ) GY (Y, T ; η, τ ) GZ (Z, T ; ζ, τ ). Notethat a similar application of the product rule may not be permitted for the concentration (Carslaw and Jaeger, 1959). The degradation factor for three-dimensionaltransport is partitioned, using somewhat arbitrary contributions for each direction,such that µ = µX + µY + µZ . In the following, mostly unidirectional Green’s
functions will be considered for the purpose of constructing a three-dimensional
solution for the concentration.
2.1.2. Green’s Functions
The following general concentration for the longitudinal direction, resulting fromone-dimensional advective–dispersive transport for equilibrium conditions, can beinferred from (6):
 X2
Z T
GX ∂GXCX ∂GXCX (X, T ) =dτ +−− GX CX
PX ∂ξ
PX ∂ξ
0
X1
Z X2
+
RGX (X, T ; ξ, 0)CX (ξ, 0) dξ,(7)
X1

where the boundaries X1 and X2 are specified depending on the solution domain. A
unique solution may be obtained by applying the boundary conditions for CX , andsubsequently imposing conditions on GX to be used in the solution of the adjoint
problem:L∗ (GX ) = −R

∂GX
∂GX
1 ∂ 2 GX−−
+ µX GX = δ(X − ξ, T − τ ).
2
∂τ
PX ∂ξ
∂ξ
(8)

Solutions for GX can often be readily inferred from existing solutions for the ADE(cf. Van Genuchten and Alves, 1982; Leij et al., 1999).
The concentration in one of the transversal directions, say Y, as a result of onedimensional dispersion may be expressed according to the GFM as (cf. Equation
(6)):

Z T
GY ∂CYCY ∂GY Y2
CY (Y, T ) =
dτ +
−
PY ∂η
PY ∂η
0
Y1
Z Y2
+
RGY (Y, T ; η, 0)CY (η, 0) dη.
(9)
Y1

Explicit expressions for the concentration in the transversal Y direction are obtainedby imposing a zero-gradient condition on CY at a finite or infinite distance and by
substituting the appropriate Green’s function, GY . The latter function is derived bysolving the adjoint problem given byL∗ (GY ) = −R∂GY
1 ∂ 2 GY
+ µY GY = δ(Y − η)δ(T − τ ),−∂τPY ∂η2(10)

147

ANALYTICAL MODELING OF NAPL DISSOLUTION

subject to the conditions that eliminate any unknown terms in Equation (9). Expressions for the transverse functions, GY and GZ , can be adapted from the literature
on Green’s functions for diffusion (Beck et al., 1992).
Table I lists longitudinal (GL ) and transversal (GT ) Green’s functions for equilibrium transport – and the conditions on G that they satisfy – that can be used
to derive solutions for a variety of boundary and initial value problems in infinite
and semi-infinite porous media. The superscripts f, s, and t denote a first-, secondand third-type inlet condition for transport in semi-infinite media. The subscript
1 is used to indicate nonequilibrium conditions. Note that the transverse Green’s
function contains the nonequilibrium parameter β, which should be set to unity forequilibrium transport.
2.2.

NONEQUILIBRIUM TRANSPORT

The porous medium is assumed to consist of two macroscopic continua: one inwhich transport is always an equilibrium process and the other where access to immobile water or nonequilibrium sorption sites is governed by a first-order rate equation. The dimensionless model for nonequilibrium transport is given by (Toride etal., 1993; Leij et al., 1993)
βR∂C1
1 ∂ 2 C1
1 ∂ 2 C1 ∂C1
1 ∂ 2 C1
−
+
+
=
+
∂T
PX ∂X 2∂X
PY ∂Y 2PZ ∂Z 2
+ ω(C2 − C1 ) − µ1 C1 ,

(1 − β)R∂C2= ω(C1 − C2 ) − µ2 C2 ,
∂T

(11)
(12)

where C1 is the dimensionless solution for the part of the aqueous phase where
solute transport is described as an equilibrium process, C2 is the concentration
for the rest of the porous medium, β is a coefficient to quantify sorption in the‘equilibrium’ phase relative to total sorption, ω is a mass transfer coefficient andµ1 and µ2 are dimensionless constants for first-order degradation in the mobile
and immobile regions of the medium. The coefficient ω, which quantifies the rateof nonequilibrium exchange to the bulk flow rate, is actually a Damköhler number(Boucher and Alves, 1959). The mathematical equivalency between dimensionless and dimensional transport parameters was illustrated by Toride et al. (1993),
among others.
Before employing the GFM, the following Laplace transform is applied:Z ∞L[C(X, Y, Z, T )] = C̄(X, Y, Z, s) =C(X, Y, Z, T ) exp(−sT ) dT (13)0to Equations (11) and (12). For any arbitrary initial distribution, f (X, Y, Z), thefollowing differential equation may be obtained after some rearrangements, to148

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

Table I. Selected longitudinal and transversal Green’s functions∗
Function
GL (X, T ; ξ, τ )GfL (X, T ; ξ, τ )GtL (X, T ; ξ, τ )GL1 (X, T ; ξ, τ )GfL1 (X, T ; ξ, τ )GtL1 (X, T ; ξ, τ )3(T − τ, τ 0 )

GT (Y, T ; η, τ )GfT (Y, T ; η, τ )GsT (Y, T ; η, τ )ExpressionsRPx×4π(T − τ )#"P [R(X − ξ ) − (T − τ )]2µ (T − τ )− X×exp − X
R
4R(T − τ )sµ (T − τ )RPxexp − X×4π(T − τ )R#"P [R(X − ξ ) − (T − τ )]2−
× exp − X
4R(T − τ )#"P [R(X + ξ ) + (T − τ )]2−exp PX X − X
4R(T − τ ) sµ (T − τ )RPX×exp − XR
4π(T − τ )!
hP [R(X − ξ ) − (T − τ )]2
× exp − X+
4R(T − τ )!P [R(X + ξ ) + (T − τ )]2−+ exp PX X − X
4R(T − τ )#"PR(X + ξ ) + (T − τ )p− X exp(PX X) erfc
2
4R(T − τ )/PX
Z T −ττ003(T − τ, τ )GL (X, ; ξ, 0) dτ 0β0
Z T −τ
0Z T −τ

3(T − τ, τ 0 )GfL (X,τ0; ξ, 0) dτ 0βτ0; ξ, 0) dτ 0β1ω + µX1 0µX τ 0
exp
−
τ ×
βR
βR
βR
s

ω
τ0
×
× δ(T − τ − τ 0 ) +R β(1 − β)(T − τ − τ 0 )


s


ω + µ22ω τ 0 (T − τ − τ 0 ) (T − τ − τ 0 )
×I1 
exp −
R
β(1 − β)(1 − β)Rs!
µ (T − τ )βRPY (Y − η)2
βRPY
exp − Y−
4π(T − τ )R
4(T − τ )sµ (T − τ )βRPYexp − Y×4π(T − τ )R!
!#"
βRPY (Y + η)2
βRPY (Y − η)2− exp −
× exp −
4(T − τ )4(T − τ )sµY (T − τ )
βRPY
exp −
×
4π(T − τ )R!
!#"
βRPY (Y + η)2
βRPY (Y − η)2+ exp −
× exp −
4(T − τ )4(T − τ )03(T − τ, τ 0 )GfL (X,

Conditions∂G= 0, ξ → ±∞∂ξG = 0, ξ = 0∂G= 0, ξ → ∞∂ξ

∂G= 0, ξ = 0, ξ → ∞∂ξ

∂G1= 0, ξ → ±∞∂ξG1 = 0, ξ = 0∂G1= 0, ξ → ∞∂ξ
∂G1= 0, ξ = 0, ξ → ∞∂ξ

∂G= 0, η → ±∞∂ηG1 = 0, η = 0∂G1= 0, η → ∞∂η

∂G1= 0, η = 0, η → ∞
∂η

149

ANALYTICAL MODELING OF NAPL DISSOLUTION

Table I. continued
Function

Expression


GtT (Y, T ; η, τ )Conditions sµ (T − τ )βRPY×exp − YR
4π(T − τ )!
!#"βRPY (Y + η)2
βRPY (Y − η)2+ exp −+
× exp −
4(T − τ )4(T − τ )#"k 2 (T − τ )+ k(Y + η) ×
−k expβRPY
×erfc

βRPY (Y + η) + 2k(T − τ )p
4βRPY (T − τ )

!

kG1 −

∂G1
= 0, η = 0∂η∂G= 0, η → ∞∂η

∗ G = 0 if τ > T .

express the problem in a single dependent variable, C̄1 :1 ∂ 2 C̄1 ∂ C̄1
1 ∂ 2 C̄1
1 ∂ 2 C̄1−+−+
PX ∂X 2∂X
PY ∂Y 2PZ ∂Z 2




aω
ω
− βRs + ω + µ1 −C̄1 + βR +f = 0,
s+b
s+b

(14)

where
a=

ω
,
1−β

(15a)

b=

ω + µ2
.
1−β

(15b)

Appendix A shows how solutions for C1 may be derived for different boundary andinitial value problems from the solution for the equilibrium problem (C). Solutions
for C2 will not be presented since they directly follow from those for C1 (cf. Toride
et al., 1993).
2.2.1. Green’s Function Method
The same procedure is followed as for the equilibrium problem except that the
GFM is applied to the transport equation in the Laplace domain given by Equation(14) rather than the regular time domain. The corresponding differential operator,
L1 , is
1 ∂2
∂
1 ∂2
1 ∂2
+
−
+
−
PX ∂X 2 ∂X PY ∂Y 2 PZ ∂Z 2aω
+ βRs + ω + µ1 −.
s+b

L1 = −

The concentration in the Laplace domain is found by solvingω
L1 (C̄1 ) = βR +
f.
s+b

(16)

(17)

150

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

The following adjoint differential operator can be derived after multiplying thisequality by the Green’s function Ḡ1 and subsequently integrating with respect tothe spatial coordinates:


1 ∂2
∂
1 ∂2
aω
1 ∂2
∗
L1 = −−−+ βRs + ω + µ1 −−.
PX ∂ξ 2 ∂ξ
PY ∂η2 PZ ∂ζ 2
s+b
(18)
The adjoint differential operator on Ḡ1 should now satisfy the problem
L∗1 (Ḡ1 ) = δ(X − ξ, Y − η, Z − ζ ).

(19)

The GFM yields the following expression for the transformed concentration in caseof unspecified initial and boundary conditions:
Z Y2 Z Z2  ξ →X2Ḡ1 ∂ C̄1C̄1 ∂ Ḡ1C̄1 (X, Y, Z, s) =dζ dη +−
− Ḡ1 C̄1
PX ∂ξ
PX ∂ξ
Y1
Z1
ξ →X1
 η→Y2
Z X2 Z Z 2 
Ḡ1 ∂ C̄1C̄1 ∂ Ḡ1+dζ dξ +−PY ∂η
PY ∂η η→Y1
X1
Z1
 ζ →Z2
Z X2 Z Y2 
Ḡ1 ∂ C̄1C̄1 ∂ Ḡ1+dη dξ +
−
PZ ∂ζ
PZ ∂ζ
X1
Y1
ζ →Z1

Z X2 Z Y2 Z Z2 
ω+βR +Ḡ1 f (ξ, η, ζ ) dζ dη dξ.
s+b
X1
Y1
Z1
(20)
The Green’s function, Ḡ1 (X, Y, Z, s; ξ, η, ζ ), can again be written as the productof unidirectional Green’s functions. At this stage only the longitudinal part remains of interest since this part entirely reflects the change in previously reported
equilibrium Green’s functions due to nonequilibrium effects. Rather than solving
Equation (19) for Ḡ1 , and hence ḠX1 (X, s; ξ ), expressions for ḠX1 (X, T ; ξ, τ )will be obtained from the results of Appendix A for a boundary value problem
with input function g(T), and for an initial value problem with initial distributionf(X).
2.2.2. Boundary Value Problem
The concentration in the Laplace domain for a one-dimensional semi-infinite medium which is initially solute free (f = 0) may be deduced from Equation (20):C̄X1 (X, s) =ḠX1 dC̄X1 C̄X1 dḠX1−− ḠX1 C̄X
PX dξ
dξ dξ

 ξ →∞
.
ξ →0

(21)

151

ANALYTICAL MODELING OF NAPL DISSOLUTION

For a first-type condition, the solute concentration at the inlet is given by the
Laplace transform of a time-dependent input concentration, g(T), and a zero-gradientoutlet condition:
C̄X1 (0, s) = ḡ(s),dC̄X1(∞, s) = 0.
dX

(22)

The conditions for the Green’s function, ḠX1 (X, s; ξ ) become apparent after inserting (22) into (21):ḠX1 (X, s; 0) =dḠX1
(X, s; ∞) = 0,
dξ

(23)

where it is assumed that ḠX1 also becomes zero for the infinite outlet condition.
In analogy with the equilibrium solution, the transformed concentration may be
written according to (21) and (A6):ḡ(s) dḠX1(X, s; 0) = ḡ(s)C̄u (X, s ∗ ),PX dξaωs ∗ = βRs + ω + µX1 −,a+bC̄X1 (X, s) =(24)

where, using Duhamel’s theorem, the transformed concentration for an arbitrary
input is related to the corresponding equilibrium concentration for one-dimensionaltransport resulting from a unit input concentration with no degradation, Cu (i.e.,g(T ) = 1 and µX = 0). The degradation coefficient µ1 is partitioned according
to µ1 = µX1 + µY + µZ , while µ2 only occurs in the longitudinal problem.
The Laplace variable is modified to formulate the nonequilibrium concentrationin terms of the equilibrium solution. From (24) it follows thatZ T
g(τ ) ∂GX1
CX1 (X, T ) =(X, T ; 0, τ ) dτ.
(25)
PX ∂ξ
0
f

The Green’s function, listed as GL1 (X, T ; ξ, τ ) in Table I, can be determinedfrom inspection of the solution for CX1 (X, T ) given by (A.10) using the spatial
derivative of equilibrium Green’s functions. Note that solutions for the nonequilibrium problem are written in terms of the Dirac delta function, δ( ), and the zeroand first-order modified Bessel functions, I0 and I1 .
The transformed boundary conditions for a third-type condition are
C̄X1 (0, s) −1 dC̄X1(0, s) = ḡ(s),
PX dX

dC̄X1
(∞, s) = 0.
dX

(26)

The following conditions need to be imposed on ḠX1 (X, s; ξ ) to obtain a uniquesolution according to Equation (21):
dḠX1
dḠX1
(X, s; 0) =(X, s; ∞) = 0.
dξ
dξ

(27)

152

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

The solution in the Laplace domain, again written according to the GFM, is based
upon the corresponding solution to the equilibrium problem:
C̄X1 (X, s) = ḡ(s)ḠX1 (X, s; 0) = ḡ(s)C̄u (X, s ∗ ).

(28)

The concentration in the regular time domain becomesZ T
CX1 (X, T ) =g(τ )GX1 (X, T ; 0, τ ) dτ.

(29)

0

Since the boundary condition does not affect the ‘nonequilibrium’ part of the solution, the Green’s function follows from the nonequilibrium solution and is given as
GtL1 (X, T ; ξ, τ ) in Table I.
2.2.3. Initial Value Problem
By judiciously specifying the two boundary conditions for ḠX1 (X, s; ξ ) to cancel
unknown terms, the concentration in the Laplace domain for a one-dimensionalsemi-infinite medium with an arbitrary initial distribution and a zero solute concentration or flux at the inlet (g = 0), is given according to Equation (20)
as

Z ∞
ω
C̄X1 (X, s) =
βR +
(30)
f (ξ )ḠX1 (X, s; ξ ) dξ.
s+b0
Inversion to the real time domain yields
Z ∞CX1 (X, T ) =f (ξ ) βRGX1 (X, T ; ξ, 0) +
0

+

Z T


ω exp[−b(T − τ )] GX1 (X, τ, ξ, 0) dτ dξ.

(31)

0

The Green’s functions were already derived for the boundary value problem.
In case of a first-type inlet condition, the transformed boundary conditions
are
dC̄X1
(∞, s) = 0.
(32)
dX
To arrive at (30), the Green’s function should satisfy (24). It can be shown that
inserting GfX1 (X, T ; ξ, τ ) yields the solution given by Equation (A.10).
For a third-type condition, the transformed boundary conditions for a zero flux
input areC̄X1 (0, s) = 0,

C̄X1 (0, s) −1 dC̄X1= 0,
PX1 dX

dC̄X1(∞, s) = 0
dX

(33)

with corresponding Green’s function conditions (27). The solution given by (A.10),
including the equilibrium solution for a third-type inlet condition, is obtained byusing GtX1 (X, T ; ξ, τ ) of Table I.

ANALYTICAL MODELING OF NAPL DISSOLUTION

153

3. Persistent Solute Sources
Consider solute movement from a planar source as sketched in Figure 1. This typeof problem arises when a contaminant pool is located at the top (LNAPL source)
or bottom (DNAPL source) of an aquifer. The longitudinal X-domain is infinite,
the horizontal transversal Z-coordinate is infinite while the vertical Y-coordinate issemi-infinite, although the analysis can be readily applied to a finite impermeable
layer (cf. Leijet al., 1999) or other geometries. Three different types of conditions
were applied at the NAPL pool, which has a negligible thickness.
First- and second-type conditions are used for a prescribed concentration and adiffusive mass flux, respectively, while a third-type scenario may be applicable forrate-limited diffusion from the source. Problems with mixed boundary conditions,such as a first- or third-type condition for the source and a homogeneous secondtype condition elsewhere at the surface, are more complicated and will not be
considered here.
The general initial and boundary conditions (i.e., excluding the source condition) are
C(X, Y, Z, 0) = 0 (−∞ < X < ∞, 0 < Y < ∞, −∞ < Z < ∞),

(34)

∂C
(±∞, Y, Z, T ) = 0
∂X

(35)

(0 < Y < ∞, −∞ < Z < ∞),∂C(X, ∞, Z, T ) = 0 (−∞ < X < ∞, −∞ < Z < ∞),
∂Y

(36)

∂C
(X, Y, ±∞, T ) = 0 (−∞ < X < ∞, 0 < Y < ∞).
∂Z

(37)

An intermediate result for the concentration is obtained by substituting the above
conditions in (6):Figure 1. Schematic of a planar solute source.

154

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTENC(X, Y, Z, T ) =Z T Z ∞ Z ∞ 
0−∞
Z ∞Z ∞
0

C ∂G
−
− GC
PX ∂ξ

 ξ →∞dζ dη−
ξ →−∞

C ∂G
+
−∞ −∞ PY ∂η η→∞

 #
G ∂C
C ∂G
+
dζ dξ −−PY ∂η
PY ∂η η=0)
Z ∞Z ∞C ∂G ∞−dη dξ dτ.
PZ ∂ζ ζ →−∞
−∞ 0
−

(38)

The exact formulation of the concentration will be determined in the following
by specifying appropriate conditions on G, written as the product GX (X, T ; ξ, τ )GY (Y, T ; η, τ ) GZ (Z, T ; ζ, τ ), depending on the conditions for the source boundary. In all cases the following conditions apply to G:

3.1.

∂G
(X, Y, Z, T ; ±∞, η, ζ, τ ) = 0,
∂ξ

(39)

∂G
(X, Y, Z, T ; ξ, ∞, ζ, τ ) = 0,
∂η

(40)

∂G
(X, Y, Z, T ; ξ, η, ±∞, τ ) = 0.
∂ζ

(41)

FIRST- TYPE CONDITION

The solute concentration is prescribed for the surface (Y = 0) according tog(X, Z, T ), (X1 < X < X2 , Z1 < Z < Z2 ),C(X, 0, Z, T ) =0, otherwise.

(42)

Equation (38) implies that the following additional condition be placed on theGreen’s function:
G(X, Y, Z, T ; ξ, 0, ζ, τ ) = 0.

(43)

All of the above conditions on G are met by selecting GL , GfT , and GT from
Table I for GX , GY , and GZ . After using the conditions and functions in (38),integration with respect to ξ and ζ is possible due to the simple form of GL and GT
if the prescribed concentration only depends on time. Remaining integrations withtime, for this and other solutions, were done numerically using Gauss–Chebyshev

ANALYTICAL MODELING OF NAPL DISSOLUTION

155

quadrature (Carnahan et al., 1969). The resulting solute concentration isC(X, Y, Z, T ) =Z T Z X2 Z Z 2
0

X1

Z1

g(τ )
GX (X, T ; ξ, τ ) ×
PY∂GY(Y, T ; 0, τ ) GZ (Z, T ; ζ, τ ) dζ dξ dτ
∂η
r


Z T
g(T − τ ) Y RPYµτ
RPY Y 2
=
exp −
−
×
4
τ 4π τ
R
4τ
0





R(X − X2 ) − τR(X − X1 ) − τ× erfc− erfc×√√
4Rτ/PX
4Rτ/PX
"
!#
r


RPZ
RPZ
× erfc(Z − Z2 ) − erfc
(Z − Z1 ) dτ.
4τ
4τ
×

(44)
Two obvious disadvantages of this formulation are that this condition imposes a
sharp discontinuity at the borders of the pool while it is typically not possible tospecify g(T) a priori. Furthermore, an upward diffusive mass flux may occur at
Y = 0, which seems unrealistic away from the pool.
For nonequilibrium transport, the conditions on the solute concentration, C1 andthe Green’s function are essentially the same as for the equilibrium problem, exceptthat GX needs to be replaced by GX1 . The concentration now becomesC1 (X, Y, Z, T ) =Z T Z X2 Z Z 2
0

X1

Z1

g(τ )
GX1 (X, T ; ξ, τ ) ×
PY∂GY(Y, T ; 0, τ )GZ (Z, T ; ζ, τ ) dζ dξ dτ
∂η
( r


Z T
g(T − τ ) Y βRPY
βRPY Y 2
=
exp −
×
4
τ
4π τ
4τ
0
"
!
!#
r
r
βRPZ
βRPZ
× erfc×(Z − Z2 ) − erfc
(Z − Z1 )
4τ
4τ




Z τ
βR(X − X2 ) − τ 0ω + µ1 0
×
erfc
−
exp −
√
τ
βR
4βRτ 0 /PX
0βR(X − X1 ) − τ 0
−erfc
×
√
4βRτ 0 /PX
×156

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTENs"× δ(τ − τ 0 ) +
s
2ω
×
R

ω
R

τ0
I1
β(1 − β)(τ − τ 0 )

)
!
#

τ 0 (τ − τ 0 )ω + µ2dτ 0 dτ.
exp −(τ − τ 0 )β(1 − β)(1 − β)R
(45)

3.2.

SECOND - TYPE CONDITION

Transport from the surface is in this case described as a diffusive process; theconcentration gradient can presumably be specified across the entire surface ofthe porous medium while the concentration at the surface is unknown. The surface
condition is written as∂C−0(X, Z, T ), (X1 < X < X2 , Z1 < Z < Z2 ),(46)(X, 0, Z, T ) =0, otherwise.
∂Y
The diffusive flux 0 may be defined in terms of dimensional parameters using the
Sherwood number (Sh) according to K ∗ cs L/(De c0 ) = Sh cs /c0 . In this expression,k ∗ is a mass transfer parameter (LT−1 ), cs is the solubility concentration (ML−3 ),L is a characteristic length (L), c0 is a reference concentration (ML−3 ) and De isan effective diffusion coefficient (L2 T−1 ) (cf. Chrysikopoulos, 1995). Generally, it
will be difficult to prescribe the flux, 0, since it will normally be controlled by thea priori unknown NAPL concentration of the aqueous phase. If an excessively highflux is specified, the solution may predict concentrations exceeding the solubilityof the NAPL. As will be discussed shortly, the diffusive flux can be constrainedwith a third-type condition.
To eliminate the unknown concentrations in the solution according to Equation(37), the remaining condition on the Green’s function should be∂G(X, Y, Z, T ; ξ, 0, ζ, τ ) = 0∂η(−∞ < X < ∞, −∞ < Z < ∞).

(47)

The Green’s function is the product of GX , GY , and GZ given by GL , GsT , and GT
of Table 1. The concentration for the equilibrium problem – for a flux that only
depends on time – is given by
Z T Z X2 Z Z 20(τ )C(X, Y, Z, T ) =GX (X, T ; ξ, τ ) ×
PY
X1
Z1
0
×GY (Y, T ; 0, τ )GZ (Z, T ; ζ, τ ) dζ dξ dτZ
1 T 0(T − τ )RPY Y 2
µτ
=
√
−
×
exp −
4 0
R
4τ
π PY τ/R

157





R(X − X2 ) − τR(X − X1 ) − τ× erfc− erfc×√√
4Rτ/PX
4Rτ/PX
"
!
!#
r
r
RPZ
RPZ
× erfc(Z − Z2 ) − erfc
(Z − Z1 ) dτ.
4τ
4τ

ANALYTICAL MODELING OF NAPL DISSOLUTION(48)
The solution for the nonequilibrium problem according to the GFM is obtained byreplacing GX in equilibrium solution (48) by GX1 .
3.3.

THIRD - TYPE CONDITION

The need for a third-type condition arises in case of rate-limited dissolution fromthe source. The surface condition is then formulated as∂C(X, 0, Z, T ) = k(X, Z, T )[C(X, 0, Z, T ) − g(X, Z, T )]
∂Y
with

c /c , (X1 < X < X2 , Z1 < Z < Z2 ),g(X, Z, T ) = s 0(49)0, otherwise,
where k is a dimensionless transfer coefficient, which is equivalent to k ∗ L/De ,and cs denotes an upper limit for the (dimensional) solute concentration (e.g., asdetermined by solubility). To simplify the mathematical analysis, it is assumed thatk is constant. Away from the source, a homogeneous second-type condition would
be more accurate since there is no mass transfer in this case (k = 0). The additional
condition for the Green’s function is given by


∂G
kG −
= 0.
(50)
∂η η=0The concentration may now be written as
Z T Z X2 Z Z 2kg(τ )
C(X, Y, Z, T ) =GX (X, T ; ξ ; τ ) ×
PY
X1
Z1
0
×GY (Y, T ; 0, τ )GZ (Z, T ; ζ, τ ) dζ dξ dτ
"s
!

Z T
 µτ 
kg(T − τ )R
RPY Y 2
=
exp −
exp −
−
4
R
π PY τ
4τ
0
 2



k
k τ
RPY Y + 2kτ
−exp+ kY erfc
√
×
PY
RPY
4RPY τR(X − X2 ) − τR(X − X1 ) − τ× erfc− erfc×√√
4Rτ/PY
4Rτ/PX

158

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN"
× erfc

r

!
!#
r
RPZ
RPZ
dτ.
(Z − Z2 ) − erfc
(Z − Z1 )
4τ
4τ
(51)

with GY equal to GtT as defined in Table I. The solution for the nonequilibrium
problem follows directly from Equation (51) by substituting GX1 for GX .

4. Examples
The previously derived solutions will be used to illustrate the effect of the surface condition on the solute distribution, and to explore differences between theequilibrium and nonequilibrium transport models. Parameter values were selected
for illustrative purposes, they reflect laboratory conditions. In all examples it isassumed that there is neither solute retardation (R = 1) nor degradation.

4.1.

SURFACE CONDITION

Consider transport from a rectangular source at the surface Y = 0 with X1 = Z1 =−1 and X2 = Z2 = 1 with PX = 20 (e.g., a pore-water velocity of 10 cm/d anddispersivity of 0.5 cm) and PY = PZ = 200. Figure 2 shows steady-state contoursof the dimensionless solute concentration in the YX-plane at Z = 0 resulting fromadvection and dispersion as calculated with Equations (45), (48), and (51). These
results were obtained with transient Green’s functions, which tend to be easier toapply than steady Green’s functions (Beck et al., 1992). For a first-type condition,
the surface concentration, g, is equal to unity at the source and zero otherwise (Figure 2(a)). The contours for a second-type condition are plotted in Figure 2(b). The
parameter 0 was set to 3.3 at the source – this value was based on a dimensional
concentration gradient of 0.33 mg/cm4 that was normalized using a characteristic
length of 10 cm and a concentration of 1 mg/cm3 – and zero for the rest of thesurface. Contours for a third-type surface condition are given in Figure 2(c) fork = 3.3.
The effects of the boundary conditions can be clearly observed. The first-typecondition, with a zero surface concentration away from the source, is probably notrealistic (Figure 2(a)). The contours obtained with the second-type condition revealthat the highest concentration for a particular longitudinal distance occurs at thesurface; the solute does not move as deep as for the first-type condition. Intermediate surface concentrations were obtained for a third-type condition (Figure 2(c)).
Less solute moves from the pool into the medium for rate-limited diffusion (Figure 2(c)) than for a fixed concentration gradient (Figure 2(b)) or instantaneous
dissolution (Figure 2(a)).

ANALYTICAL MODELING OF NAPL DISSOLUTION

159

Figure 2. Steady-state contours of the dimensionless concentration (C = c/c0 ) in the
XY-plane at Z = 0 as the result of equilibrium transport from a planar source with
X1 = Z1 = −1, X2 = Z2 = 1, PX = 20 and PY = PZ = 200: (a) first-type condition
with g(T ) = 1, (b) second-type condition with 0 = 3.3, and (c) third-type condition due torate-limited dissolution with g(T ) = 1 and k = 3.3.

4.2.

NONEQUILIBRIUM TRANSPORT

Figure 3 shows contours of the relative concentration as a function of depth, Y,
and dimensionless time, T, at X = 10 and Z = 0 for a first-, second-, and thirdtype condition. The development of the concentration profile over time is shownin Figures 3(a–c) for the same conditions as used for Figures 2(a–c). Additionally,corresponding profiles were determined for nonequilibrium transport (Figures 3(d–f). These only pertain to the ‘equilibrium’ concentration. The solute concentration
is calculated according to Equations (45), and modifications of (48) and (51), respectively, using β = 0.5, ω = 1.0, and µ1 = µ2 = 0. Nonequilibrium conditions
may decrease the total amount of solute entering the medium and increase thetime required to reach the eventual peak solute concentration. Furthermore, steeper
concentration gradients may occur because transport is confined to the mobile
zone as compared to transport described by the conventional ADE. Both factors
increase the opportunity for solute spreading, including downward movement, andthe maximum concentrations will consequently be lower.
Figure 4 shows breakthrough curves for C1 at X = 2, Y = 0.1 and Z = 0
for a third-type condition assuming five different degrees of nonequilibrium. The
curves pertain to the same scenario as the contours shown in Figure 3. For a firsttype condition, the highest maximum concentration is obtained for the equilibrium

160

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

Figure 3. Relative concentration (C or C1 ) contours in the YT-plane at X = 2 and Z = 0
resulting from a planar source with X1 = Z1 = −1, X2 = Z2 = 1, PX = 20 and
PY = PZ = 200 according to the equilibrium model (a,b and c) and nonequilibrium modelwith β = 0.5 and ω = 1.0 (d,e,f): (a) and (d) first-type condition with g( T ) = 1, (b) and(e) second-type condition with 0 = 3.3, (c) and (f) third-type condition with g (T ) = 1 andk = 3.3.

scenario (Figure 4(a)). The breakthrough curve is apparently more sensitive to thevalue of the transfer parameter, ω, than the partition coefficient, β. The higher
maximum for ω = 0.1 compared to ω = 1 may be due to the fact that less solute
has yet moved to the nonequilibrium region in the former case. If a third-typecondition is used (Figure 4(b)), the maximum concentration is mostly affected by βwhile ω determines how soon this maximum is reached. The earliest breakthroughoccurs in both cases for the most pronounced case of nonequilibrium (β = 0.25
and ω = 0.1).
5. Summary and Conclusions
The first part of this paper contains a brief review of the use of the GFM to analytically model solute transport as described by the ADE for uniform and steady water
flow. The GFM has not been widely employed for this purpose despite some of its
advantages. Very concise expressions can be written for the solute concentration.
The GFM is well suited to be employed for different boundary and initial conditions. Solutions for multidimensional problems may be obtained by multiplying
the Green’s functions for separate one-dimensional problems. Of particular interest

ANALYTICAL MODELING OF NAPL DISSOLUTION

161

Figure 4. Relative concentration as a function of time (T) at X = 2, Y = 0.1 and Z = 0
resulting from a planar source with X1 = Z1 = −1, X2 = Z2 = 1, PX = 20 and
PY = PZ = 200 predicted for equilibrium conditions and four sets of parameters fornonequilibrium conditions with the following conditions at the source boundary: (a) first-typecondition (g(T ) = c/c0 = 1) and (b) third-type condition (g(T ) = 1 and k = 3.3).

162

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

are the solutions for nonequilibrium transport, which follow in a relatively straightforward manner from those for equilibrium transport. Table I lists longitudinal andtransversal Green’s functions that can be used to formulate the solution for a wide
variety of transport problems in infinite and semi-infinite.
In the second part of this study the GFM was applied to solve transport froma persistent contamination source. Transport took place from a rectangular source
located at the bottom or top of an aquifer during one-dimensional equilibrium andnonequilibrium flow for a first-, second-, and third-type condition at the source.
The first-type and, particularly, the second-type conditions have the disadvantagethat the theoretical solution cannot account for changes in the aqueous phase concentration, which will affect the actual solute flux from the pool. The third-typecondition appears more plausible, but leads to errors away from the pool in case ofa homogenous third-type condition.
The first example demonstrated the surface condition for transport from a planarsource. Solute contours were shown in the XY- and YT-planes. The first-type condition unrealistically mandates that the surface concentration be zero away fromthe source; the maximum concentration was relatively low. For a second-type condition the problem is that the gradient at the source needs to be specified a priori.
A relatively high maximum concentration was predicted, but it occurred at thesurface and possible contamination did not extend far into the medium. The second
example compared equilibrium and nonequilibrium transport from a planar source
using a first-, second-, and third-type condition. Nonequilibrium conditions tended
to lower the maximum concentration but enhanced the spreading of the solute overa larger area. Breakthrough curves were shown for a first-and third-type condition
involving different degrees of nonequilibrium. These curves exhibited the feature that nonequilibrium promotes earlier breakthrough, as is well known for onedimensional transport. However, the maximum concentration may be considerably
lower due to reduced mass transfer into the medium and increased opportunity forsolute spreading.
Appendix A. Derivation of Solutions for Nonequilibrium Transport
The analytical solution procedure will be illustrated for longitudinal transport,
which accounts for all nonequilibrium effects. Solutions for the complete problem
can be readily obtained with the product rule using appropriate transversal Green’s
functions for the diffusion problem (Beck et al., 1992). The temporal Laplace transform of the one-dimensional ADE for equilibrium transport yields the following
ordinary differential equation for an arbitrary initial distribution, f:
1 d2 C̄
dC̄
−− (Rs + µ)C̄ + f = 0.
2
PX dX
dX

(A.1)

The Laplace transform of the corresponding nonequilibrium Equations (11) and(12) can be combined to yield the following transport problem in terms of one

ANALYTICAL MODELING OF NAPL DISSOLUTION

163

dependent variable (Toride et al., 1993):

1 d2 C̄1 dC̄1aωω−
− βRs + ω + µ1 −C̄1 + βR +f = 0.
PX dX 2
dX
s+b
s+b
(A.2)
This equation has a similar form as the equilibrium problem given by (A.1); thisfeature facilitates the derivation of solutions for C1 (X, T ) as shown below.
A.1.

BOUNDARY VALUE PROBLEM

The influent concentration is given by g(T) while the initial concentration f (X) =0. If the solution to the equilibrium problem for g(T ) = 1 and µ = 0 is given
by Cu (X, T ), the solutions of the equilibrium and nonequilibrium problems in the
Laplace domain are related according toaωC̄1 (X, s) = ḡ(s)C̄u X, βRs + ω + µ1 −= ḡ(s)C̄u (X, s ∗ ), (A.3)
s+bwhere s ∗ is a modified Laplace variable obtained from s without spatially or temporally dependent parameters (cf. de Smedt and Wierenga, 1979; Walker, 1987).
The inversion with respect to time is based on the property that the inverse of
the iterated Laplace transform of a function is equal to the generalized convolution integral of the function (Sneddon, 1995). The Laplace transform is applied
sequentially; in this case, t is replaced by either t1 or t2 in the expression for theconcentration. The corresponding transformation variable s is split into s1 and s2 .
The sequential inversion of a function, u, which is transformed in this manner may
be written as−1
h(t1 , t2 ) = L−1
t1 {Lt2 [ū(s1 , s2 )]}.

(A.4)

The function in the regular time domain is given by the convolution integral
Z tu(x, t) =h(τ, t − τ )dτ.
(A.5)
0

Application of this theorem and using the properties of the Laplace transform(Spiegel, 1965) yields

 

aω
−1
−1
h(t1 , t2 ) = Lt1 Lt2 C̄u X, βRs2 + ω + µ1 −
s +b

 


1
aω
t2
−1
−1
= Lt1 exp − ω + µ1 −Lt2 [C̄u (X, βs2 )]s1 + b βR




1
t2
t2
=
Cu X,
exp −(ω + µ1 )
×
βR
βR
βR
r
r




aωt2
aωt2
× δ(t1 ) +
I1 2
(A.6)
t1 exp(−bt1 ) .
βRt1
βR

164

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

From (A.6) the following concentration can be written down:




τ0
1
τ0
C1 (X, T ) =g(T − τ )Cu X,
exp −(ω + µ1 )
×
βR
βR
0
0 βR
s
s
"
!
0
0 (τ − τ 0 )aωτaωτ× δ(τ − τ 0 ) +I1 2×βR(τ − τ 0 )βR
× exp[−b(τ − τ 0 )] dτ 0 ,
(A.7)
Z T

Z τ

where δ( ) is the Dirac delta function and I1 is the first-order modified Bessel
function.

A.2.

INITIAL VALUE PROBLEM

The initial solute concentration is given by f(X), while the influent concentration,g(T), is equal to zero. This problem has not been studied as widely as the boundary
value problem, but inspection of the ordinary differential equations given by (A.1)
and (A.2) suggests that C̄1 (X, s) can be written in terms of the solution of thetransformed equilibrium problem without degradation according to

 

ω
aω
βR +
C̄ X, βRs + ω + µ1 −
s+b
s+b


ω
= βR +C̄(X, s ∗ ).
s+b

C̄1 (X, s) =(A.8)

Again, C1 (X, T ) is obtained with the iterated Laplace transformation. It may beverified that (cf. Spiegel, 1965; Walker, 1987)


 

ω
aω
−1
Lt 2
βR +
C̄ X, βRs2 + ω + µ1 −s1 + bs1 + b


ω
= L−1
βR +
×
t1
s1 + b
 



aω
t2
−1
exp − ω + µ1 −Lt2 [C̄(X, βRs2 )]s1 + b βR




t2
t2
{δ(t1 ) + exp(−bt1 ) ×= C X,
exp −(ω + µ1 )
βR
βR
r
r



 r

ω
aωt2
aωt2
aωt2
×
.
(A.9)
I1 2
I0 2
t1 +
t1
βR
βR
βRt1
βR

h(t1 , t2 ) = L−1
t1

ANALYTICAL MODELING OF NAPL DISSOLUTION165

The solution in terms of the equilibrium concentration therefore becomes



Z T 
τ
(ω + µ1 )T
C1 (X, T ) =C X,
δ(T − τ ) exp −
+
βR
βR
0


ω
ω + µ1
ω + µ2+exp −τ−(T − τ ) ×βR
βR(1 − β)Rs"!
2ω τ (T − τ )× I0+
R β(1 − β)ss!#)2ω τ (T − τ )
βτ
+
I1
dτ.
(A.10)(1 − β)(T − τ )R β(1 − β)with I0 as the zero-order modified Bessel function.
References
Beck, J. V., Cole, K. D., Haji-Sheikh, A. and Litkouhi, B.: 1992, Heat Conduction Using Green’s
Functions, Hemisphere Publ., Washington, DC.
Boucher, D. F. and Alves, G. E.: 1959, Dimensionless numbers, Chem. Engng Progr. 55, 55–64.
Bradford, S. A., Abriola, L. M. and Rathfelder, K.: 1998, Flow and entrapment of dense nonaqueous
phase liquids in physically and chemically heterogeneous aquifer formations, Adv. Water. Resour.
22, 117–132.
Carslaw, H. S., Jaeger, J. C.: 1959, Conduction of Heat in Solids, Clarendon Press, Oxford.
Carnahan, B., Luther, H. A. and Wilkes, J. O.: 1969, Applied Numerical Models, Wiley, New York.
Chrysikopoulos, C. V.: 1995, Three-dimensional analytical models of contaminant transport from
nonaqueous phase liquid pool dissolution in saturated subsurface formations, Water Resour. Res.
31, 1137–1145.
Chrysikopoulos, C. V., Voudrias, E. A. and Fyrillas, M. M.: 1994, Modeling of contaminant transport from dissolution of nonaqueous phase liquid pools in saturated porous media, Transport in
Porous Media 16, 125–145.
Cleary, R. W. and Adrian, D. D.: 1973, Analytical solution of the convective–dispersive equation for
cation transport in soils, Soil Sci. Soc. Am. Proc. 37, 197–199.
Coats, K. H. and Smith, B. D.: 1964, Dead-end pore volume and dispersion in porous media, Soc.
Petrol. Engng J. 4, 73–84.
De Smedt, F. and Wierenga, P. J.: 1979, A generalized solution for solute flow in soils with mobile
and immobile water, Water Resour. Res. 15, 1137–1141.
Ellsworth, T. R. and Butters, G. L.: 1993, Three-dimensional analytical solutions to the advection–
dispersion equation in arbitrary Cartesian coordinates, Water Resour. Res.29, 3215–3225.
Galya, D. P.: 1987, A horizontal plane source model for ground-water transport, Ground Water 25,
733–739.
Van Genuchten, M. Th. and Alves, W. J.: 1982, Analytical solutions of the one-dimensional
convective–dispersive solute transport equation, Tech. Bull., U.S. Dep. Agric., 1661.
Van Genuchten, M. Th. and Wierenga, P. J.: 1976, Mass transfer studies in sorbing porous media. I.
Analytical solutions, Soil Sci. Soc. Am. J. 53, 1303–1310.
Greenberg, M. D.: 1971, Application of Green’s Functions in Science and Engineering, Prentice
Hall, Englewood Cliffs, NJ.
Güven, O., Molz, F. J. and Melville, J. G.: 1984, An analysis of dispersion in a stratified aquifer,
Water Resour. Res. 20, 1337–1354.

166

FEIKE J. LEIJ AND MARTINUS TH. VAN GENUCHTEN

Holman, H.-Y. N. and Javandel, I.: 1996, Evaluation of transient dissolution of slightly water-solublecompounds from a light nonaqueous phase liquid pool, Water Resour. Res. 32, 915–923.
Leij, F. J., Toride, N. and van Genuchten, M. Th.: 1993, Analytical solutions for nonequilibrium
solute transport in three-dimensional porous media, J. Hydrol. 151, 193–228.
Leij, F. J., Priesack, E. and Schaap, M. G.: 1999, Transport from persistent solute sources modeledwith Green’s functions, J. Contam. Hydrol. (in press).
Prakash, A.: 1984, Ground-water contamination due to transient sources of pollution, J. Hydraul.
Engng 110, 1642–1658.
Renardy, M. and Rogers, R. C.: 1993. An Introduction to Partial Differential Equations, SpringerVerlag, New York.
Roach, G. F.: 1982, Green’s Functions, Cambridge Univ. Press.
Shan, C. and Javandel, I.: 1997, Analytical solutions for solute transport in a vertical aquifer section,
J. Contam. Hydrol. 27, 63–82.
Sneddon, I. H.: 1995, Fourier Transforms, Dover, New York.
Spiegel, M. R.: 1965, Theory and Problems of Laplace Transforms, Schaum’s Outline Ser., McGrawHill, New York.
Stakgold, I.: 1979, Green’s Functions and Boundary Value Problems, Wiley, New York.
Toride, N., Leij, F. J. and Van Genuchten, M. Th.: 1993, A comprehensive set of analytical solutions for nonequilibrium solute transport with first-order decay and zero-order production, Water
Resour. Res. 29, 2167–2182.
Walker, G. R.: 1987, Solution to a class of coupled linear partial differential equations, IMA J. Appl.
Math. 38, 35–48.
Yeh, G. T.: 1981, AT123D: Analytical transient one-, two-, and three-dimensional simulation ofwaste transport in an aquifer system, Rep. ORNL-5602. Oak Ridge National Laboratory, Oak
Ridge, TN.
Yeh, G. T. and Tsai, Y.-J.: 1976, Analytical three-dimensional transient modeling of effluent
discharges, Water Resour. Res. 12, 533–540.