WATER RESOURCES RESEARCH, VOL. 25, NO.

6, PAGES 1141-1148, JUNE 1989

~,0
Contaminant Transport and Biodegradation
A Numerical Model for Reactive Transport in Porous Media
MICHAEL A. CELIA

Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge

J. SCOTT KINDRED

AssociatesInc., Redmond,Washington

ISMAEL HERRERA

lnstituto de Geofisica, National Autonomous University of Mexico, Mexico City

A new numerical solution procedure is presented for simulation of reactive transport in porous
media. The new procedure, which is referred to as an optimal test function (OTF) method, is
formulated so that it systematically adapts to the changing character of the governing partial
differential equation. Relative importance of diffusion, advection, and reaction are directly incorpo-
rated into the numerical approximation by judicious choice of the test, or weight, function that appears
in the weak form of the equation. Specific algorithms are presented to solve a general class of
nonlinear, multispecies transport equations. This includes a variety of models of subsurface contam-
inant transport with biodegradation.

1. INTRODUCTION The recent work of Herrera [1984, 1985a,b], Herrera et

A proper description of contaminant transport is impor- al. [1985], and Celia et al. [1989a] provides a systematic
tant in many aqueous systems. For subsurface flows, de- framework in which numerical solutions to general second-
scription of the flow physics must often be augmented by order equations can be obtained. The procedure leads to
chemical and/or biological considerations. This generally numerical approximations that automatically change as the
leads to advective-di1Iusive-reactive transport equations. governing partial differential equation changes. This is
When multiple species are present in the aqueous phase, the achieved by choosing test functions, which are present in a
governing equations form a set of partial differential equa- weak form of the original governing equation, to satisfy
tions that are coupled through reaction terms. These equa- specific, rigorously derived criteria. No arbitrary parameters
tions generally need to be solved numerically. are involved. It is the purpose of this paper to develop a
The transport equation is one for which numerical solution solution algorithm, based on this procedure, for systems of
procedures continue to exhibit significant limitations for advective-ditIusive-reactive transport equations. In a com-
certain problems of physical interest. The most widely cited panion paper [Kindred and Celia, this issue] the methodol-
example is the case of advection domination. In this case, ogy is applied to a specific model of transport in biologically
one is usually forced to choose between nonphysical oscil- reactive porous media.
lations or excessive (numerical) diffusion. While a variety of The presentation in this paper begins by exposing the
methods have been developed specifically for advection- underlying numerical algorithm for the case of a single
dominated flows [Leonard, 1979; Hughes and Brooks, 1982; advective-ditIusive-reactive transport equation with con-
Tezduyar and Ganjoo, 1986; Baptista, 1987], only partial stant coefficients. This allows the salient features of the
success can be reported to date. method to be explained in a simple mathematical setting.
A key to providing reliable numerical simulations is rec- The procedure is then extended to the case of spatially
ognition of the changing nature of the governing equation. variable coefficients, and finally to the case of multiple
f
Diffusion domination implies behavior analogous to that equations that are coupled through nonlinear reaction terms.
predicted by model parabolic partial differential equations; Example calculations are presented to demonstrate the
advection domination implies behavior analogous to first- numerical procedure and to provide a direct link to the
order hyperbolic partial differential equations; and reaction biodegradation model that is developed in the companion
domination implies first-order ordinary differential equation paper [Kindred and Celia, this issue].
(in time) behavior. A numerical procedure that fails to
accommodate these disparate behaviors cannot be expected
to produce reliable numerical solutions.

NUMERICAL ALGORITHM-SINGLE SPECIES

TRANSPORTWITH CONSTANT COEFFICIENTS
Copyright 1989 by the American Geophysical Union.
Paper number 89WRO0271. The numericaldevelopmentbeginsby examininga single
0043-1397/89/89WR-00271
$05.00 species, advective-difIusive-reactive transport equation
1141

,
Golder
 1142 CELIA ETAL.: CONTAMINANT TRANSPORT AND BIODEGRADATION, 1

R au/at + v au/ax -D a2u/ax2+ Ku = f(x, t) (1) The first key to the numerical procedure is application of
integration by parts twice to each term of the summation in
O~x~L t>O (4). For each element integration, the first integration by
u(O, t) = go t> 0 parts produces
au/ax(L, t) = gL t> 0
i~j+I(:£xU)W(X) dx
u(x, 0) = uo(x) 0~ x ~ L
defined over the finite spatial domain n = [0, L]. Any
boundary conditions can be accommodated in the numerical
= (Xj+ I a2u
-;;::z -V
ax
au
--Ku
ax
) D
w dx =
au
D-w-
ax
Vuw
algorithm; the conditions given above provide a typical JXj
example. Nomenclature is such that R is a retardation au dw dw
coefficient (dimensionless); V is fluid velocity (L/1); D is a
+ (XI + I I
-D--+
axdx
Vu--Kuw
dx
) dx
diffusion coefficient (L2/1); K is a first-order reaction coeffi- JX} \ /

