On the Use and Error of Approximation

in the Domenico (1987) Solution

Michael R. West 1, Bernard H. Kueper2, and Michael J. Ungs3

Abstract
A mathematical solution for solute transport in a three-dimensional porous medium with a patch source under
steady-state, uniform ground water flow conditions was developed by Domenico (1987). The solution derivation
strategy used an approximate approach to solve the boundary value problem, resulting in a nonexact solution.
Variations of the Domenico (1987) solution are incorporated into the software programs BIOSCREEN and
BIOCHLOR, which are frequently used to evaluate subsurface contaminant transport problems. This article
mathematically elucidates the error in the approximation and presents simulations that compare different versions
of the Domenico (1987) solution to an exact analytical solution to demonstrate the potential error inherent in the
approximate expressions. Results suggest that the accuracy of the approximate solutions is highly variable and
dependent on the selection of input parameters. For solute transport in a medium-grained sand aquifer, the Dome­
nico (1987) solution underpredicts solute concentrations along the centerline of the plume by as much as 80% de­
pending on the case of interest. Increasing the dispersivity, time, or dimensionality of the system leads to increased
error. Because more accurate exact analytical solutions exist, we suggest that the Domenico (1987) solution, and
its predecessor and successor approximate solutions, need not be employed as the basis for screening tools at
contaminated sites.

Introduction integral transforms). Although the configuration of the