cient (1/1); f is a given forcing function; and u(x, t) is the Application of an additional integration by parts then yields
dependent scalar of interest, which in this case is a measure
of concentration of a dissolved substance. The forcing L X}+I (:£ u)w dx = [ D -wau -Du --Vuw
dw

x ax dx
function f(x, t) is a source/sink term that may include x) -
reaction terms that are not dependent on the unknown
function u. To begin the numerical development, let the + (Xj+ I d2wDp+Vdw
coefficients R, V, D, and K be assumed constant. The dx -KW)UdX
numerical procedure that is to be applied to (1) consists of JXj
application of the algebraic theory of Herrera [1984, 1985a, = [ D a;
b] in space, thereby producing a semidiscrete system of
au w -Du ~ L
dw -Vuw ] XJ+l + XJ+' (:;E~w)udx (5)

ordinary differential equations in time. This resulting set of XJ Xl

ordinary differential equations is then integrated in time where :;Exis the original spatial operator and :;E':is its formal
using standard methods.
adjoint.
Let (1) be rewritten in the form The second key to the numerical procedure is to recognize
:£xu = D a2u/ax2-V au/ax -Ku = R au/at -f(x, t) (2) that the original integral (on the left side) of (5) can be
replaced by nodal evaluations by choosing w(x) such that
where the operator :£x incorporates both the spatial deriva- :;E*xw = 0 within each [xi' xi+I]' That is to say, proper choice
tives and the reaction term. The weak form of (2) is then
of the test function w(x) effectively concentrates information
formed by multiplying the equation by a weight, or test,
at node points and eliminates interior element integrations.
function, w(x), and integrating over the domain [0, L] In addition, the concentration of information at node points
-L is accomplished in the absenceof any trial function definition
(D a2u/ax2-V au/ax -Ku)w(x) dx (as would be the case in a standard finite element formula-
tion). Both the special choice of test function and the lack of
a trial function distinguish this numerical formulation from
= iL(R au/at -f)w(x) dx (3) standard finite element or Petrov-Galerkin methods.
Consider a choice of w(x) that satisfies :;E*xw= 0 within
Next, let the spatial domain [0, L] be subdivided into E each[xi' xi+I],j=O, 1, ..., £-1. Furthermore, allow w(x) the
subintervals {[XO,XI]' [XI' X2]' ..., [XE-I' XE]}' with Xo = 0, ability to exhibit discontinuities at node points. Given such a
XE = L. These subintervals, or elements, are separated by definition ofw(x), (4)and (5) can be combined and restated in
the E + 1 node points {XO,XI' ..., XE}' If the solution u is terms of known coefficients and unknown nodal values of the
assumed to be at least (:;1continuous, u E (:;1[0,L](that is, u dependentvariable. In light of (5) and the fact that :;E*xw = 0,
and au/ax are continuous functions over Os Xs L), and w(x) (4) can be rewritten as
is at least (:;-1 continuous, wE (:;-1 [0, L] (that is, J~ w(x')
dx' is a continuous function over 0 ~ X ~ L), then the Eil (Xj'
(:£xu)w dx
integral on the left side of (3) can be written equivalently as j = 0 J x)

iL(D a2u/axZ-V au/ax -Ku)w(x) dx

£-1
= j~O i:j+l(D iJ2U/iJx2-V iJu/iJx -Ku)w(x) dx (4)
E-)
dw au
= L D ~ + Vw UJ" -[[Dw ]1-;- j
} uX
The equality (4), which introduces elementwise integration, j= )
is permissible due to the continuity constraints on u and w.
The case of lower continuity on u has been treated by
Herrera [1984, 1985a, b]; for the present development, u E
-Vw ) u + Dw au
-(6)]L
ax 0
1[1will suffice. Also, let any discontinuities that occur in w(x) +[(-D~
be restricted to node points. Within each [Xp Xj+I]' it is In (6), the double bracket denotes a jump operator, which is
assumed that w(x) E 1[2. defined by [[o]]x = lim.,--.o[(o)x+e-(o)x-e]. Since D, V, and
J J J

r
 CELIA ET AL.: CONTAMINANT TRANSPORT AND BIODEGRADATION, 1 1143

ware are all known, the only unknowns in (6) are the 2£ + boundary conditions, 2E+2linear algebraic equations result
2 nodal values {Uj, (iJu/iJx)j}.f~o'These nodal values corre- for the 2E+2 nodal unknowns. Solution of this set of
spond to the function U and its spatial derivative. These equations produces exact nodal values, since no approxima-
values are unique at each node by the continuity assumption tion has been introduced, and the set of test functions is T
U E C1[O,L]. Equation (6) can be written more succinctly as complete [see Herrera, 1984]. Detailed algorithms for the
a simple linear combination of known coefficients and un- case of constant coefficient, ordinary differential equations
known nodal values are presented by Herrera et al. [1985].

i L E au
(:£xu)w dx = 2: Ajuj + Bj a;j
When au/at e~ 0, the time dimension must be included.
Examination of (3) indicates that a spatial integral of the
(6') product of (au/at)(x,t) and w(x) must be evaluated. To
0 j=O perform this integration, au/at is approximated using a
polynomial expansion that involves the nodal values appear-
For the operator :£x of (2), the formal adjoint operator is,,ccording
to (5), ing in the expression on the right side of (6). The natural
choice for the approximation developed above is a piecewise
:£1 = Dd2/dx2+ Vd/dx K (7) cubic Hermite polynomial interpolation in space, since this
involves nodal values of u and au/ax, namely,
The homogeneoussolutionscorrespondingto :£; are given
by E
iJu/iJt= iJu/iJt= L {Vj(t)c/Joj(x)+ (iJV/iJx)J{t)<P\j(x)} (10)
'l'\(x) = exp [(a + /3)x] (8a)
j=O
'l'2(X) = exp [(a -/3)x] (8b) In (10),{Vi' (iJV/iJx)i}f~oare time-dependent nodal values of
where a = -VIW and /3 = (1/2D) (V2 + 4KD)I/2. Equations function and spatial derivative, and <Poi,<P1jare standard
cubic Hermite polynomials (see, for example, Carey and
(8a) and (8b) represent the two fundamental solutions of the
aden [1983, pp. 63-65]). Substitution of expansion (10) into
homogeneous adjoint equation. Any linear combination of
the integral of interest yields
these solutions also satisfies the homogeneous adjoint equa-
tion. The computational procedure proposed herein uses as
test fullctions two nonzero solutions to ':£:W = 0 defined in
each element and formed as linear combinations of solutions
(8a) and (8b). These functions are defined such that they are
nonzero only within the element of interest, and zero in all
other elements. As such, each w(x) E (:-1[0, L]. For any
d
di [ -;;:;- ] (L cf>ij(X)w(X)dx}+
iJUj (1) Jo
(11)
element e, defined by [xi' Xj+1]' the test functions are defined
by Given that cPoj{X),cPv{x),and w(x) are each well defined and
w1(x) = C11 exp [(a + f3)x) + C12 exp [(a -f3)x) known functions, the integrals can be evaluated directly and
(9a) (11) can be written, with inclusion of the coefficient R, as

wj(x) = 0 x < Xj X> Xj+ R ~ax) (II'

wz(X) = C21 exp [(a + f3)x] + C22 exp [(a -f3)x] (9b) Thus the resulting approximation that derives from (11), (6),
and (3) is
J

w~(x)=O X<Xj x>Xj+1 "E { ao-(U)+

d ,B°-d ( -1aUo ) -AoUo-Bo-1 dUo}
L., :I dt J :/ dt ax :I J J axj=O
The constants CII, C12' C21' and C22, are chosen, for
convenience, to satisfy the conditions w~(x) = 1, W~(Xj+J =
0, W2(X) = 0, and W2(Xj+J = 1. Therefore CII = -exp
f(x, t)w(x) dx (12)
[-(a+ f3)xj-2f3(llx)]/{1-exp
[-(a-f3)xj]/{l-exp
[-2f3(llx)]},
[-2f3(llx)]},
C 12 = exp
C21 = -exp
=iL.
[-(a-f3)xj+I-2f3x)/{1-exp [2f3(llx)]}, and C22 = exp An equation of the form of(12) is written for each of the 2£
[-(a-f3)xj+J/{l-exp [2f3(1lx)]), with Ilx = Xj+I-Xj' Any choices of w(x) given in (9). The two boundary conditions
other linear combination of the fundamental solutions (8) required for the second-order equation (I) provide two
would be equally acceptable. Since there are E elements, 2E additional equations, so that 2£+2 equations result for the
equations are produced corresponding to the 2E linearly 2£+2 nodal values. This provides the semidiscrete system of
independent test functions, {w~(x), W2(X)}:=I' The undeter- equations
mined nodal values that appear in (5) are the function u and P .dU/dt -Q .U = F (13)
the spatial derivative au/ax, forming the set of 2E+2 un-
knowns {Uj, (au/ax)j}f~o' Two boundary conditions provide where matrix P contains the ap f3jcoefficients of (11), Q
the two additional equations needed to close the system. contains the coefficients Aj, Bj of (6), and the vector U
For steady state conditions, au/at = 0 and (1)reduces to an contains nodal values of U and au/ax. The structure of the
ordinary differential equation. After evaluation of the right coefficient matrices P and Q is exactly analogous to that of
side forcing term I~x)w(x) dx and imposition of the two standard cubic Hermite collocation. This structure is illus-

Xi<X<Xi+
X'<X<Xj+
 1144 CELIA ET AL.: CONTAMINANT TRANSPORT AND BIODEGRADATION, I

p=

1. Typical five-diagonal matrix structure. Matrix P of (13) is shown. Matrix Q has the same structure. First and
last rows correspond to boundary conditions.

trated in Figure 1. Computational requirements for this functions exhibit an upstream bias. This corresponds to the
optimal test function (OTF) method are very similar to those physically asymmetric process of advective transport. The
for Hermite collocation algorithms. amount of the upstream bias, or "upstream weighting," is
Equation (13) can be solved by any time integrator of directly related to the parameters of the governing equation
choice. A simple scheme is the variably weighted Euler and involves no arbitrary parameters or coefficients. This is
method. in contrast to traditional upstream finite difference or finite
element methods, wherein the degree of upstreamingis
p. (un+l-un)/llt-Q. [8Un+l+(1-8)un]=Fn+6
generally related to an arbitrary parameter. Finally, when
(14a)
reaction dominates, the functions are weighted toward the
or nodes, and spatial coupling is diminished. This is consistent
[(l/~t)P-IJQ] .Un + I = [(II Ilt)P+(l- 8)Q] .Un + Fn+ 6 with the spatial derivatives in the governing equation being
unimportant, and the point process (in space) of the first-
(14b) order reaction being dominant.
where 8 is a weighting parameter, usually taken as 0 :5 8 :5 This automatic shifting of the test functions in accordance
I. Imposition of the initial conditions and subsequentmarch- with the physics inherent in the governing partial differential
ing through time yields the discrete approximation of inter- equation implies a certain optimality of the test functions.
est. Examination of the equality of (5) reinforces the optimality
This development, for both steady state and transient claim in the sense that the differential equation is written in
cases,requires the equation to be second order in space, that
is, D # O. If the equation is purely advective, then solution

X
of the homogeneous adjoint operator equation is no longer
given by (8), and the entire deyelopment must be reformu- A
lated. This is consistent with the formal change from a
"o e e
parabolic to a hyperbolic equation. While the procedure for
0.6
0.6 r W1 w2 J
pure advection is analogous to that presented above, the
present formulation only considers the case of D # O. 0.4
j
0.2
J
J
3. OPTIMAL TEST FUNCTIONS nl
The choice w(x) given by (9) derives from a rigorous and 0.4 (J.tS O.a 1.0

exact mathematical treatment of the spatial operator in the X

governing equation. The test functions reflect all of the
physical processes that are described by the governing
differential operator. Furthermore, the formulas of (9) indi-
cate that the test functions respond automatically to changes
in the coefficients D, Y, and K. Figure 2 illustrates three
different sets offunctions, one each for the casesof diffusion
domination, advection domination, and reaction domination.
When diffusion dominates, the functions approachpiecewise
linear forms (analogous to finite element "hat" functions,
although the current functions are fully discontinuous).
Since diffusion is a symmetric process that has spatial Fig. 2. Typical test functions, plotted over one element for the
coupling, these functions are seen to correspond to the cases of (a) diffusion domination, (b) reaction domination, and (c)
dominant physical process. When advection dominates, the advection domination.

Fig.
 CELIA ET AL.: CONTAMINANT TRANSPORT AND BIODEGRADATION. I 1145

terms of nodal values in a mathematically precise way. equation of the form of (1), where the coefficients are
Furthermore, optimal accuracy occurs for steady state equa- functions of the unknown u, linearization involves evalua-
tions in both one dimension (exact nodal solutions for tion of the coefficients using previous (known) values of the
equation (1) when R = 0) and two dimensions (see Discus- solution. These can be values of previous iterations or values
sion in section 8). The numerical algorithm presented above at previous time steps. In either case, the linearized govern-
is therefore referred to as an OTF method. ing equation is analogous to a variable coefficient problem.
Techniques discussed in the previous section therefore per-
4. TREATMENT OF NONCONSTANT COEFFICIENTS tain.
When the coefficients in (1) exhibit spatial variability, it is Linearization can take one of several forms. For time
generally not possible to exactly satisfy the homogeneous marching problems, nonlinear coefficients can be estimated
adjoint equation, ~:w = O. In this case, an approximation from solutions at previous time steps, allowing for solution
procedure is required to provide good estimates to w, so that at the new time level exclusive of iteration. Conversely,
~:w = O. progressively improved solutions at the new time level can
Two approaches have been developed to treat the variable be obtained through iterative solution procedures. These
coefficient case. The first, which is used herein, involves iterative procedures include simple (Picard) iteration, New-
replacement of the coefficients in (1) with piecewise La- ton-Raphson iteration, and various modifications of each.
grange polynomials. The approximate coefficients have dif- The linearization that is chosen herein for the case of
ferent polynomial definitions in each element and are gener- nonlinear reactive transport uses a noniterative time march-
ally discontinuous across element boundaries. The basic ing algorithm coupled with a second-order projection tech-
idea is that with these approximate coefficients, the adjoint nique for estimating the nonlinear coefficients. Estimates of
equation is simply an ordinary differential equation with nonlinear coefficients are based on solutions at the previous
polynomial coefficients. As such, two independent series two time levels. For a time weighting parameter 8, as used in
solutions can be developed for the homogeneous adjoint (14), the appropriate projected value is given by
equation using standard techniques (see, for example, Hilde-
brand [1976]). Celia and Herrera [1987] have shown that un+8=(1+8)un-8un-l (15)
through judicious choice of interpolation points within each This is a second-order (in time) estimate for Un+B.This very
element,piecewise nth degreeLagrangeinterpolants can simple linearization was chosen because of the relatively
provide O«l1xfn+2) accuracy in the numerical solution. The mild nonlinearities and the very good numerical results that
arguments used to develop this theory are analogous to those were subsequently produced. As stated previously, applica-
used to choose collocation point locations in the collocation tion of the optimal test function method is indifferent to the
finite element method (see, for example, Prenter [1975]). The linearization chosen; OTF is applicable to any linearization
interpolation points in the OTF method, and the collocation
procedure.
points in the collocation finite element method, both turn out
to be the Gauss-Legendre integration points within each
element [Celia and Herrera, 1987]. An alternative approach 6. SETS OF COUPLED EQUATIONS-REACTIVE
for developing test functions uses a local (elementwise) TRANSPORTOF MULTIPLE SPECIES
collocation solution for the homogeneous adjoint equation.
Mathematical description of transport of multiple species
This again yields two independent solutions for w(x) within
in flowing groundwater involves a set of coupled, advective-
each element. Details of this local collocation procedure are
diffusive-reactive transport equations. Depending on the
provided by Herrera [1987].
nature of both the flow system and the chemical or biological
For the current treatment, consider approximation of the
reactions, different forms of equations pertain. Rubin [1983]
coefficients by piecewise constants. That is, the coefficients
provided a detailed overview of various types of chemical
are constant within each element and can change from one
element to the next. Since ~:w = 0 within each element reactions and their concomitant mathematical descriptions.
Kindred and Celia [this issue], in a companion to this paper,
separately, the adjoint equation in any element involves
discuss biologically reactive media and its related mathemat-
constant coefficients. Therefore si~ple analytical solutions
ical description.
can be obtained for the test functions. By the error estimate
The set of transport equations chosen for the current
cited above, this provides O«l1xf) accuracy in the approxi-
presentation has the following form:
mation. If D(x), V(x), and K(x) are replaced by D(x), V(x),
and R(x), with overbarred quantities referring to piecewise 2
iJc) iJc) iJ c)
constant approximations, then the test functions remain as at + V) ~ -D) -a?"+ K)(c), ..., cN, XI, ...,XM)c)
defined in (9) with D, V, and K, interpreted to be the values
within the element of interest. = !i(CI, , CN,X\,
2
5. TREATMENT OF NONLINEAR EQUATIONS aC2 aC2 a C2
-+
at V 2 --D
ax 2 ~ax- + K2(cI ..,....C N XI ., XM)Cz
When the governing partial differential equation is nonlin-
ear, some linearization technique must be applied to allow
for tractable numerical solution. This is true of any numer- =fic ...,CN'Xi,...,XM) (16a)
ical approximation, including the optimal test function
method. The linearization mayor may not involve iteration,
2
depending on the severity of the nonlinearity and the pref- OCN OCN 0 CN
-+ VN--DN-;:-T+ KN(Ct,' " CN, XI, ..., XM)CN
erence of the analyst. ot ox or
For the optimal test function method, a linearization is
performed prior to derivation of the adjoint operator. For an = fN<c\, , CN,X\, ,XM)

~
 1146 CELIA ET AL.: CONTAMINANT TRANSPORT AND BIODEGRADATION,

G\(c\, , CN,X"
iJt

=F,(c" , CN,X\, ,XM)

(16b)

~ + GM(c\, CN,X\, XM)XM

at
= FM(c CN,X1, ,XM)

In (16), Cl' CZ'..., CNare measuresof aqueousconcentrations

of N dissolved species, each of which is hydrodynamically
transported and has the ability to react with other constituents.
The variables Xl, Xz, ..., XM are concentration measuresof
reactive species that are stationary. This set of equations is 0 10 20 30 40 50 60 70 80 90 100
taken to generally correspond to various models of biodegra-
Distance
dation discussedby Kindred and Celia [this issue]. The present
objective is to describe the overall numerical algorithm pro- Fig. 3. Comparisonof numericalsolutionsto analytical solu-
tions. Solid curves indicate OTF numerical solutions; symbols
posed to solve the set (16). Kindred and Celia [this issue] indicateanalyticalsolutionvalues.
provide derivation and physical motivation for the equation
forms as well as additional computational details.
The solution procedure adopted herein utilizes the optimal transport systems. A variety of problems related specificallyto
test function concepts developed in the earlier sections of this biodegradationare solved in the companion paper [Kindred
paper. To proceed from one time step to the next, a simple, and Celia, this issue] using the techniques developed herein.
noniterative linearization is employed. Equation (15)is used to Additional test problems for the specific case of advection
provide required estimatesof the solution for evaluationof the dominated, nonreactive transport are reported by Celia et al.
nonlinear reaction coefficients and forcing functions. The pro- [1989a]. Bouloutas and Celia [1988] provide an in-depth nu-
jection is applied to each Ci' i= 1,2, ..., N, as well as each Xi, merical analysis of a related OTF algorithm for nonreactive
i= 1, 2, ..., M. Let the projected values be denoted with an transport of a single species.
asterisk. The linearized version of (16)is thus written
7.1. Comparison to Analytical Solutions
Equation (1), with constant coefficients, is solved for
, c~,xf,
several different combinations of the coefficients. Three
casesare considered: (1)R = I, K = 0; (2)R = 1, K = 0.01;
= Ji(cf, cN,Xf, XXI) and (3) R = 2, K = 0.01. Each case uses V = 1.0, D = 0.2,
f(x,t) = 0.0, and 6x = 2.0, implying a grid Peclet number of
ax/at + Gj(cf, CN,Xj, 10. Initial conditions are fixed as uo(x) = 0, and boundary
conditions are go = 10, gL = o. Figure 3 shows numerical
=Fj(cf,"',CN,Xf,"',XXI) (17b) results for each case, as well as the respective analytical
solutions. Agreement is excellent in each case. Numerical
These equations are used to moye from the level n to n + 1. mass balance errors are less than 1% for the first two
Given known values c;*,XJ',each Kj, Gj,fj, and Fj are known solutions and approximately 1.5% for the third solution.
functions of space. These spatially variable coefficients are
next approximated by piecewise constants, with each coef- Two-Species Transport
ficient taken as constant within any element. Given these
A second example is presented that involves transport of
piecewise constant values, the adjoint operator is obtained two species, coupled through the reaction terms. One sta-
for each element. Homogeneous solutions of the adjoint tionary component is included, but it is assumed to be
equation are subsequently obtained and used as test function temporally static. One physical analogy for this system is
in the optimal test function algorithm. Each transport equa-
transport of two substrates in the presence of a stationary
tion has its own set of test functions, and each equation is biological population. Dissolved species might be an organic
solvedseparately. contaminant and dissolved oxygen, leading to an aerobic
The equations for X,{X, t) have no spatial derivatives, imply- biodegradation transport problem. Simulations that include
ing no coupling in space. As such, these equations are solved temporal dynamics of the stationary species are included in
for each node in space, with no coupling between nodes. the work by Kindred and Celia [this issue].
Linearization is againbasedon the projection techniqueof(15).
Governing equations for this system take the form
Use of constant values of Gj and Fj over a time step, based on
the values c;, X;, allows (17b)to be solved analytically. iJc)/iJt+ V iJc)/iJx-D iJZclliJxZ
+ K)(c), XJc) = h(cz, Xu

7. EXAMPLE CALCULATIONS
iJc7/iJt+ViJc7/iJx -D iJzc7/iJxz+ Kz(cz, XJcz =!z(c\, XJ
Several example calculations are presented to demonstrate
the viability of OTF methods for advective-diffusive-reactive

+
(17a)
7.2.
XM)X1
XM)Xi
iJX\
xXt)Ci
 CELIA ET AL.: CONTAMINANT TRANSPORT AND BIODEGRADATION, I

DISTANCE(m)
Fig. 4. Numerical solution of (18) and (19) using OTF. Solutions are plotted at time t = 68 days.

Specific functional forms for the nonlinear coefficients are species 2. This error is not a consequence of the time
chosen, based on the analogy to biodegradation, as follows: stepping algorithm, because the solution was insensitive to
time step reductions. Reduction of the space step Ax from 2
K1(CI, XJ = (V~xl/(Kl + cJ)81 (19a)
m to I m reduced the mass balance errors to less than 3% for
Kz(cz, XJ = (V~XI/(K~ + cz))az (19b) species I and less than 1% for species 2. The actual solutions
for Ax = 2 and Ax = 1 are virtually indistinguishable at t =
II (cz, XJ = -KIZ(V;"xI/(K~+ cz))azcz = -KIZKzcz (19c) 68 days, with only a slight change evident at the interacting
concentration fronts (the fronts are displaced by approxi-
h(cl,XJ= -KZI(V:nXI/(Kl + cJ)8lcl = -KZIK1cl (19d) mately 1 m). This appears to be a consequence of the
where physical interpretations are that V:" is the maximum piecewise constant approximation of the spatially variable
uptake rate for species i, K~, the half-saturation constant for (nonlinear) reaction coefficient, which is a relatively poor
species i, Kij is the yield ratio coefficient for species i when
approximation in the vicinity of the sharp concentration
speciesj is limiting, and 8jis equal to 1 if species i is limiting fronts. Implementation of the high-order spatial procedures
the reactions and zero otherwise. Species 1 and 2 are of Celia and Herrera [1987] or reduction of the space step (as
assumed to react in a fixed ratio; this is reflected in the indicated above) will alleviate the mass balance discrepancy.
coefficients Kij. Equations (18) allow either of the species to
While the mass balance approaches unity as Ax decreases,
be limiting. There is no bacterial growth in this equation (GI the numerical solutions for concentrations remain virtually
= 0), so the stationary species has a fixed value X I. Deriva- indistinguishable from those presented in Figure 4.
tion and discussion of these equations is provided in the
companion paper [Kindred and Celia, this issue].
DISCUSSION AND CONCLUSIONS
Parameter values are assigned as follows: V:" = 1.0
days-I, i = 1, 2; K~ = 0.1 mg/L, i = 1,2; KI2 = 2.0, K21 = This paper presents a new numerical method for solution
0.5. Initial conditions are Cj(X, 0) = 3.0, C2(X, 0) = 0.0. of reactive transport problems. The technique derives fromjudicious
Imposed boundary conditions are CI(O, t) = 3.0, C2(0, t) = mathematical manipulation of the weak form of the
10., iJCj/iJx (L, t) = iJC2/iJX(L, t) = 0, where L = 100 m. The governing equations. Special choice of the test functions,
stationary species is fixed at XI = 0.2 mgiL. Flow parame- based on solution of the homogeneous adjoint equation,
ters are V = 1.0 rn/day and D = 0.2 m2/day. The system leads to an optimal approximation technique. The method
described by these parameters corresponds to a step intro- has the desirable property of automatic adjustment of the
duction of species C2 at the left boundary. This propagates approximation to accommodate varying degreesof diffusion,
into the domain by advection and diffusion. Its presence advection, and reaction domination.
instigates uptake of both species 1 and 2. Figure 4 shows For one-dimensional, steady state problems (ordinary
plots of concentration as a function of distance for time t = differential equations with constant coefficients), the method
68 days, for both species 1 and 2, calculated using the OTF produces exact nodal values for the unknown function and,
algorithm with grid spacing f).x = 2 m and time step ilt = 0.2 if desired, its first derivative. Procedures of arbitrarily high
days. The reaction occurs until species I, which is limiting to order, on fixed grids, have been developed for the case of
the left of the invading front, disappears, at which point variable coefficients [Celia and Herrera, 1987; Herrera,
reaction ceases. To the right of the front, species 2 is 1987]. In addition, tensor product test functions applied to
limiting, due to its initial concentration of zero. The grid the two-dimensional Laplace equation lead again to an
Peclet number for this example is 10. This small to interme- optimal approximation. This latter case can be shown to
diate value of Grid Peclet number is easily accommodated by produce a sixth"order approximation using only nine node
the OTF method. This example is expanded in the compan- points [Celia et al., 1989b]. Collatz [1960] has shown that
ion paper of Kindred and Celia [this issue] to include both this is the highest possible order of approximation for a nine
inhibited and uninhibited bacterial growth. point approximation to the Laplace equation. The OTF
Mass balance errors for this two-species example were approachtherefore appearsto be a very promising numerical
approximately 5% for species 1 and approximately 3% for method.

8.
1147
1148 CELIA ET AL.: CONTAMINANT TRANSPORT AND BIODEGRADATION,

Current efforts are focusing on improved treatment of the Herrera, I., Unified approach to numerical methods, I, Green's
time domain as well as extensions to multiple spatial dimen- formula for operators in discontinuous fields, Numer. Methods
Partial Differential Equations, 1(1), 25-44, 1985a.
sions. As indi~ated above, initial efforts in multiple dimen-
Herrera, I., Unified approach to numerical methods, II, Finite
sions using tensor product test functions have produced very elements, boundary elements, and their coupling, Numer. Meth-
encouraging results. We plan to extend these results to the ods Partial Differential Equations, 1(3), 159-186, 1985b.
case of multidimensional reactive transport. We are also Herrer~, I., The algebraic theory approach for ordinary differential
investigating several formulations that apply the OTF con- equations: Highly accurate finite differences, Numer. Methods
Partial Differential Equations, 3(3), 199-218, 1987.
cept to the full space-time equation. The companion paper to Herrera, I., L. Chargoy, and G. Alducin, Unified approach to
this one [Kindred and Celia, this issue] uses the OTF numerical methods, III, Finite differences and ordinary differen-
simulator developed herein as a numerical tool to demon- tial equations, Numer. Methods Partial Differential Equations,
strate behavior of various biologically reactive subsurface 1(4), 241-258, 1985.
Hildebrand, F. B., Advanced Calculus for Applications, 2nd ed.,
transport systems. Prentice-Hall, Englewood Cliffs, N.J., 1976.
Hughes, T. J. R., and A. Brooks, A theoretical framework for
Acknowledgments. This work was supported in part by the Petrov-Gaierkin methods with discontinuous weighting functions:
National Science Foundation under grant 8657419-CESand by the Applications to the streamlin~-upwind procedure, in Finite Ele-
Massachusetts Institute of Technology under Sloan Foundation mentsin Fluids, vol. 4, edited by R. H. Gallagheretal., pp. 47-65,
grant 26950. John Wiley and Sons, New York, 1982.
Kindred, J. S., and M. A. Celia, Contaminant transport and biodeg-
REFERENCES radation, 2, Conceptual model and test simulations, Water Re-
sour. Res., this issue.
Baptista, A. M., Solution of advection-dominatedtransport by Leonard, B. P., A stable and accurate convective modelling proce-
Eulerian-Lagrangianmethodsusing the backward methods of dure based on quadratic upstream interpolation, Compo Meth.
characteristics,Ph.D. thesis,Dep. of Civ. Eng., Mass. Inst. of Appl. Mech. Eng., 19, 55-98, 1979.
Technol.,Cambridge,Mass.,1987. Prenter, P. M., Splines and Variational Methods, 323 pp., John
Bouloutas,E. T., and M. A. Celia, An analysis of a class of Wiley and Sons, San Diego, Calif., 1975.
Petrov-Galerkinand Optimal Test Functions methods,in Proc. Rubin, J., Transport of reacting solutes in porous media: Relation
SeventhInt. Coni ComputationalMethodsin WaterResources, between mathematical nature of problem formulation and chem-
edited by M. A. Celia et al., pp. 15-20,ElsevierScience,New ical nature of reactions, Water Resour. Res., 19(5), 1231-1252,
York, 1988. 1983.
Carey,G. F., andJ. T. Oden,Finite Elements:A SecondCourse, Tezduyar, T. E., and D. K. Ganjoo, Petrov-Gaierkin formulations
volume2, The TexasFinite ElementSeries,Prentice-Hall,Engle- with weighting functions dependent upon spatial and temporal
wood Cliffs, N.J., 1983. discretization: Application to transient convective-diffusion prob-
Celia,M. A., and I. Herrera,Solutionof generalordinarydifferen- lems, Compo Method Appl. Mech. Eng., 59, 49-71, 1986.
tial equations by a unified theory approach,Num. Methods
Partial Differential Equations,3(2), 117-129,1987. M. A. Celia,ParsonsLaboratory, Room48-207,Departmentof
Celia,M. A., I. Herrera,E. T. Bouloutas,andJ. S. Kindred,A new Civil Engineering,MassachusettsInstitute of Technology,Cam-
numerical approachfor the advective-diffusivetransport equa- bridge,MA 02139.
tion, Num. Methods Partial Differential Equations,in press, I. Herrera, Instituto de Geofisica, National Autonomous Univer-
1989a. sity of Mexico, Apdo. Postal 22-582, 14000 Mexico D.F., Mexico.
Celia..M. A., I. Herrera, and E. T. Bouloutas,Adjoint Petrov- J. S. Kindred,GolderAssociatesInc., 4104148thAvenueNorth-
Galerkinmethodsfor multi-dimensionalflow problems,in Proc. east,Redmond,WA 98052.
SeventhInt. Coni on Finite ElementMethodsin Flow Problems,
pp. 953-958,UAH Press,Huntsville, Ala., 1989b.
Collatz, L., Numerical Treatment ofDifferential Equations, 3rd ed.,
Springer, New York, 1960. (Received April 21, 1988;
Herrera, I., Boundary Methods: An Algebraic Theory,Pitman, revised January 16, 1989;
London, 1984. accepted January 31, 1989.)