Analytical solutions are frequently used as screening domain and source zone geometry is often simplified to
tools to evaluate contaminant fate and transport in the facilitate a tractable solution, if a correct mathematical
subsurface environment. Analytical and semianalytical approach is employed, the solution will be an exact ana­
solutions are exact solutions to a boundary value problem lytical or semianalytical solution to the specified bound­
that comprises a governing differential equation and asso­ ary value problem.
ciated boundary and initial conditions. Many approaches The governing equation for three-dimensional (3D)
have been used in the published literature to derive exact solute transport in saturated porous media under uniform
analytical and semi analytical solutions (e.g., separation steady-state flow subject to advection, dispersion, sorp­
of variables, Laplace integral transforms, and Fourier tion, and first-order decay is (Bear 1979):

lCorresponding author: Department of Civil Engineering,

ac vac a2 c
-a
z2 + ).c=o
+-­
Queen's University, Kingston, Ontario, Canada K7L 3N6; (613) at R ax R
533-2130; fax (613) 533-2128; mwest@civil.queensu.ca
(1)
2Department of Civil Engineering, Queen's University, King­
ston, Ontario, Canada K7L 3N6; (613) 533-6834; fax (613) 533­
2128; kueper@civil.queensu.ca where C C(x,y,z,t) is the concentration of the solute
-­
3Tetra Tech Inc., 3746 Mt. Diablo Boulevard, Suite 300, [M/L3], y is the horizontal transverse spatial coordinate
Lafayette, CA 94549-3681; (510) 283-3771; fax (510) 283-0780.
[L], z is the vertical transverse spatial coordinate [L],
Received April 2006, accepted September 2006.
Copyright © 2007 The Author(s)
x is the longitudinal spatial coordinate [Ll, t is time [T], v
Journal compilation © 2007 National Ground Water Association. is the average linear ground water velocity unidirectional
doi: 10. 1111/j. 1745-6584.2006.00280.x in x [Ur], R is the retardation coefficient, Dx (= CtxVx) is
126 Vol. 45, No.2-GROUND WATER-March-April 2007 (pages 126-135)
the coefficient of longitudinal dispersion [L2/T], Dy The appeal of the Domenico (1987) solution is its util-
(¼ ayvx) is the coefficient of horizontal transverse disper- ity. Since its derivation and resulting expression are rela-
sion [L2/T], Dz (¼ azvx) is the coefficient of vertical trans- tively simple, modifications to boundary conditions and
verse dispersion [L2/T], a is the coefficient of dispersivity, the governing equations (e.g., sorption or multiple species
and k is the aqueous phase decay constant [1/T]. In decay) are implemented with relative ease when compared
Equation 1, both aqueous and sorbed phases undergo to some corresponding exact analytical solutions. In addi-
first-order decay. Note that the boundary value problem tion, it can be incorporated in spreadsheet programs with
assumes that diffusion has a negligible contribution to little effort. Despite its frequent use, little information is
transport and is ignored in the definition of D. available about its inherent error. Domenico and Schwartz
Several exact analytical solutions to Equation 1 have (1998, 381) presented a limited comparison of the Dome-
been presented in the published literature. Cleary and nico (1987) solution to the exact analytical solution by Leij
Ungs (1978) derived a series of exact analytic solutions, et al. (1991) and noted ‘‘some error of approximation.’’ It is
including a patch shaped source zone in a 3D domain that somewhat surprising that the first investigations of error
was semi-infinite in x and finite in y and z. The boundary (beyond the original authors) were reported in research
value problem specified no flux boundaries at the hori- literature only recently (Guyonnet and Neville 2004;
zontal and vertical limits of the domain, and the solution West and Kueper 2004; Falta et al. 2005).
comprises two summation series. Wexler (1992) also It is the experience of the authors that the Domenico
derived an exact analytical solution to the same boundary (1987) solution and its predecessors and successors have
value problem as Cleary and Ungs (1978) but considered been adopted in some regulatory settings and readily used
a domain that was infinite in y and z. The solution by by the contaminant hydrogeology community as not only
Wexler (1992) was an extension of an earlier exact ana- a screening tool but also a decision-making tool for quan-
lytical solution by Sagar (1982), who solved Equation 1 tifying monitored natural attenuation and carrying out
for a conservative solute. The solutions by Sagar (1982) risk analyses. For example, the use of BIOSCREEN and/
and Wexler (1992) comprise an integral that must be or BIOCHLOR at contaminated sites is recommended by
numerically evaluated. Huyakorn et al. (1987), Leij et al. the U.S. EPA (1998, 2005), the United States Air Force
(1991), and Batu (1996) also derived exact analytical sol- Environmental Restoration Program (1998), and the
utions to Equation 1 but used more sophisticated source United States Naval Facilities Engineering Command
boundary conditions than those of the foregoing authors (2001). It has also been relied upon in the research litera-
(e.g., Gaussian distributions, multiple sources of variable ture. Huntley and Beckett (2002) and Falta et al. (2005),
concentrations). Leij et al. (1991) solve Equation 1 for the for example, applied the Domenico (1987) solution in
case of decay of the aqueous phase only. various studies of plume and source zone behavior. While
An alternate approach to the aforementioned bound- Falta et al. (2005) conducted a limited assessment of the
ary value problem solutions was developed by Domenico inherent error for their particular work, no such effort was
(1987), who extended other earlier, less complex sol- undertaken by Huntley and Beckett (2002).
utions proposed by Domenico and Robbins (1985) and It could be argued that the error of approximation
Domenico and Palciauskas (1982). Variations of the arising from the Domenico (1987) solution is irrelevant
Domenico (1987) solution were incorporated into the com- given the uncertainties surrounding issues such as source
puter software programs BIOCHLOR (U.S. EPA 2000, zone delineation and aquifer heterogeneity. However, the
2002) and BIOSCREEN (U.S. EPA 1996), which are use of the Domenico (1987) solution only exacerbates the
available to users through the U.S. EPA. These versions uncertainty with approximate mathematics. The ‘‘uncer-
of the Domenico (1987) solution consider linear isother- tainty’’ due to mathematical approximation can be elimi-
mal sorption and a variety of patch source boundary con- nated by choice, a luxury that cannot always be extended
ditions. Unfortunately, a common misconception among to complex subsurface domains. In addition, mathemati-
many users of BIOCHLOR, BIOSCREEN, and other cal error does not fall within the definition of scientific
versions of the Domenico (1987) solution is that these uncertainty. Uncertainty should only be attributed to phe-
programs use an exact analytical solution. nomena that are uncontrollable, either due to a low state
The Domenico (1987) solution is an approximate in scientific knowledge or a natural complexity that is
solution to Equation 1 and not an exact analytical solu- unobservable. The errors in the Domenico (1987) solution
tion. The solution approach used by Domenico (1987) do not satisfy this definition of uncertainty. In summary,
entailed the assembly of different analytical solutions into why employ an inexact mathematical model when exact
one equation. In particular, a one-dimensional (1D) exact versions are available?
analytical solution for advection and longitudinal disper- In this paper, the approximate nature of the Domenico
sion in a semiinfinite domain (e.g., Bear 1979) was used (1987) solution, and its predecessor and successor sol-
to replace the longitudinal component of a modified exact utions, is investigated both mathematically and with
analytical solution for 3D diffusion (Crank 1975). In modeling. A review of the derivations and mathematics
a preceding article (Domenico and Robbins 1985) that used in these approximate solutions is first presented to
discussed the approach in detail, the authors identified elucidate errors due to solution technique, truncation, and
their solution methodology as an ‘‘extended pulse approx- parameter substitution. To our knowledge, such an analy-
imation’’ and not a rigorous superposition model. The sis of the Domenico (1987) solution, and associated sol-
intent of the authors was to derive an approximate solu- utions, has not been published. Several simulations are
tion suitable for parameter estimation of field data. then conducted to compare the output generated from the
M.R. West et al. GROUND WATER 45, no. 2: 126–135 127
approximate solutions to the exact analytical solutions by for a patch source located between –Y/2
y
Y/2 and
Cleary and Ungs (1978) and Wexler (1992). A full sensi- –Z/2
z
Z/2, where Y is the source width and Z is the
tivity analysis is not the objective of this article. Guyonnet source thickness. The Domenico and Robbins (1985)
and Neville (2004) recently conducted a thorough dimen- technique approximates the 3D transient solution to
sionless comparison between the Domenico (1987) solu- Equation 2 using the product of 1D solutions in each
tion and an exact analytic solution by Wexler (1992) in 3D direction (Domenico and Robbins 1985):
domains. It is our opinion that the use of dimensionless pa-
rameters by Guyonnet and Neville (2004) masks the errors Cðx; y; z; tÞ ’ Cðx; tÞ 3 Cðy; tÞ 3 Cðz; tÞ ð3Þ
inherent in the BIOSCREEN and BIOCHLOR models.
Our approach is pragmatic, clearly illustrating the poten- The methodology reflected in Equation 3 is not
tial for error arising from the use of the Domenico (1987) mathematically rigorous and introduces errors of approxi-
solution in both two-dimensional (2D) and 3D contami- mation relative to exact techniques. Furthermore, Equa-
nant transport applications of practical interest. tion 3 is not a solution to the original mass transport
equation and, therefore, does not conserve mass. This
Approximations in Solution Formulation approach reduces to solving Equation 2 using the follow-
ing three separate partial differential equations:
The historical progression of the various forms of the
Domenico (1987) solution and the different approxi-
mations in solution formulation are discussed in this sec- @C @C @2C
¼ 2v 1 Dx 2 ; 0
x
N ð4Þ
tion and highlighted in Table 1. The approximations are @t @x @x
multiple and varied, including incorrect mathematics,
truncations of exact analytic solutions, and parameter @C @2C
substitutions. ¼ Dy 2 ; 2N
y
N ð5Þ
@t @y
Throughout this section, comparisons are made
between the approximate Domenico (1987) solution and
other exact analytical solutions. In all cases, the compar- @C @2C
isons reference the same boundary value problem. The ob- ¼ Dz 2 ; 2N
z
N ð6Þ
@t @z
jective here is to examine the impact of mathematical
techniques used to derive the Domenico (1987) solution For a patch source, the exact analytic solution to
(and other related solutions), not the underlying boundary Equations 5 and 6 is given by Crank (1975), which, for
value problem assumptions or formulations. Equation 5 as an example, is:
Approximate Expression Development " ! !#
Co y 1 Y=2 y 2 Y=2
The approximate technique employed by Domenico Cðy; tÞ ¼ erf pffiffiffiffiffiffiffi 2 erf pffiffiffiffiffiffiffi ð7Þ
(1987) was developed earlier by Domenico and Robbins 2 2 Dy t 2 Dy t
(1985), who approximated a solution for the case of advec-
tive-dispersive transport in the absence of sorption and where Co is the concentration of the patch source. The
decay. The corresponding governing equation is given by: exact analytical solution to Equation 4 is the Ogata and
Banks (1961) solution:
0
x
N
@C @C @2C @2C @2C     vx   
¼ 2v 1Dx 2 1Dy 2 1Dz 2 ; 2N
y
N Co x 2 vt x 1 vt
@t @x @x @y @z
2N
z
N Cðx;tÞ ¼ erfc pffiffiffiffiffiffiffi 1 exp erfc pffiffiffiffiffiffiffi
2 2 Dx t Dx 2 Dx t
ð2Þ ð8Þ

Table 1
Summary of Solutions

C(x,y,z,t) w
Advection Dispersion Decay Decay C(x,t) 3 Truncation of
Reference (1D) (3D) (aqueous) (sorbed) Sorption C(y,t) 3 C(z,t) C(x,t)

Approximate expressions used in this study

Domenico and Robbins (1985) d d d d

Domenico (1987) d d d d d

U.S. EPA (1996)–BIOSCREEN d d d d d d

U.S. EPA (2000)–BIOCHLOR d d d d d

Guyonnet and Neville (2004) d d d d d d d

Exact analytic solution used in this study

Wexler (1992) d d d d d

128 M.R. West et al. GROUND WATER 45, no. 2: 126–135

 Domenico and Robbins (1985) truncated Equation 8 Equation 3 was employed, which is an approach that
in their work, thus omitting the product of the exp and reduces to solving Equation 10 using the following three
erfc functions. Some discussion on the inherent error partial differential equations:
due to this truncation is found in Bear (1979) and Mar-
@C @C @2C
tin-Hayden and Robbins (1997). When Equation 3 is em- ¼ 2v 1 Dx 2 2 kC; 0
x
N ð11Þ
@t @x @x
ployed, the final truncated solution after substituting t ¼
x/v in the transverse dispersion terms is (Domenico and
Robbins 1985): @C @2C
¼ Dy 2 ; 2N
y
N ð12Þ
  @t @y
Co x 2 vt
Cðx;y;z;tÞ ’ erfc pffiffiffiffiffiffiffi
8 2 Dx t
2 0 1 0 13
@C @2C
6 By 1 Y=2 C By 2 Y=2 C7 ¼ Dz 2 ; 2N
z
N ð13Þ
@t @z
36 B ffiC
4erf @ rffiffiffiffiffiffiffiffi A 2 erf B ffi C7
@ rffiffiffiffiffiffiffiffi
x x A5
2 Dy 2 Dy Again, the exact analytical solution to Equations 12
v v
2 0 1 0 13 and 13 is given by Crank (1975) as presented in Equation
7. The exact analytical solution to Equation 11 is given by
6 B z 1 Z=2 C B z 2 Z=2 C7 Bear (1979):
36 B ffiC
4erf @ rffiffiffiffiffiffiffiffi A 2 erf B ffi C7
@ rffiffiffiffiffiffiffiffi
x x A5
2 Dz
v
2 Dz
v

The substitution of t ¼ x/v warrants some mention.

ð9Þ Cðx; tÞ ¼
Co
2 f exp

0
(
xv
2Dx
" rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#)
12 11 2

4kDx
4kDx
v
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1
Domenico and Palciauskas (1982) proposed that if time
Bx 2 vt 1 1 v2 C
(t) is approximately equal to x/vc, where vc is the contami- 3 erfcB
@ pffiffiffiffiffiffiffi C
A
nant velocity, then the contaminant will spread from Y to 2 Dx t
Y 1 (Dyt)1/2 and Z to Z 1 (Dzt)1/2. While this substitution (
" rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#)
may have some conceptual utility for cases with Dx ¼ 0, xv 4kDx
Dy > 0, and Dz > 0 (i.e., longitudinal plug flow), it is 1 exp 11 11 2
2Dx v
unclear what influence on accuracy interpreting plug flow 0 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1

g
has for cases with Dx > 0. The substitution was re- 4kDx
addressed in a subsequent paper as well. Domenico and Bx 1 vt 1 1 v2 C
3 erfcB
@ pffiffiffiffiffiffiffi C
A ð14Þ
Robbins (1985) described time (t) as ‘‘running time.’’ 2 Dt x
They continued to state that ‘‘reinterpreting this time as x/
v for a moving coordinate system . has the effect of
maintaining the original source dimensions at x ¼ 0 so where Equation 14 is the nontruncated (full) 1D analytic
that the condition C ffi Co is maintained at x ¼ 0 for solution to Equation 11. Both the truncated and non-
t > 0.’’ This interpretation of the substitution is conceptu- truncated solutions are presented in Bear (1979), where the
ally more awkward. A rigorous mathematical solution to a truncated solution omits the second product of exp and
boundary value problem would require no such parameter erfc. Domenico (1987) employed the truncated solution.
substitution to maintain the source boundary condition. Using Equation 3 and the substitution t ¼ x/vc (with vc ¼ v),
The application of Equation 3 presents another con- the completed Domenico (1987) solution is:
flict that is evident in Equation 9. The product of the sepa- ( " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#)
rate solutions for Equations 4 through 6 yields Co3 , not Co. Co xv 4kDx
This issue is not resolved by either Domenico and Robbins Cðx; y; z; tÞ ’ exp 12 11 2
8 2Dx v
(1985) or Domenico (1987). However, Domenico and 0 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1
Schwartz (1998) provide a slightly different formulation of 4kDx
Equation 3 that addresses this conflict, where C(x,y,z,t) is Bx 2 vt 1 1 v2 C
3 erfcB
@ pffiffiffiffiffiffiffi C
A
actually treated as the relative concentration C(x,y,z,t)/Co, 2 Dx t
and is equal to the product of relative solutions for x, y, and
z; thus, Co3 is replaced with Co. 2 0 1 0 13
Similar to Domenico and Robbins (1985), Domenico 6 By 1 Y=2 C By 2 Y=2 C7
36 B ffiC
4erf @ rffiffiffiffiffiffiffiffi A 2 erf B ffi C7
@ rffiffiffiffiffiffiffiffi
(1987) used Equation 3 to approximate a solution for the x x A5
case of advective-dispersive transport subject to aqueous 2 Dy 2 Dy
v v
phase decay in the absence of sorption. The correspond- 2 0 1 0 13
ing governing equation is given by:
6 Bz 1 Z=2 C Bz 2 Z=2 C7
36 B ffiC
4erf @ rffiffiffiffiffiffiffiffi A 2 erf B ffi C7
@ rffiffiffiffiffiffiffiffi
@C @C @ C 2
@ C 2
@ C 2 0
x
N
x x A5
¼ 2v 1 Dx 2 1 Dy 2 1 Dz 2 2 kC; 2N
y
N 2 Dz 2 Dz
@t @x @x @y @z v v
2N
z
N
(10) ð15Þ

M.R. West et al. GROUND WATER 45, no. 2: 126–135 129

 An interesting issue arises in Equation 15 with the The solutions to Equations 17 and 18 can be
substitution of t ¼ x/vc. Since vc is the contaminant veloc- derived using the techniques described by Crank (1975),
ity for this substitution, one can rationalize vc ¼ v/R (for resulting in:
Dx ¼ 0), but the treatment of decay (k) is more awkward. 2 0 1 0 13
Decay in the transverse direction has been excluded in all
Co 6 B 1 Y=2 C
6erf By r
By 2 Y=2 C7
derivations to date and thus introduces an additional error Cðy; tÞ ¼ ffiffiffiffiffiffiffiffi C 2 erf B
@ rffiffiffiffiffiffiffiffi
C7 ð20Þ
of approximation to Domenico (1987) type solutions that 2 4 @ Dy A Dy A5
2 t 2 t
incorporate decay. Furthermore, the application of Equa- R R
tion 3 reduces the 3D decay stated in Equation 10 to 1D
decay in only the x direction (e.g., Bear 1979). 2 0 1 0 13
Co 6 B
6erf Bz 1 Z=2 C Bz 2 Z=2 C7
Adoption of Domenico and Robbins (1985) and Cðz; tÞ ¼ 4 @ rffiffiffiffiffiffiffiffi C
A 2 erf B
@ rD ffiffiffiffiffiffiffiffi C 7
A5 ð21Þ
2 Dz z
Domenico (1987) in subsequent studies 2 t 2 t
R R
The Domenico and Robbins (1985) technique of em-
ploying Equation 3 was implemented in BIOCHLOR (U.S.
EPA 2000, 2002), BIOSCREEN (U.S. EPA 1996), and Employing Equation 3 yields the final approximate
Guyonnet and Neville (2004) to approximate 1D advection nontruncated solution:
with 3D dispersion, including linear isothermal sorption

f
and first-order decay (see Table 1); both Equation 3 and the ( " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#)
truncation of Equations 8 and 14 were implemented in the Co xv 4kRDx
latter two solutions. The original Domenico and Robbins Cðx; y; z; tÞ ¼ exp 12 11
8 2Dx v2
(1985) and Domenico (1987) solutions differ from those
in BIOSCREEN, BIOCHLOR, and Guyonnet and Neville 0 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1
(2004) in that the latter solutions incorporate sorption. vt 4kRDx
Bx 2 11 C
BIOCHLOR and BIOSCREEN consider decay in only B R rffiffiffiffiffiffiffiffi v2 C
3 erfcB C
the aqueous phase, while Guyonnet and Neville (2004) @ Dx A
2 t
examine decay in both the aqueous and sorbed phases. R
( " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#)
The governing equation describing advection, dispersion, xv 4kRDx
linear isothermal sorption, and first-order decay in both the 1 exp 11 11
2Dx v2
aqueous and sorbed phases is given by Equation 1. If Equation 0 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1

g
3 is employed, Equation 1 is reduced to the following three vt 4kRDx
partial differential equations that must be solved individually: Bx 1 11 C
B R 2 C
3 erfcB rffiffiffiffiffiffiffiffi v C
@C v @C Dx @ 2 C @ Dx A
¼ 2 1 2 kC; 0
x
N ð16Þ 2 t
@t R @x R @x2 R
2 0 1 0 13
@C Dy @ 2 C
¼ ; 2N
y
N ð17Þ 6 By 1 Y=2 C By 2 Y=2 C7
@t R @y2
36 B
4erf @ rDffiffiffiffiffiffiffiffi C B
A 2 erf @ rDffiffiffiffiffiffiffiffi C
A5
7
y y
@C Dz @ 2 C 2 t 2 t
¼ ; 2N
z
N ð18Þ R R
@t R @z2 2 0 1 0 13
The exact analytical solution to Equation 16 is given 6 Bz 1 Z=2 C Bz 2 Z=2 C7
by Bear (1979): 36 B
4erf @ rDffiffiffiffiffiffiffiffi C B ffi C7
A 2 erf @ rffiffiffiffiffiffiffi
z D z A5

f
2 t 2 t
( " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#) R R
Co xv 4kRDx
Cðx; tÞ ¼ exp 12 11 ð22Þ
2 2Dx v2
0 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Equation 22 is the nontruncated form of the approx-
vt 4kRDx imate solution to Equation 1 using the Domenico and
Bx 2 11 C
B R 2 C
3 erfcB rffiffiffiffiffiffiffiffi v C Robbins (1985) technique (i.e., Equation 3). If the sub-
@ Dx A stitution of t ¼ x/vc (with vc ¼ v/R) is made in the
2 t
R erf terms, Equation 22 becomes the solution used in
( " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#) BIOCHLOR.
xv 4kRDx
1 exp 11 11 The exact analytical solution to Equation 1 was orig-
2Dx v2 inally derived by Cleary and Ungs (1978) for a finite
0 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 transverse domain. By considering an infinite transverse

g
vt 4kRDx
Bx 1 11 C domain, Sagar (1982) derived an exact solution to Equa-
B R rffiffiffiffiffiffiffiffi v2 C
3 erfcB C ð19Þ tion 1 for the case of a conservative solute. Based on the
@ Dx A
2 t work by Sagar (1982), Wexler (1992) derived an exact
R analytical solution to Equation 1. The exact solution pre-
sented by Wexler (1992) can be used to further illustrate
130 M.R. West et al. GROUND WATER 45, no. 2: 126–135
error in the Domenico (1987) solution and the approach function values close to 1, and the transverse components
by Domenico and Robbins (1985). The exact analytical of Equation 24 can become negligible. In relative terms,
solution to Equation 1 presented by Wexler (1992) is: large source geometries (Y and Z), large values of R, large
rffiffiffiffiffiffiffiffiffi  Z t values of vc, small values of Dy and Dz, and small values
x R xv of x should reduce the influence of the transverse compo-
Cðx; y; z; tÞ ¼ Co exp s23=2
8 pDx 2Dx 0 nents when examining concentrations along the plume
   centerline (y ¼ z ¼ 0). Thus, under these conditions, a
v2 x2 R
3 exp 2 k 1 s2 Domenico (1987) type solution approaches the exact 1D
4Dx R 4Dx s
( " sffiffiffiffiffiffiffiffiffiffiffi # analytical solution of Bear (1979) presented in Equation
Y R 19. As a consequence, a practitioner conducting a suite of
3 erf y 1
2 4Dy s screening calculations with a Domenico (1987) type solu-
" sffiffiffiffiffiffiffiffiffiffiffi #) tion will encounter increasing or decreasing errors of
Y R approximation depending on the choice of parameter.
2 erf y 2
2 4Dy s
( " sffiffiffiffiffiffiffiffiffiffiffi # Quantification of the Error of Approximation
Z R
3 erf z 1 The relative error (r) is evaluated using:
2 4Dz s
" sffiffiffiffiffiffiffiffiffiffiffi #) ðOD 2 OW Þ
Z R rð%Þ ¼ 3 100 ð25Þ
2 erf z 2 ds OW
2 4Dz s
ð23Þ where OD is the output from the approximate Domenico
(1987) solution, and OW is the output from the exact ana-
where s is the dummy variable of integration. The lytical solutions by Cleary and Ungs (1978) and Wexler
Domenico and Robbins (1985) technique (i.e., Equation 3) (1992); For all cases, the solutions by Cleary and Ungs
essentially removes the transverse error functions (erf) (1978) and Wexler (1992) produced the same results. For
from the integration and treats the s as a constant t, giving: convenience, reference will be made to only Wexler
rffiffiffiffiffiffiffiffiffi  Z t (1992) hereafter. The solution by Wexler (1992) (Equa-
x R xv tion 23) requires numerical integration techniques to
Cðx; y; z; tÞ ¼ Co exp s23=2
8 pDx 2Dx 0 evaluate the integral. This need for numerical integration
   prevents analytical analysis of the relative error between
v2 x2 R
3 exp 2 k 1 s2 ds the exact analytical solution by Wexler (1992) and the
4Dx R 4Dx s
( "  s ffiffiffiffiffiffiffiffiffiffi # approximate solutions presented by Domenico and Rob-
Y R bins (1987), Domenico (1987), U.S. EPA (1996, 2000,
3 erf y 1
2 4Dy t 2002), and Guyonnet and Neville (2004). Thus, a compar-
" sffiffiffiffiffiffiffiffiffiffi #) ative evaluation, as presented Equation 25, is necessary.
Y R
2 erf y 2
2 4Dy t
( " sffiffiffiffiffiffiffiffiffiffi # Comparison of Solutions and Input Parameters
Z R Simulations are executed here for both 2D and 3D
3 erf z 1
2 4Dz t porous media with a constant concentration patch source
" sffiffiffiffiffiffiffiffiffiffi #) by comparing Domenico (1987) to Wexler (1992) for
Z R
2 erf z 2 both transient and steady-state conditions. Exact analy-
2 4Dz t tical solutions by Huyakorn et al. (1987), Leij et al.
ð24Þ (1991), and Batu (1996) could also have been used for
this purpose. The 3D cases simulate an aquifer with only
The x component of Equation 24 is actually the solu- a partially penetrating source zone, while the 2D case
tion by Bear (1979) presented in Equation 19; hence, considers a source zone distributed through the full aqui-
once the integral is analytically evaluated, Equation 24 fer thickness. The consideration of dimensionality not
reduces to Equation 22. There is no mathematical basis only highlights the influence of domain configuration on
to justify Equation 3 or the approximation of s ¼ t in the error of approximation in Domenico (1987) but also
Equation 24; it is a matter of convenience. illustrates the differences in prediction for the two cases,
Based on the aforementioned approximations, it is which can be overlooked while screening a site. We
anticipated that source width, source depth, transverse assume that our example site is fairly well characterized,
dispersion, retardation, transverse dimensionality, and such that field solute concentrations, hydraulic gradient,
time (or the substitution of t ¼ x/vc) will contribute to the and geological conditions are known. We assume that the
‘‘error of approximation.’’ In some instances, the con- contaminant is trichloroethene with a source concentra-
ditions where the error is minimized can be surmised tion (Co) of 11 mg/L, and all concentration values are
from the properties of the error function (erf). The erf is assessed along the plume centerline, C(x,0,0,t), for all
an asymptotic exponential integral where erf(N) / 1. simulations. In addition, the plume length (Lp), defined by
Therefore, large values of the argument (i.e., the ratio of the position of the 5 parts per billion contour, is assumed
the numerator to the denominator) will yield error to be 100 m for the purposes of calculating a constant
M.R. West et al. GROUND WATER 45, no. 2: 126–135 131
dispersivity. The parameters are summarized in Table 2.
Note that R ¼ 1 for all simulations as the Domenico Table 3
(1987) solution does not incorporate sorption. The depth Summary of Steady-State Simulations
of the source zone below the water table is 2.5 m.
Run Solution Configuration ax
One would normally conduct a suite of screening
simulations to investigate the influence of various param- 1 Wexler (1992) 2D 10 m
eters (e.g., velocity, source width) on the solution output. 2 Domenico (1987) 2D 10 m
A subset of simulations is conducted herein, where only 3 Wexler (1992) 2D 0.1 L
the dispersivity values are adjusted. A full dimensionless 4 Domenico (1987) 2D 0.1 L
analysis was conducted by Guyonnet and Neville (2004). 5 Wexler (1992) 3D 10 m
As mentioned in the previous section, transverse disper- 6 Domenico (1987) 2D 10 m
sion and dimensionality will potentially introduce signifi- 7 Wexler (1992) 3D 0.1 L
8 Domenico (1987) 3D 0.1 L
cant error due to the mathematical approximations
employed in the Domenico and Robbins (1985) tech-
nique. The three approaches suggested in the BIO-
CHLOR manual (U.S. EPA 2000, 2002) are implemented solution for C(x,t). These simulations are summarized in
to obtain reasonable dispersivities for the subject example Table 5 and use the parameters listed in Table 2 where
site. In the first approach, the longitudinal dispersivity appropriate. Relative error is evaluated along the plume
(ax) is 10% of the Lp; using this approach, a constant ax ¼ centerline while employing dispersivity values calculated
10 m is calculated for all values of x. The second and using ax ¼ 0.1 L.
third approaches recognize the scale dependent nature of
dispersivity (Gelhar et al. 1992). The method of Pickens
and Grisak (1981) is adopted where ax ¼ 0.1 L, where L Steady-State Simulations
is the distance along the domain to the point of reference. Steady-state concentration profiles for solute transport
The third method uses the findings of Xu and Eckstein in 2D and 3D porous media are presented in Figure 1,
(1995) where ax ¼ 0.83(log10L)2.414; in the latter two which compares the output generated from solutions by
approaches, ax increases with x (i.e., L). For this ex- Domenico (1987) and Wexler (1992). The simulations are
ample, the horizontal transverse dispersivity is given by conducted using either a constant ax ¼ 10 m or a variable
ay ¼ 0.1ax and, where appropriate, the vertical transverse ax ¼ 0.1 L. Clearly, the solutions are sensitive to the
dispersivity is given by az ¼ 0.01ax. Both transient choice of dispersivity, with the variable ax producing
and steady-state simulations were executed in the 2D study greater concentrations near the source zone and reduced
by assuming that the vertical transverse dispersivity was concentrations downstream relative to the constant ax.
negligible (i.e., az ¼ 0 and Dz ¼ 0). Descriptions of the Furthermore, the dimensionality of the problem has a
steady-state and transient simulations are provided in significant impact on the concentration profiles, which
Tables 3 and 4, respectively. The relative error between underscores the need for dimensional screening at sites.
solution outputs for each of the simulations is evaluated The relative error between the two solutions is pre-
using Equation 25. sented in Figure 2 for the four combinations of dis-
Finally, a series of transient 3D simulations are con- persivity and dimension. The Domenico (1987) solution
ducted to demonstrate the error due to truncation of the

Table 4
Table 2 Summary of Transient Simulations
Summary of Input Parameters for All Simulations
Distance to
Description Parameter Value
Monitoring Method of
Source concentration Co 11 mg/L Run Case Solution Well (m) Calculating ax
Source width Y 10 m 1 A Wexler (1992) 10 Xu and
Source depth below Z/2 2.5 m Eckstein (1995)
the water table (3D case) 2 A Domenico (1987) 10 Xu and
Source depth below Z/2 2.5 m Eckstein (1995)
the water table (2D case) 3 A Wexler (1992) 100 Xu and
Source decay half-life t(c)1/2 N Eckstein (1995)
Solute plume decay half-life t(k)1/2 5 years 4 A Domenico (1987) 100 Xu and
Hydraulic gradient rh 0.008 Eckstein (1995)
5 B Wexler (1992) 10 Pickens and
Hydraulic conductivity K 1 3 1024 m/s
Grisak (1981)
Longitudinal dispersivity ax Various
6 B Domenico (1987) 10 Pickens and
Transverse horizontal dispersivity ay 0.1ax Grisak (1981)
Transverse vertical dispersivity az 0.01ax 7 B Wexler (1992) 100 Pickens and
Porosity h 0.25 Grisak (1981)
Ground water velocity v 0.277 m/d 8 B Domenico (1987) 100 Pickens and
Solute retardation factor R 1.0 Grisak (1981)

132 M.R. West et al. GROUND WATER 45, no. 2: 126–135

 Table 5
Simulation Notation for Figure 5

Label Description

Full 1 year Relative error between Wexler

σ( )
(1992) and the full (nontruncated)
Domenico (1987) solution 2D αx = 10 m
2D αx = 0.1L
(e.g., BIOCHLOR) for 1 year 3D αx = 10 m
Truncated 1 year Relative error between Wexler 3D αx = 0.1L
(1992) and the truncated
Domenico (1987) solution (i.e.,
Domenico [1987] and
BIOSCREEN) for 1 year 1 10 100 1000
Full 2 years Relative error between Wexler Distance along domain (m)
(1992) and the full (nontruncated)
Domenico (1987) solution Figure 2. Relative error between Wexler (1992) and Dome-
nico (1987) for 2D and 3D steady-state simulations.
(e.g., BIOCHLOR) for 2 years
Truncated 2 years Relative error between Wexler
(1992) and the truncated Domenico
(1987) solution (i.e., Domenico variable ax yields an error of 0.5% near the source zone
[1987] and BIOSCREEN) that continually increases with distance, such that an
for 2 years error of 232% is realized at a distance of approximately
Full 5 years Relative error between Wexler 1000 m.
(1992) and the full (nontruncated)
Domenico (1987) solution
(e.g., BIOCHLOR) for 5 years Transient Simulations
Truncated 5 years Relative error between Wexler (1992) Breakthough curves are presented in Figure 3 for
and the truncated Domenico (1987) monitoring wells located at 10 and 100 m downstream of
solution (i.e., Domenico [1987] and the source. Case A uses dispersivities calculated from Xu
BIOSCREEN) for 5 years and Eckstein (1995), while the method of Pickens and
Grisak (1981) is used in case B. The case A ax values for
the 10- and 100-m wells are 0.83 and 4.42 m, respec-
underpredicts the steady-state concentrations through the tively. The ax values of 1 and 10 m were applied to the
majority of the domain, with the exception of the constant 10- and 100-m wells for case B, respectively.
ax for distances close to the source zone. As suggested Figure 3 illustrates that the Domenico (1987) solu-
earlier from mathematical considerations, the error in the tion can underpredict the concentrations at each well for
Domenico (1987) solution is exacerbated by increasing all times. The relative error associated with the under-
the dimensionality of the domain from 2D to 3D. The prediction is presented in Figure 4. The error between
choice of either a constant or variable ax has a significant Wexler (1992) and Domenico (1987) is significant at
impact on the behavior of the error. For the 3D case, the early time for all wells, with values of 230% for the 10-m
constant ax produces relative errors ranging from 2.5% well and 280% for the 100-m well. For late time, the
to 224%, peaking near the source zone. However, the well near the source zone (10 m) exhibits negligible
error. However, the error at the 100 m well achieves
a steady-state relative error of 27% for case A but
12 216% for case B. It can be observed that increasing the
Wexler (1992)
2D αx = 0.1L dispersion in the system increases the error between the
Domenico (1987)
10
2D αx = 0.1L
solutions.
Concentration (mg/L)

Wexler (1992) Wexler (1992)

8 3D αx = 0.1L
2D αx = 10 m
Domenico (1987) Domenico (1987)
Influence of Truncation in C (x,t)
2D αx = 10 m 3D αx = 0.1L
6 Profile curves of relative error between the truncated
and nontruncated Domenico (1987) solutions (see Table 5)
4 and Wexler (1992) are presented in Figure 5. A total of
Wexler (1992)
3D αx = 10 m six simulations were conducted to examine the relative
2 Domenico (1987)
3D αx = 10 m error at t ¼ 1, 2, or 5 years for Equation 25.
The influence of truncating C(x,t) can be observed by
0
1 10 100 1000 comparing the relative error calculations for each simula-
Distance along domain (m) tion. The truncated cases (e.g., Domenico [1987] and
BIOSCREEN) yield additional error of approximation when
Figure 1. Comparison of solutions for steady-state condi- compared to the nontruncated cases (e.g., BIOCHLOR);
tions in 2D and 3D porous media. Simulations are summa-
rized in Table 3. however, all cases produce significant error with increas-
ing x. At a time of 1 year and a distance of 100 m, the
M.R. West et al. GROUND WATER 45, no. 2: 126–135 133
 12 12
A B
10 10 Both solutions
Both solutions
x= 10 m

Concentration (mg/L)
x = 10 m

Concentration (mg/L)
8 8
Wexler (1992)
Wexler (1992) x = 100 m
6 x = 100 m 6

Domenico (1987)
Domenico (1987)
4 x = 100 m 4
x = 100 m

2 2

0 0
0 1 2 3 4 0 1 2 3 4
Time (yrs) Time (yrs)

Figure 3. Comparison of breakthrough curves at x = 10 and 100 m between Wexler (1992) and Domenico (1987). (Case A)
Simulations with dispersivities calculated using Xu and Eckstein (1995). (Case B) Dispersivities calculated using Pickens and
Grisak (1981). Simulations are summarized in Table 4.

nontruncated Domenico (1987) solution produced a rela- behavior and magnitude of the relative error cannot be
tive error of 232%, but the truncated Domenico (1987) predicted a priori. Thus, users of the Domenico (1987)
solution yields a relative error of 241%. The additional solution and other associated versions cannot ascertain
error of approximation diminishes with increasing time as the degree of underprediction without a comparative
the two solutions approach steady-state conditions. study (as conducted herein) or potentially by examining
the type curves of Guyonnet and Neville (2004).
Additional relative error is observed in expressions
Conclusions (i.e., BIOSCREEN and Domenico [1987]) that truncate
The accuracy of the Domenico (1987) solution is the analytical solution for 1D solute transport. When per-
parameter dependent. Depending on the choice of para- forming simulations to evaluate the effects of solution
meters, the magnitude of relative error can range from formulation and truncation, the relative error between the
significant (e.g., 280%) to negligible. When comparing exact analytical solutions by Cleary and Ungs (1978) and
Domenico (1987) to the exact analytical solutions of Wexler (1992) and the various approximate expressions
Cleary and Ungs (1978) and Wexler (1992), the Domeni- (i.e., Domenico [1987]; BIOSCREEN; BIOCHLOR) rang-
co (1987) solution typically underpredicts concentrations. ed from 2% to 280%, depending on the case of interest.
The magnitude of relative error is a function of parameter Guyonnet and Neville (2004) stated that along the
value, time, and dimensionality; the application of dis- plume centerline, and for ground water flow regimes
persivity is particularly awkward as different techniques dominated by advection and mechanical dispersion rather
yield variable error behavior. Most important, the than by molecular diffusion, the discrepancies between

10
Cases A & B 5
0 x = 10 m
0 Full & Trunc - 5 yrs
-10
-5 Full - 2 yrs
-20 Trunc - 2 yrs
-10
-30
-15
Case B
σ( )

σ( )

-40 x = 100 m -20

-50 -25 Full - 1 yr

-60 -30 Trunc - 1 yr

Case A
-70 x = 100 m -35

-80 -40

-90 -45
0 1 2 3 4 0 10 20 30 40 50 60 70 80 90 100
Time (yrs) Distance (m)

Figure 4. Relative error for breakthrough curves at x = 10 Figure 5. Transient 3D comparison of relative error
and 100 m between Wexler (1992) and Domenico (1987). between different approximate expressions and Wexler
Case A denotes simulations with dispersivities calculated (1992) at t = 1, 2, and 5 years. Longitudinal dispersivity
using Xu and Eckstein (1995) while dispersivities calculated was estimated using ax = 0.1 L. Refer to Table 5 for label
using Pickens and Grisak (1981) are presented in case B. descriptions.

134 M.R. West et al. GROUND WATER 45, no. 2: 126–135

the Domenico (1987) solution and the exact Sagar (1982) Guyonnet, D., and C. Neville. 2004. Dimensionless analysis
solution can be considered negligible for practical pur- of two analytical solutions for 3-D solute transport
in groundwater. Journal of Contaminant Hydrology 75,
poses. Based on Figures 1 through 5 of this work, we take
141–153.
a more conservative stance and suggest that in the context Huntley, D., and G. Beckett. 2002. Persistence of LNAPL sour-
of maximum allowable concentrations at contaminated ces: Relationship between risk reduction and LNAPL
sites, the solutions developed using the technique by Do- recovery. Journal of Contaminant Hydrology 59, 3–26.
menico and Robbins (1985) can potentially significantly Huyakorn, P., M. Ungs, L. Mulkey, and E. Sudicky. 1987. A
three-dimensional analytical method for predicting leachate
underpredict solute concentrations along the plume cen-
migration. Ground Water 25, no. 5: 588–598.
terline when compared to other exact analytical solutions Leij, F., T. Skaggs, and M. van Genuchten. 1991. Analytical
for both transient and steady-state conditions in 2D and solutions for solute transport in three-dimensional semi-
3D domains. infinite porous media. Water Resources Research 27, no.
Given the aforementioned magnitude of relative error 10: 2719–2733.
Martin-Hayden, J., and G. Robbins. 1997. Plume distortion and
and the abundance of exact analytical solutions available to
apparent attenuation due to concentration averaging in
practitioners, the Domenico (1987) solution and its associ- monitoring wells. Ground Water 35, no. 2: 339–346.
ated solutions need not be employed as screening tools. We Ogata, A., and R. Banks. 1961. A solution of the differential
fully encourage practitioners to use peer-reviewed pub- equation of longitudinal dispersion in porous media. USGS
lished exact analytic solutions, preferably with open source Prof. Pap. No. 411-A. Reston, Virginia: USGS.
Pickens, J., and G. Grisak. 1981. Scale-dependent dispersion in
code. In some cases, a compiled program with a graphic
a stratified granular aquifer. Water Resources Research 17,
user interface is more desirable. There are many such soft- no. 4: 1191–1211.
ware packages available in industry. While not a compre- Sagar, B. 1982. Dispersion in three dimensions: Approximate
hensive list, some examples of commercial software that analytical solutions. ASCE, Journal of Hydraulics Division
use exact analytical solutions include PRINCE (distributed 108, no. HY1: 47–62.
United States Air Force Environmental Restoration Program.
by Waterloo Hydrogeologic Inc.), ATRANS (S.S. Papado-
1998. Handbook for Remediation of Petroleum-Contami-
pulos Inc.), and SOLUTRANS (Fitts Geosolutions). nated Sites (A Risk Based Strategy). Prepared for Air Force
Center for Environmental Excellence, Technology Transfer
Division.
Acknowledgments United States Naval Facilities Engineering Command. 2001.
Funding for this work was provided by Texas Instru- Guidance for Optimizing Remedial Action Operation
ments Inc. Additional funding was provided by the Ontario (RAO), Special Report SR-2101-ENV. Prepared for
Department of the Navy RAO/LTM Optimization Working
Graduate Scholarship Program, the Natural Sciences and Group.
Engineering Research Council of Canada, and Queen’s U.S. EPA. 2005. Monitored Natural Attenuation of MTBE as
University in the form of student scholarships to the a Risk Management Option at Leaking Underground Stor-
senior author. The authors would also like to acknowl- age Tank Sites, EPA/600/R-04/1790. Washington, DC: U.S.
edge and thank the four anonymous reviewers for their EPA, Office of Research and Development.
U.S. EPA. 2002. BIOCHLOR Natural Attenuation Decision
thorough and insightful comments. Support System. Version 2.2 User’s Manual Addendum.
Washington, DC: U.S. EPA.
U.S. EPA. 2000. BIOCHLOR Natural Attenuation Decision
References Support System. User’s Manual Version 1.0, EPA/600/R-00/
Batu, V. 1996. A generalized three-dimensional analytical solute 008. Washington, DC: US EPA, Office of Research and
transport model for multiple rectangular first-type sources. Development.
Journal of Hydrology 174, 57–82. U.S. EPA. 1998. Technical Protocol for Evaluating Natural
Bear, J. 1979. Hydraulics of Ground Water. New York: McGraw- Attenuation of Chlorinated Solvents in Ground Water,
Hill. EPA/600/R-98/128. Washington, DC: US EPA, Office of
Cleary, R., and M. Ungs. 1978. Report No. 78-WR-15, Water Research and Development.
Resources Program. Princeton, New Jersey: Princeton U.S. EPA. 1996. BIOSCREEN—Natural Attenuation Decision
University. Support System. User’s Manual Version 1.3, EPA/600/
Crank, J. 1975. The Mathematics of Diffusion. New York: R-96/087. Washington, DC: National Risk Management Re-
Oxford University Press. search Laboratory, Office of Research and Development.
Domenico, P.A. 1987. An analytical model for multidimensional West, M., and B. Kueper. 2004. Natural attenuation of solute
transport of a decaying contaminant species. Journal of plumes in bedded fractured rock. In Proceedings of the U.S.
Hydrology 91, 49–58. EPA/NGWA Fractured Rock Conference: State of Science
Domenico, P.A., and F.W. Schwartz. 1998. Physical and Chem- and Measuring Success in Remediation. National Ground
ical Hydrology, 2nd ed. New York: John Wiley & Sons Inc. Water Association, NGWA Press.
Domenico, P.A., and G.A. Robbins. 1985. A new method of con- Wexler, E. 1992. Analytical solutions for one-, two-, and three-
taminant plume analysis. Ground Water 23, no. 4: 476–485. dimensional solute transport in groundwater systems with
Domenico, P.A., and V.V. Palciauskas. 1982. Alterative boundaries uniform flow. In Techniques of Water Resources Inves-
in solid waste management. Ground Water 20, no. 3: 303–311. tigations of the United States Geological Survey, chap. B7,
Falta, R., N. Basu, and R. Rao. 2005. Assessing impacts of par- Book 3, Applications of Hydraulics, 53–55.
tial mass depletion in DNAPL source zones: II. Coupling Xu, M., and Y. Eckstein. 1995. Use of weighted least-squares
source strength functions to plume evolution. Journal of method in evaluation of the relationship between dis-
Contaminant Hydrology 79, 45–66. persivity and field scale. Ground Water 33, no. 6: 905–908.

M.R. West et al. GROUND WATER 45, no. 2: 126–135 135