EULERIANLAGRANGIAN LOCALIZED ADJOINT
METHOD AND SMOOTHED AGGREGATIONS
ALGEBRAIC MULTIGRID
by
Caroline Heberton
B.A., Grinnell College, 1978
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
2000
This thesis for the Doctor of Philosophy
degree by
Caroline Heberton
has been approved
by
Andrew Knyazev
Haboo
Date
Heberton, Caroline (Ph.D., Applied Mathematics)
EulerianLagrangian Localized Adjoint Method and Smoothed Aggrega
tions Algebraic Multigrid
Thesis directed by Professor Thomas Russell
ABSTRACT
A threedimensional implementation of an EulerianLagrangian Lo
calized Adjoint Method (ELLAM) for the solution of an advectiondiffusion
equation is described. A system of integral equations is derived from a second
order partial differential equation, and the numerical treatment of each term
in the integral equation is discussed. Computational results for a variety of
test problems are presented. Analysis of a onedimensional problem is devel
oped. Stability is demonstrated in the purely diffusive and purely advective
limits, under assumptions on boundary conditions, and in the case of an infi
nite spatial domain with no diffusion. The effect of repeated application of the
operator on a concentration peak is investigated.
A theoretical result for the convergence of a variant of a smoothed ag
gregations algebraic multigrid (AMG) method for the solution of a discretized
secondorder scalar problem with jumps in coefficients is presented. An al
gorithm is developed with components shown to satisfy abstract convergence
criteria. The algorithm uses a coarsening strategy which respects boundaries of
m
highcoefficient subdomains, and introduces a tentative prolongator and pro
longator smoother designed to treat problems with discontinuous coefficients.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Thomas Russell
ACKNOWLEDGMENTS
Portions of Chapter One of this thesis are being published as a part of a U.S.
Geological Survey Water Resources Report by C. I. Heberton, T. F. Russell,
L. F. Konikow, and G. Z. Hornberger. Funding was also provided by the
USGS. Research was also supported in part by National Science Foundation
Grant No. DMS9706866 and Army Research Office Grant No. 37119GSAAS.
Chapter Two represents joint work with Petr Vanek and Marian Brezina.
CONTENTS
Figures ............................................................. ix
Tables.............................................................. xvi
Chapter
1. EulerianLagrangian Localized Adjoint Method ...................... 1
1.1 Introduction..................................................... 1
1.2 ThreeDimensional ELL AM....................................... 4
1.2.1 Governing Equation for Solute Transport........................ 5
1.3 Formulation of ELL AM equations................................ 8
1.3.1 Cell Integral Equations........................................ 9
1.3.2 Outflow Boundary Equations..................................... 9
1.3.3 Assumptions .................................................. 11
1.4 Numerical Methods............................................. 11
1.4.1 Mass Tracking ................................................ 12
1.4.2 Numerical Integration......................................... 14
1.4.3 Dispersion.................................................... 14
1.4.4 Storage at New Time Level..................................... 14
1.4.5 Mass Storage at Old Time Level ............................... 17
1.4.6 Approximate Test Functions.................................... 18
1.4.7 Source Integral............................................... 21
1.4.8 Sink Integral................................................. 22
vi
1.4.9 Inflow Boundary Integral ...................................... 22
1.4.10 Outflow Integrals............................................. 23
1.4.11 Decay......................................................... 25
1.4.12 Assumptions .................................................. 25
1.5 Test Problems................................................... 26
1.5.1 OneDimensional Flow........................................... 26
1.5.2 Uniform, ThreeDimensional Flow ............................... 29
1.5.3 TwoDimensional Radial Flow.................................... 33
1.5.4 Initial Condition in Uniform Flow.............................. 40
1.5.5 Constant Source in Nonuniform Flow............................. 42
1.6 OneDimensional ELL AM.......................................... 51
1.6.1 Limiting Case of No Advection.................................. 54
1.6.2 Limiting Case of No Diffusion.................................. 57
1.6.3 Infinite Spatial Domain ....................................... 58
1.6.4 Further Considerations for Finite Spatial Domain............... 63
1.6.5 Numerical Dispersion and Oscillations ......................... 65
1.7 Conclusion...................................................... 77
1.8 Further Research................................................ 77
2. Algebraic Multigrid for Problems with Discontinuous Coefficients . 78
2.1 Introduction.................................................... 78
2.2 Abstract Convergence Theory..................................... 82
2.3 Model Problem................................................... 86
2.4 Algorithm....................................................... 87
2.4.1 Stages of Processing........................................... 90
vii
2.5 Smoother properties................................................ 92
2.6 Tentative Prolongator ............................................. 99
2.6.1 Assumptions on aggregates........................................ 99
2.6.2 Kinds of aggregates............................................. 100
2.7 Single subdomain with highcoefficient............................ 102
2.8 Multiple subdomains with highcoefficients........................ 106
2.9 Computational Experiments......................................... 112
2.10 Conclusion....................................................... 114
References............................................................. 116
viii
FIGURES
Figure
1.1 Cell 7 receives advected mass from inflow, source in cell 1, and
storage in cell 1........................................... 10
1.2 Outflow boundary receives advected mass from inflow, source in
cell 1, and storage in all cells............................ 11
1.3 Preimage of cell may be irregularly shaped and not easily de
limited by backtracking. Instead, the known mass distribution
at time tn is tracked forward along streamlines of the advective
flow to time tn+l.............................................. 14
1.4 Approximate test functions in one direction on a uniform grid,
with NS=2...................................................... 19
1.5 Approximate test functions in one direction on a uniform grid,
with NS=4...................................................... 20
1.6 Approximate test functions in one direction on a nonuniform
grid, with NS=2................................................ 20
IX
1.7 Plots of concentration as a function of cell node for onedim
ensional flow with a constant velocity field and low dispersion.
Shown are the analytical solution (lowest graph), ELLAM re
sults using CELDIS = 1 (121 time steps), NSC = 32, NSR =
NSL = 2, NT = 128, and ELLAM results using CELDIS = 10.1
(12 time steps), NSC = 4, NSR = NSL = 2, NT = 128 (upper
graph). Results for CELDIS = 1 are virtually identical to the
analytical....................................................... 29
1.8 Plots of concentration as a function of cell node for onedim
ensional flow with a constant velocity field and high dispersion.
Shown are the analytical solution (lowest graph), ELLAM re
sults using CELDIS = 1 (121 time steps), NSC = 32, NSR =
NSL = 2, NT = 128, and ELLAM results using CELDIS = 10.1
(12 time steps), NSC = 4, NSR = NSL = 2, NT = 128 (upper
graph). Results using CELDIS = 1 are virtually identical to the
analytical....................................................... 30
1.9 Concentration vs. cell node plot with CELDIS = 61 (two time
steps), NSC = 4, NSR = NSL = 2, NT = 128......................... 32
1.10 Plots of concentration as a function of cell node for decay con
stant A = 0.01 sec1. Shown are the analytical solution (lower
graph) and ELLAM results using CELDIS = 1 (121 time steps),
NSC 32, NSR = NSL = 2, NT = 128......................... 32
x
1.11 Concentration contours of analytical solution in the horizontal
plane containing the solute source (layer 1) for threedimensional
solute transport in a uniform steady flow problem......... 35
1.12 Concentration contours in the horizontal plane containing the
solute source (layer 1) for threedimensional solute transport in
a uniform steady flow problem with CELDIS = 7 (two time
steps), NSC = NSR = NSL = 4, NT = 16............... 35
1.13 Concentration contours in the horizontal plane containing the
solute source (layer 1) for threedimensional solute transport in
a uniform steady flow problem with CELDIS = 1 (14 time steps),
NSC = NSR = NSL = 4, NT = 4................................. 36
1.14 Concentration contours in the horizontal plane containing the
solute source (layer 1) for threedimensional solute transport in
a uniform steady flow problem with CELDIS = 0.1 (134 time
steps), NSC = NSL = 4, NSR = 8, NT = 16............ 36
1.15 Contour plot of analytical solution for two dimensional radial
flow problem................................................ 37
1.16 Contour plot of concentration for run with two time steps using
CELDIS = 75, NSC = NSR = 4, NSL = 2, NT = 16................ 38
1.17 Contour plot of concentration for run with 29 time steps using
CELDIS = 5, NSC = NSR = 4, NSL = 2, NT = 4.................. 40
1.18 Contour plot of concentration for run with 563 time steps using
CELDIS = 0.25, NSC = NSR = 4, NSL = 2, NT = 4. Concen
tration maximum is 1.019.................................. 40
xi
1.19 Contour plot of concentration for run with 563 time steps using
CELDIS = 5, NSC = NSR = 8, NSL = 2, NT = 4. Concentra
tion maximum is 1.0056...................................... 41
1.20 Contour plot showing log of concentrations in analytical soluion
of Dirac problem at i = 90. Concentration maximum is 25195. . 44
1.21 Contour plot showing log of concentrations for spike initial con
dition, and CELDIS = 5, NSC = NSR = NSL = 4, NT = 2.
Concentration maximum is 15160........................... 44
1.22 Contour plot showing log of concentrations in analytical soluion
of Dirac problem at t = 130. Concentration maximum is 14539. 45
1.23 Contour plot showing log of concentrations for dispersed initial
condition, and CELDIS = 5, NSC = NSR = NSL = 4, NT = 2.
Concentration maximum is 16909................................ 45
1.24 Contour plot showing log of concentrations of analytical solution
to Dirac problem at t = 130. Concentration maximum is 8645. . 46
1.25 Contour plot showing log of concentrations for dispersed initial
condition, and CELDIS = 5, NSC = NSR = NSL = 4, NT = 2
at t 130. Concentration maximum is 8167..................... 46
1.26 Contour plot showing concentrations of a twodimensional finite
element solution of the Burnett and Frind problem. Contours
shown are 0.1 to 0.9........................................ 47
xii
1.27 Contour plot showing concentrations of ELL AM solution to two
dimensional Burnett and Frind problem using CELDIS = 30,
NSC = NSR = NSL = 4, NT = 32. Contours shown are 0.1 to
0.9............................................................ 49
1.28 Contour plot showing concentrations of a threedimensional fi
nite element solution of the Burnett and Frind problem. Con
tours shown are 0.1 to 0.9................................... 50
1.29 Contour plot showing concentrations of ELLAM solution to low
dispersion Burnett and Frind problem using CELDIS = 30, NSC
= NSR = NSL = 4, NT = 32. Contours shown are 0.1 to 0.9. . 51
1.30 Contour plot showing concentrations of a threedimensional fi
nite element solution of the Burnett and Frind problem with
high vertical transverse dispersivity. Contours shown are 0.1 to
0.9............................................................ 51
1.31 Contour plot showing concentrations of ELLAM solution to high
dispersion Burnett and Frind problem using CELDIS = 21, NSC
= NSR = NSL = 4, NT = 32. Contours shown are 0.1 to 0.9. . 52
1.32 Mass at old time level is integrated exactly using cell centers and
preimages of cell boundaries as integration points, and applying
the trapezoidal integration rule............................... 54
1.33 Initial and advected peaks with Cr = 10......................... 68
1.34 Initial and advected peaks with Cr = \.......................... 69
1.35 Initial and advected peaks with Cr = 1.......................... 69
1.36 Initial and advected peaks with Cr = ~.......................... 70
xiii
1.37 Initial and advected peaks with Cr = 10^................... 70
1.38 Eigenvectors of 10 x 10 S(A~1Bwith Cr = .................. 72
1.39 Initial condition and results using 10 x 10 (S(A~1B)^:)T with
T = 1,10,100,1000 and Cr = \............................... 73
1.40 Initial condition and results using 40 x 40 {S{A~lB)^)T with
T = 1,10,100,1000,10000 and Cr = \......................... 75
1.41 Initial condition and results using 200 x 200 (S(A~lB)^)T with
T = 1,10,100,1000,10000 and Cr = ......................... 75
1.42 Initial condition with 2 nodes on a front, and results using 40 x 40
(S(A15)
1.43 Initial condition with 3 nodes on a front, and results using 40x40
(S'(yl1B)^):r with T = 1,100 and Cr = \................... 76
1.44 Initial condition with 4 nodes on a front, and results using 40 x 40
{S{A~1B)^)T with T = 1,100 and Cr = ...................... 77
1.45 Initial condition with 5 nodes on a front, and results using 40x40
(,S{A~lB)&)T with T = 1,100 and Cr = \..................... 77
1.46 Superposition of results showing initial condition and advected
peak using Cr \. Grids are 5x5, and 4 refinements, each
by a factor of three. Initial condition is a spike on 5 x 5 grid.
Advected peaks show increasing height with decreasing Ax. . 79
1.47 Superposition of results showing initial condition and advected
peak using Cr = Grids are 5x5, and 4 refinements, each
by a factor of three. Initial condition is a spike on 5 x 5 grid.
Advected peaks show increasing height with decreasing Ax. . 79
xiv
2.1 Configuration with two highcoefficient subregions touching at
a single node........................................................ 118
xv
TABLES
Table
1.1 Parameters used in ELLAM simulation of solute transport in a
onedimensional, steadystate flow system.................... 28
1.2 Parameters used in ELLAM simulation of transport from a con
tinuous point source in a threedimensional, uniform, steady
state flow system............................................ 33
1.3 Parameters used in ELLAM simulation of....................... 37
1.4 Parameters used in ELLAM simulation of threedimensional
transport from a point source with flow in the xdirection and
flow at 45
degrees to x and yaxes................................... . 43
1.5 Parameters used for ELLAM simulation of transport in a vertical
plane from a continuous point source in a nonuniform, steady
state, twodimensional flow system (described by Burnett and
Frind, 1987).................................................... 48
1.6 Sections represent Courant numbers, Cr = 1/2, 1/3, 1/64,
1/500, 1/2000, 99/200, respectively............................. 66
1.7 Eigenvalue and coefficient of the respective eigenvector in eigen
vector decomposition of spike initial condition at node seven on
10 x 10 grid.................................................... 72
xvi
2.1 Checkerboard pattern, coefficients lO0, 10_tr................... 116
2.2 Two cubes touching at a node. Coefficient=10<7 in dark subre
gions (see Figure 2.1)............................................. 117
2.3 Element coefficients random in (10_
xvii
1.
ELLAM
1.1 Introduction
Solute transport in flowing groundwater has been intensively studied
in recent years in an effort to predict longterm effects of contamination, and
to design strategies for remediation. Convective and diffusive processes both
move contaminant through saturated subsurface media. Therefore, advection
diffusion transport equations, important in many areas of applied science and
engineering, are of particular concern to groundwater hydrologists active in the
assessment of aquifer contamination.
Solution of an advectiondiffusion transport equation may pose diffi
cult challenges to a numerical method. While diffusion equations tend to be
tractable using standard finite difference or finite element methods with a level
of discretization costeffective for application to large systems over extensive
time periods, admitting an advective component to the equation may require
prohibitively fine discretizations. Consider a onedimensional prototypical so
lute transport equation,
dc dc d2c
dt J*~Vdx dxx
= 0,
(1.1)
where c is concentration, v is velocity, and D is diffusion. Define the dimen
sionless Peclet number by
vL
D
where L is a characteristic length of the system, and the grid Peclet number
by
Pen
vAx
'9 D
where Ax is the mesh size of the grid. For mesh density such that Ax satisfies
Peg < 2
1
standard Eulerian numerical methods, such as centered finite differences or
Galerkin finite elements, perform well. These methods produce nonphysical os
cillations in concentration with larger grid Peclet number. Such oscillations are
often avoided by using an upstream finite difference or finite element method,
which effectively increases D enough to make Peg < 2. The resulting concen
tration profiles, however, are damped relative to their physical counterparts,
and also involve a grid orientation effect due to the introduction of a large
numerical dispersion term which is not rotationally invariant [83]. A solution
with neither spurious oscillations nor numerical dispersion indeed requires that
Peg < 2, meaning
n
X
A A
Ax ~ 2
. 2Z
(1.2)
(1.3)
where nx is the number of grid cells in the length L. Enforcement of (1.2) in
three dimensions is often impractical, since Peclet numbers in the hundreds are
common in field problems [].
The need for the scale (1.3) is actually an artifact of the numerical
methods and not the result of a physical process. A moving front solution
of (1.1) with a step function as the initial concentration profile, has the form
c(x, t) = erf
x vt
Vwt
J
where erf is the error function. The width of this front is proportional to VD,
hence to ^=. Eulerian methods need a number of cells proportional to \/Pe
on a front to satisfy Ax proportional to Thus the number of cells on a
moving front must increase as a problem becomes more advection dominated.
A method requiring Ax proportional to would require a constant number
of cells on a front, regardless of its width. This would be practical in three
dimensions, even with large Pe.
Eulerian methods also tend to suffer from large time truncation errors
when applied to advection dominated problems. These errors will be propor
tional to a power of At multiplied by a higher time derivative of the analytic
2
solution c, depending on the particular time stepping procedure. Errors will be
large when a steep front passes, unless time steps are very small. Time step
ping in a way which follows the flow, would permit larger time steps without
degradation of accuracy.
Subsurface contaminant transport problems involve steep, but not
discontinuous fronts because of the combined effects of advection and disper
sion. Problems may be advection or dispersiondominated, often varying in
different parts of the spacetime computational domain. Methods are sought
which are effective with these combined effects.
Eulerian and characteristic methods are both being advanced to ac
curately and efficiently solve advectiondiffusion transport equations. Eulerian
methods are characterized by the use of a fixed spatial grid. PetrovGalerkin
techniques [77, 16, 98, 22], optimal test functions [15, 25, 28, 10, 8], total varia
tion diminishing scheme [29], stabilized finite elements [46], and the streamline
diffusion finite element method [39, 63, 20, 64, 62, 51, 67, 66, 69, 68, 74, 102,
103], are alternatives to classical difference or element methods on fixed grids.
Each incorporates some strategy to minimize or to tune space or time trunca
tion error in order to avoid oscillations, while introducing as little numerical
dispersion as possible. They tend to need restricted size of the time steps for
good accuracy.
Characteristic methods treat separately the advective and diffusive
components of the transport equation, tracking along characteristics of the
flow, then solving a diffusion equation on a fixed grid. These methods in
clude the method of characteristics [48, 78, 11, 59], modified method of char
acteristics [65, 42], characteristic Galerkin method [32, 93], transportdiffusion
method [79], characteristicmixed finiteelement method [2, 101, 3], operator
splitting method of [40, 100, 30], the LagrangianGalerkin method [76], and
EulerianLagrangian methods as discussed below. These methods offer reduced
time truncation error compared to Eulerian methods, and can use larger time
steps. Primary drawbacks are difficulties in mass conservation and in rigorous
formulation of boundary flux.
3
The EulerianLagrangian localized adjoint method (ELLAM) main
tains mass conservation and provides a framework for treatment of general
boundary conditions. ELLAM was introduced by Celia et al. [26], Russell [81],
and Herrera et al. [58] for constantcoefficient onedimensional equations, then
extended to onedimensional variablecoefficient cases by Russell and Tru
jillo [82], Wang [94], and Wang, Ewing, and Russell [96]. Celia and Fer
rand [24] and Healy and Russell [55] extended ELLAM to a onedimensional
finitevolume setting. Onedimensional nonlinear advectiondiffusion equations
have been treated by Ewing [33] and Dahle, Ewing, and Russell [31]. Addi
tional work has been done on ELLAM schemes for onedimensional transport
equations involving reaction [34, 43, 35, 27, 36, 37, 85, 41, 95, 38].
Implementation of ELLAM schemes in multiple spatial dimensions
involves additional concerns relative to the onedimensional case [82], Various
twodimensional ELLAM schemes have been developed [94, 43, 35, 97, 12,
13, 56], including one with optimal order error estimates proved by Wang.
Celia [23] explores a threedimensional ELLAM in a framework within which
all characteristic methods can be viewed.
In this thesis, a threedimensional ELLAM implementation is de
scribed, and numerical results are provided. A forward tracking approach is
taken to advection, and dispersion is treated by an implicit formulation in time.
The methods apparent robustness is discussed in a onedimensional context.
1.2 ThreeDimensional ELLAM
A threedimensional finitevolume EulerianLagrangian localized ad
joint method (FVELLAM, [55]) code has been developed to simulate three
dimensional solute transport in flowing ground water for a single dissolved
chemical constituent and represents the processes of advective transport, hy
drodynamic dispersion (including both mechanical dispersion and diffusion),
mixing (or dilution) from fluid sources, and simple chemical reactions (in
cluding linear sorption and decay). This code has been implemented as an
alternative algorithm within the U.S. Geological Survey (MOC3D) transport
model. It is integrated into the MOC3D model, which is itself integrated with
4
MODFLOW96, the USGS threedimensional, finitedifference, groundwater
flow model [52], [54], [53]. The code is written in FORTRAN77.
ELLAM [26] solves an integral form of the solutetransport equation.
The code uses an implicit method for dispersion calculations, which allows
for large time steps without stability constraints. The EulerianLagrangian
approach involves tracking mass through time, then solving a dispersion equa
tion on a fixedinspace grid. This is particularly advantageous for advection
dominated problems, as are typical of many field problems involving ground
water contamination, since EulerianLagrangian approaches tend to generate
less numerical dispersion than standard Eulerian finitedifference or finite
element methods.
1.2.1 Governing Equation for Solute Transport
The solute transport equation is the same as presented in Konikow
and others ([71], equation 6):
cfc (eD>Tj^) Zc'w (1.4)
+A (eC f QbC)
= 0
(summation over repeated indices is understood), where C is volumetric con
centration (mass of solute per unit volume of fluid, ML~3), g& is the bulk
density of the aquifer material (mass of unit solids per unit volume of aquifer,
ML~3), C is the mass concentration of solute sorbed on or contained within
the solid aquifer material (mass of solute per mass of aquifer material, MM1),
Â£ is the effective porosity (dimensionless), V is a vector of interstitial fluid ve
locity components (LT1), D is a secondrank tensor of dispersion coefficients
(L2T1), W is a volumetric fluid sink (W < 0) or fluid source (W > 0) rate
per unit volume of aquifer (T1), C' is the volumetric concentration in the
sink/source fluid (ML3), A is the reactive decay rate (T1), t is time (T), and
Xi are the Cartesian coordinates (L). Porosity is the ratio of pore volume to
5
the bulk volume of the porous medium. Interstitial velocity is the velocity that
fluid moves through the pores of the subsurface medium to achieve a given flow
rate across the faces of a volume of the saturated medium.
The terms controlling sorption are combined into a single parameter,
the retardation factor (Rf), assumed to be constant in time, since we consider
a linear phaseequilibrium relationship in which C is proportional to C. The
retardation factor is defined as:
Rf = 1 +
6bC
eC '
An integral form of the solutetransport equation, which is a state
ment of conservation of mass over the domain of integration, is the governing
equation for this finitevolume ELLAM approach. Equation (1.4) is integrated
against a spacetime test function to provide the formulation, including treat
ment of (cell or subdomain) boundary conditions and solute decay.
The test function effectively specifies the domain of integration for the
transport equation by the portion of the spacetime domain where its value is
nonzero. On a subdomain of integration, the test function can be seen as an
integration weight at each point. Varying the weight along streamlines of the
flow is a convenient mechanism to provide solute (growth or) decay.
Divide (1.4) by Rf, multiply by a test function u(x,t) and integrate
over time and space. Assuming Rf is constant in time, we have:
l r (^+
(1.5)
where is the entire spatial transport subdomain, and T is the end of the sim
ulation time period starting at time zero. Integrate equation (1.5) is integrated
by parts using
d(eC) d(ueC) du ^
u dt = ~Jk msC
and
or 1 1
V (eCV sT>VC) = V (ueCV usDVC) Vu (eCV sBVC)
lif Hf Rf
6
to yield the global equation,
f fT (9(nsC)
Jn Jo \ dt
+ ^V (ueCV ueDVC) + JVu eDVC
Rf Rf
u
Rf
yew
,du V
Â£C{m + R/Vu~Xu)
dtdx = 0.
(1.6)
The EulerianLagrangian aspects of the method derive from the re
quirement that the test function satisfy the adjoint equation, ^ Vw
Xu = 0 Thus, for the time step from tn to tn+1, we choose u of the form
u(x,i) = /(x, t)e_A(tn+1_t), where ^ + jq V/ = 0, so that / is constant
along characteristics of the retarded interstitial velocity field. Note that with
u = eHtn+1t) (that is, / = 1) in the following, we arrive at a statement of
global conservation of mass over the time step:
Jn(ueC)n+1dx fa(ueC)ndx
+ fa /<"+ 57 V (CV eDVC) + Â£7 Vu eDVC
g'ECWdtdx.
Kf ^
= 0.
To obtain a system of equations, each representing mass conservation
on a cell, we use u = J2iUi where l is the index for the finitedifference cells
Oi in the transport subdomain, and local spacetime test functions are defined
with /;(x, tn+1) = 1 on Oi and //(x, tn+1) = 0 elsewhere:
e\(in+1t) on characteristics from any (x,t) Gflx (tn,tn+1)
< into Oi at time level n + 1 (1.7)
0 otherwise.
7
We thus arrive at a system of local ELLAM equations,
/n,(Â£Cr+1dxeJA1/.(<) Ac
+ S (,CV eDVC) n dtds
suppuiC\Tn+1 ^
 C' /
IJ *w'*;T.C'W dtdx
supp uiHsupp W 1
= 0
where d signifies the spatial boundary of the argument; supp denotes the
support of a function, that is, the closure of the part of its domain where a
function assumes a nonzero value; n is the unit outward normal vector on the
specified boundary; is the preimage at time tn under advection of cell fl; at
time tn+1; Tn+1 = dQx (tn, tn+1) is the spacetime boundary at time step n +1;
t denotes time, and dx and ds signify differential volume and area, respectively.
Note that equations (1.8) appear as spacetime integrals of diffusion
equations. ELLAM can be viewed as a method of characteristics, tracking
mass along streamlines of the flow to accumulate data to the righthand side
of the system of equations.
1.3 Formulation of ELLAM equations
ELLAM approaches the hyperbolicparabolic nature of the solute
transport equation by combining a method of characteristics technique for
advection, with a backward Euler in time and centered differences in space
solution to a diffusion equation. The details of these approximations depend
on the evaluations of the integrals in (1.8).
ELLAM requires knowledge of fluid velocity everywhere within the
domain of the transport problem. Consequently, a flow equation is solved on
a potentially larger domain, and output is used to solve the solute transport
equation on a subdomain. The flow equation is an elliptic or parabolic pressure
(groundwater flow) equation, which is solved using the USGS finitedifference
code, MODFLOW.
8
This ELL AM method approximates total solute flux across the trans
port subdomain boundary by advective flux. It should be noted that this
approximation is not demanded by ELLAM methods, in general, but is a fea
ture of this particular implementation. This approximation means that bound
ary face concentrations are not coupled to cell center concentrations through
the numerical derivative (concentration gradient). All mass moving into and
out of the transport subdomain can be tracked using the advective algorithm.
Mass tracked across outflow boundaries provides data for a system of outflow
boundary equations decoupled from the cell equations. Userinput inflow con
centrations, together with the outflow solutions, then appear on the righthand
side of the system of cell equations representing local statements of mass con
servation. These cell integral equations are solved for Cn+1, concentration at
the new n + 1 time level at each cell center.
1.3.1 Cell Integral Equations
Take the conservation of mass equations for each cell (1.8), approxi
mate the total boundary flux with advective flux, the dispersion time integral
with a backward Euler formulation, and then rearrange terms. The system of
equations to be solved is then,
So,(eC)"+1 dx At Jdai A.(eUVC)"+1 n da = e fn; (eC)" dx
It e n dtds ,,
suppulnvn+l 3 \ m J
+ If eA(tn+1t) Y'C'jf dtdx.,
supp u; flsupp w 3
where 12* means the preimage in the spatial domain at tn of 12/ at tn+1, Tn+1 =
<912 x (tn,tn+1) is the spacetime boundary at time step n + 1, n is the unit
outward normal vector, and supp f = {x\f{x) ^0}. Note that the right
hand side of (1.9) consists of advective mass contributions from storage, inflow
boundary, and sources. (The term storage will refer to solute mass already
in the transport system at the end of the last time step.) Figure 1.1 illustrates
the possibility of all of these mass contributions being nonzero.
9
Figure 1.1. Cell 7 receives advected mass from inflow, source in cell 1, and
storage in cell 1.
1.3.2 Outflow Boundary Equations
The term in (1.7) expressing mass crossing the transport subdomain
boundary during a time step is:
In ftn+1 V (ueCV ueDVC) dtdx =
C1 Jda ~{ueCV ueDVC) n dsdt.
Considering just the outflow portion of the boundary, this becomes,
jT f(dau,, ueDVC) n dsdt
kd^)outflow (UÂ£Coutflowfij) n ^S>
(1.10)
where total flux across the boundary is now approximated by advective flux.
We index the outflow boundary faces with ll and define the following
test functions:
{ea(tn+1t) on characteristics from O at time level n
into boundary area (dCi)u at any time during time step
0 otherwise.
(1.11)
The mass across outflow boundary face ll is the mass stored at the previous
time level which flows across the face, together with any inflow and source
mass that both enters the transport subdomain, and leaves through cell face ll
during the time step, as illustrated in figure 1.2, is given by,
^dtds = (1.12)
10
Figure 1.2. Outflow boundary receives advected mass from inflow, source in
cell 1, and storage in all cells.
e_AAt/(ao)*(^)nrfx
If e^'^eCinfio^ndtds
supp uundnin fiow x (tn ,tn+1)
+ If gA(tn+1t) Y^C'jjr dtdx.,
supp un C\supp Wx (tn ,tn+1)
where (dVt)u is a boundary face, and (dÂ£l)*u is the preimage at time tn of
(dQ,)u x [tn,tn+l]. This system of equations will be solved for Coutfi0w
1.3.3 Assumptions
As described by Konikow and others [71], a number of assumptions
have been made in the development of the governing equations for a coupled
system modeling flow and transport. For the derivation of the transport equa
tion, we assume:
Chemical reactions do not affect the fluid or aquifer properties.
The dispersivity coefficients are constant over a flow time step, and the
aquifer isisotropic with respect to longitudinal dispersivity.
1.4 Numerical Methods
The groundwater velocity field needed for solution of the transport
equation is obtained using the USGS MODFLOW code. This implicit finite
difference code yields a head (pressure) distribution which is calculated for a
given time step or steadystate flow condition. The specific discharge, or flow
11
rate per unit area, across each face of every finitedifference cell within the
transport subgrid is calculated and used to find the interstitial velocity at any
point using an interpolated value of the specific discharge, divided by porosity.
1.4.1 Mass Tracking
For each cell in the fixed finitedifference grid, the integrals on the
righthand side of equation (1.9) represent solute mass advected into the cell
during the time step from storage, the transport subdomain boundary, or a
fluid source, respectively. These integrals are formed by accumulating mass
tracked forward along characteristic curves of the velocity field, determined
from the MODFLOW solution as described above. Velocity at any point is
determined by linearly interpolating between cell faces in the direction of the
component of interest, and piecewiseconstant interpolation in the other two
directions.
. A system of three ordinary differential equations is solved to find the
characteristic curves [x = x(t), y = y(t); z = z{t)\ along which the fluid is
advected:
dx Vx
dt Rf
dy=Vy_
dt Rf
(1.13)
(1.14)
dz K
dt Rf
(1.15)
This is accomplished by introducing a set of moving points that can be
traced within the stationary coordinates of a finitedifference grid. Each point
corresponds to one characteristic curve, and values of x, y, and z are obtained
as functions of t for each characteristic [47]. Each point moves through the
flow field by the flow velocity acting along its trajectory.
The ELLAM equations, (1.9) and (1.12) suggest that mass is tracked
backwards along characteristics to the preimage of each cell or boundary face.
It is not possible to exactly locate all of the mass at the previous time level
by backtracking a finite number of points, however. In order to achieve mass
balance, the known mass distribution at the old time level is tracked forward
12
to the new time level. (See figure 1.3.) This approach ensures that mass will
not be lost or gained in advection from one time level to the next, but doesnt
guarantee an accurate mass distribution at the new time.
ELLAM tracks points which are the centers of volumes of fluid. Thus
mass in a fluid volume is tracked under advection during a time step, dis
tributed among destination cells, and accumulated to the righthand side stor
age, inflow, or source integral for each cell.
1.4.2 Numerical Integration
We will now discuss the numerical treatment of (1.9) and (1.12). The
j,i,k subscripts for a cell will denote the spatial finitedifference grid in
dexing, with this index order indicating an x,y,z sequencing. A righthanded
coordinate system, with the vertical index increasing from top to bottom is
used.
The equations are first divided through by porosity, which is repre
sented by a piecewise constant function in space and time. Each individual
term will now be discussed.
1.4.3 Dispersion
Time integration is accomplished using a one point in time backward
Euler rule. Spatially, a one point integration rule with a seven point stencil is
used:
where m = 1,2, 3 is the summation index for the dispersion term. Finite
difference approximations to the space derivatives in the dispersion integral
(1.16)
13
Figure 1.3. Preimage of cell may be irregularly shaped and not easily de
limited by backtracking. Instead, the known mass distribution at time tn is
tracked forward along streamlines of the advective flow to time tn+l.
14
are calculated using centered differences as in M0C3D, generalized for varying
grid dimensions. (See [71], page 64, for the form of the ebDmS^ expansions.)
UXm
1.4.4 Storage at New Time Level
The quantity mass/porosity in a cell at the new time level tn+1 is
expressed using the trapezoidal rule for integration, formulated over each cell
octant. Concentrations at octant corners are weighted averages of neighboring
node concentrations, determined by trilinear interpolation.
For each octant,
mass _ X Ar Anhn+l Y'8 C*
porosity ~ ^corner=l ^ corner ('ll?!
^2corner=l ^2nbr=l(WeW^)nbrC'nbr,
where for an interior octant, nbr is one of the eight neighboring grid nodes
between which concentration varies trilinearly. In case of a boundary octant, a
boundary face value is needed for calculation, and is taken to be the following:
Inflow user input;
No flow same as associated interior node;
Outflow calculated using cell parameters, boundary flow rate, and
mass tracked across boundary during transport time step.
Coefficients calculated by (1/8) octant volume nodal weight for all nodes
neighboring a cell comprise the storage matrix entries for the equation for each
cell. Boundary terms are put on the righthand side of the equation because all
boundary face concentrations are determined before the solution of the interior
equations.
It should be noted that linear interpolation is approximate in the
case where adjacent cells in the same layer of the the transport subdomain
have varying thicknesses, as is allowed by MODFLOW. Extreme variations
could adversely affect accuracy of the solution.
For an interior cell with all neighbors active and of the same thickness,
using b values at time n+ 1:
((
Axj
A;Cj_i + Ax.
) + (
A xj
Axj+i + Ax.
)4)
15
1 Ayii + A yi Ayi+1 + A y{
x ((
b>i_i
j,i,k
) + (j
bi
j,i,k
bj,i,k + bjt itki bjtitk + bjjijk+1
 ((*^rr)(T^r)(;
) 4)CJM
bj,i,k
A.Xj1  AXj Aj/j_x I Aj/j bj^fe I bj^^i
__ ^ Axj ^ ^ Ayi j ^ &j,i,k
)Cji,ii,kl
Axjx f A Xj A?/j_i + Ayi bj+ bjjy
i,i,k+1
)Cii ,i
i 1,6+1
+ L. 6"'*
'Ax^i + A xj Ayi+i + Ay{ bjtitk + b
j,i,k+1
(M+1,6+1
+ (
+ (
Ax,
)(
A yi
j yi bjtitk
Axj_i + A Xj A y,+i + Ayi bjt frj,i,6i
Ayi / Ax, s b
)C:
)(
)(]
j,i,k
Ayii Ayi Ax^+i A Xj bjj^ bjj^
jl,i+l,kl
Ax,
)(l
b
'j,i,k
1,6+1
Ay,_i I Ayi Axj+i ) Axj bjj^ + ^++6+1
i / A^j w Ayi bjtitk ^
+ , A_ Ar^rrMin ;w+m+i,6i
Axj+i + Axj Ayi+i + Ayi fr+i,6 ~f~ ^j,i,6i
+ (
Ax,
)(
Aj/i
)(]
j,i,k
Axj+i  A Xj A yi\\ ~\~ Ayi bjj^ ^,,,6+1
)&j+l,i+l,k+l)
+ ((
x ((
A Xj , . A Xj .
J ) + (xir^r)4)
Axj_i + A Xj Axj+i + A Xj
A Vi bjjtk
+ (
Ay,_i I Ayi bjj^ I bj^^
AVi ^ ^ bjjfc
)Cj,i1,61
Ay,i H Ayi bjjtk + &j,i,6+i
;,i1,6+1
, ( A yi ^ bj^k xsi
+ 1TTTTAi
Ayi+i I Ayi bjjfi f" bjtitki
+ (
A yi
)(
bi_i
j,i,k
Ayi+i Ayi bjjfe + bji k+i
)Cj,i+l,k+l)
X ((
Ay,_i + Ay/ v Ayi+l + Ayi
Ax, x, bjjtk
)(i
Axji + Axj bjjfi + bjjt&+i
)Qi,
,i,A:+l
, A Xj . bjj^k
Axj+1 + Axj bj^k + bjjfii J+1hk~1
16
+
42
+
o'
I
Si
I
O'
+
Si
+
ci
H
+
<1
+
H
<
+
^5
+
I
r
rH
+
+
o
+
I
+
fn
rH
I
o'
I
+
Si
H
<
+
I
Si
O
+
Si
535
<]
H
<
+ s
+ <
<3 .
5=3
+ S'
7 <
!S>
535
55 55
<1 <1 <1
+ 5$ + 55 +
+ <1 7 < rH 4"
535 55 sS
<1 <1 <3
<3
S
+ ^
t, <
i
>
< 
535
<
+
H
535
<
+
<
1 *i
+ <
+
535
<]
+
535
<
H
<
H
<
H
<1
H
<1
j7 <\ ^<3 ^<3
+ + + x + +
<1
+
rH
I
s
<
rH
1 *$Â£
<1 <
se *T rcT 55 + 55 + s
rC^ + <1 rH 1 <3 rH <
4e 1 >S> 1
55
rcT < <1
I
I
I
Si
o"
vf<
H
H
<1
+ ^
t <1
>
< S
+ <1
ii
+
< i
535
<1
+
rH
+
s
<1
H
+
535
<1
+
rH
t
535
<1
+
>
55
+ <1
a
Sft
<1
+
535
<1
5>5
<1
I *
?o
r1
H
5>5
<3
+ s
535
<1
+
rH
4"
P5
<1
535
<3
4
4e
fO
+
Si
+
+
H
< s
I
+ <1
rH
+
H
<1
535
o
+ H.
rH "'T
535
<1
I
42
+
X
+
+

X + + 1
t"
X(((
+ (
X ((
Axj
Axji + A Xj
A yi
A yi+i + A yi
A Xj
Axji + Ax.
) + (
A Xj
)
A Xj+i + A Xj
A yi
)4)((
A Vix + Ay,
) + (
A yi
Ayii + A yi
A yi
)Cjl,il,k + (
A Xj
Axj+i + Ax3
A yi+x + Ay,
)4)
1.4.5 Mass Storage at Old Time Level
The total mass advected into each cell during a transport time step
that was already stored within the system at the old time level is needed for
the righthand side of the ELLAM equation. Numerically, this is accomplished
by tracking mass forward from the old time level, tn, along characteristics.
Each cell is divided into subcells determined by parameters NSC, NSR, and
NSL, specifying the number of subcells in the column, row, and layer direction,
respectively. The center of each subcell is tracked through the time step under
advection. Depending on the exact location of this point in the destination cell
at the new time, all of the mass in the subcell may or may not also be found
in that destination cell. In order to mitigate the effects of unwarranted mass
lumping, subcell mass is distributed among cells neighboring the destination
cell using the approximate test functions, wi, described below. The value of
wi at the subcell center destination point is the fraction of subcell mass to be
distributed to cell Â£li.
This yields the formulation,
,aa t
In;r &
dx. = e \ \ _
h (NSC){NSR)(NSL)
subcell
center
Wi
(p')C(p)
(1.18)
where summation runs through all subcells of each cell in the transport sub
domain, and pf is the image of p under forward tracking to the new time level.
1.4.6 Approximate Test Functions
For each active cell, an approximate test function is defined for the
purpose of distributing advected mass from the old time level among neighbor
ing cells at the new time level. The designation approximate test function
18
1
Figure 1.4. Approximate test functions in one direction on a uniform grid,
with NS=2.
is given since the graph of this function looks like a characteristic (indicator)
function with slant sides extending into adjacent cells, while the test functions
described in the derivation of the governing equation are exactly character
istic functions at time tn+1. An approximate test function is determined by
NSC, NSR, and NSL, the proximity of the transport subdomain boundary, and
the active status of neighboring cells. Mass is not split across the transport
boundary or into inactive cells.
An approximate test function associated with interior cell Oj is the
product of a onedimensional approximate test function in each direction:
On a uniform grid, we define reference coordinates x, y, z with respect
to a cell Q,i with node indices j, i, k by
and similarly for y and z. For an interior cell on a uniform grid with all
surrounding cells active,
and similarly for g and h. These functions are shown graphically in figures 1.4
and 1.5.
wjik(x,y,z) = f(x)g{y)h{z).
X Xi
0
NSCx + \(NSC + 1)
2 2 NSC
2 NSC
< X < ~ \ +
2 NSC
1
/(Â£)
1
2 ^ 2NSC
2NSC
NSCx + (NSC + 1)
0
1 _ i
2 2 NSC
1 + 1 <
2 ^ 2NSC
< X < 2 + 2Ni
< X
2 NSC
19
1
Figure 1.5. Approximate test functions in one direction on a uniform grid,
with NS=4.
On a nonuniform grid, test functions are scaled by grid ratios. In
one direction on a nonuniform grid, with neighboring cells active, approximate
defined generally, on a potentially nonuniform grid by,
and similarly for g and h. There are seven more test functions which may be
nonzero at any given point in cell j, i, k. If x, y, z > 0 their values given by,
wj+i,n = (! /0=)) g{y) h(z),
Wj,i+i,h = f(x) (1 g(y)) h{z),
Wji,k+i = f(x) g{y) (1 h(z)),
Wj+i,i+i,k = (1 f{x)) (1 g(y)) h(z),
Wj.i+i.fc+1 = f(x) (1 g(y)) (1 h(z)),
wj+i,i,k+i = (l /(Â£)) g{y) (1 h(z)),
wj+W)k+i = (1 f{x)) (1 g{y)) (1 h(z)).
Figure 1.6. Approximate test functions in one direction on a nonuniform grid,
test functions corresponding to / defined above are illustrated graphically in
figure 1.6.
For a point in Qi with reference coordinates x,y,z Â£ [1,], / is
1 i
2 > 2
/(Â£)=<
with NS=2.
20
For x, y, z in any other octant of the cell, there are also eight potentially nonzero
test functions evaluated similarly. Thus equations of the form (1.19) for /, g, h,
defined explicitly for reference coordinates in the interval [ , ], are sufficient
to evaluate all test functions at any point, and serve as a global definition.
The sum of the values of the test functions at any point in the trans
port domain is one, thus conserving mass for the integral equation.
The test function component in the direction normal to the boundary
extends from the center of a boundary cell to the boundary face with the value
of one. Thus there is no splitting of mass across the boundary.
There is no test function associated with an inactive cell. In a cell
adjacent to an inactive cell, the value of the test function that would normally
be assigned to the inactive cell is distributed proportionally to other test func
tions that are nonzero at that point. All test functions are zero in inactive
cells.
Extreme variation in cell thickness among neighboring cells in a layer
may adversely affect model results. In this case, linear interpolation of concen
tration or approximate distribution of advected mass may be inaccurate.
1.4.7 Source Integral
The last term in (1.9) pertains to sources and sinks within the trans
port subdomain. ELLAM assumes a source or sink acts uniformly over the
' entire cell surrounding a source or sink node.
To treat a source, a single time step is discretized into a number of
subtime steps determined by parameter NT, and the composite trapezoidal
integration rule is applied in time. This time discretization evens the spatial
distribution of incoming mass: Mass is tracked to varying locations within
the transport! subgrid depending on when in the transport time step the mass
enters the system. At each subtime step, inflow mass is spatially discretized,
tracked, and accumulated just like mass already in the system at the start of
the transport time step, but for the shorter interval.
Mass from sources is accumulated to the righthand side of the local
conservation of mass equation for cell fh with the following integration. The
21
domain of integration is all sources which intersect the spacetime test function
for cell Q,i, so that all source mass that flows into Ot during a time step is
integrated. Multiple sources within a source cell are summed.
If
supp uiDsupp W
~ 1
E
all
sources source
subcell
center
dtcht
NT
Y Y eWn+1~trn)
p= m=0
A T
NSC NSR NSL
Wi
(/) E c.%
= Kf
where summation runs through all sources in all subcells of every source cell
in the transport subdomain; pf is the image of p under forward tracking to the
new time level; Qs is the flow rate of source s in the source cell; AT =
or if m = 0 or NT; and tm = tn + represents the time during time
step at which discretized source mass enters the system.
1.4.8 Sink Integral
Analytically, the domain of integration is the support of the space
time test function for a cell Oi, intersected with any sink cells. To approximate,
this term is formulated only if Ot is a sink, and sink concentration is assumed
to be the average nodal concentration for the transport time step, with the ex
ception of a sink due to evapotranspiration, where sink concentration is taken
to be zero. Integration rules are a one point in space and a one point backward
Euler in time. The averaging of concentration results in this integral approx
imatation contributing to both the left and righthand sides of the equation
for sink cell Ot with coordinates j, i, k. This gives
22
ff eA(tn+1t)i dtdkx.
suppuiDsupp W f
= sink HC^e E ^ + C'et^) dx
( \
At (C"+1+Cn)
Â£j,i ,fc 2
(1.21)
where Q is sink flow rate; et refers to evapotranspirative flux; and multiple
sinks in a sink cell are summed.
1.4.9 Inflow Boundary Integral
For an inflow boundary integral as for a source, a single time step is
discretized into a number of subtime steps determined by parameter NT. The
composite trapezoidal rule is applied in time. At each subtime step, inflow
mass is spatially discretized, tracked, and accumulated just like mass already in
the system at the start of the transport time step, but for the shorter interval.
The only difference in the treatment of the inflow boundary from the treatment
of the source is that only the twodimensional boundary face is discretized,
while for a source, the entire cell is discretized. For a cell fi/, the integration
is performed over the intersection of the spacetime test function for that cell
and the transport subdomain boundary; that is, all mass entering through the
boundary and advected to Qi during the time step is accumulated to the right
hand side of local equation l. Mass advected into 0; during a time step from
the inflow portion of the transport subdomain boundary is give by,
JJ e
suppu(nrTi+1 ^
n dtds
NT
E E E
all P m~l
inflow face
faces subarea
center
gA(tn+1tm)
Arm
ref area
Wl {p^CinflowQ
(1.22)
23
where ref area = NSC NSR, NSC NSL,..., depending on plane of face;
pf is the image of p under forward tracking to the new time level; Q is inflow
rate; tm = tn + At represents the time during time step at which discretized
source mass enters the system; and AT = or if m = 1 or NT.
Summation runs over each p on the transport subdomain inflow boundary,
with the approximate test function wi used to select mass advected to Clt.
1.4.10 Outflow Integrals
Concentration is calculated at each outflow boundary face using cell
parameters, velocity information from MODFLOW, and the amount of mass
tracked across the cell boundary determined by ELLAM.
On the left hand side of the system of boundary equations (1.24) is an
integral approximated using a one point in space, one point backward Euler in
time formulation. This time approximation eliminates the exponential factor,
and yields,
I{da)u ft+1 e~X('tn+1~^Coutfiow^j n dtds
f(dQ)u ^'outfiowjij n ds (1.23)
~ AtQn Coutfi ow
where ll is the index for boundary faces; and Qu is the outflow rate across
boundary face ll, determined using the outflow velocity calculated from MOD
FLOW output and cell parameters. The concentration on face ll is the un
known in the boundary equation.
The righthand side boundary integrals are constructed from the mass
contributions tracked across the boundary from interior cells, sources, and
inflow boundaries during the transport time step. All mass associated with
a tracked point that reaches the outflow boundary at any time during the
time step is considered to leave the transport subdomain. Test functions are
evaluated to distribute mass among neighboring boundary outflow faces.
Numerically, the righthand side of (1.12) is accomplished using a
trapezoidal in time, midpoint on subcell in space approximation:
24
ff eKtn+1t)(j2L n dtds
*A/ Kf
suppunnrn+1
AA t
E E
p=
subcell
AxjAyib
(NSC) (NSR) (NSL)
wu(pf)C(p)
NT
+ EIE
all P= m= 0
inflow face
faces subarea
center
A(tn+1tTO)_
AT
ref area
IVll ) ('inflow Q
(1.24)
Â£j,i,k
NT
E EE'
p= m=0
A
AT
A'6'C' iV^L
^(pO E C'sf
sources source
subcell
center
cell
where the approximations to the integrals of stored, inflow, and source mass
are as in (1.18), (1.22), and (1.20), except that the approximate test function
for cell Cli has been replaced by wu, the approximate test function for boundary
face ll. Boundary face test functions are defined analogously to those associated
with cells, but only using factors in the two directions parallel to the boundary
face. No mass is distributed to a neighboring boundary face which is not part
of the outflow boundary.
We thus have a system of equations represented by a diagonal matrix,
to be solved for Coutfi0w
1.4.11 Decay
When simulating linear decay, all mass in the system at the beginning
of each transport time step is decayed by a factor e~Xt, where A is the decay
rate. Inflow and source mass are decayed in the same way, where the time
interval is now not the entire time step, but the portion of it during which new
mass is in the transport subdomain. This decay algorithm has no numerical
stability restrictions associated with it. If the halflife is on the order of or
smaller than the transport time step, however, some accuracy will be lost.
25
When a solute subject to decay enters the aquifer through a fluid
source, it is assumed that the fluid source contains the solute in the concentra
tion specified by C'. The ELLAM simulator allows decay to occur only within
the groundwater system, and not within the source reservoir. In other words,
for a given stress period, C' remains constant in time.
1.4.12 Assumptions
Implicit in the above numerical treatment of the terms of the integral
equations are the following assumptions:
Concentration at an outflow boundary face at the new time level is
well approximated by the mass crossing the face during the time step
divided by the fluid volume across the face.
Mass in or out of the transport subdomain during a time step via a
source or sink cellis well approximated by the average nodal concen
tration during the time step times the fluid volume through the source
or sink.
Cell thicknesses are smoothly varying within a vertical layer.
Transport subgrid boundaries are assumed to be far enough from the
plume that any errors in the treatment of the boundaries will not have
a significant effect on the solution. The boundary condition is that
the normal component of the concentration gradient on the boundary
is zero, so there is no dispersive flux across the transport subdomain
boundary.
1.5 Test Problems
The threedimensional ELLAM code was tested by running the same
suite of test cases as was applied to the USGS threedimensional method of
characteristics code, MOC3D. These tests were first used by Konikow and
others [71] to evaluate the MOC3D Version 1 implementation incorporating an
explicit formulation of the dispersion equation, then by Kipp and others [70]
for MOC3D Version 2, an implicit method. These benchmark problems are
used to evaluate features of a solute transport code with relevance to field
26
applications.
1.5.1 OneDimensional Flow
The first test case is a simple onedimensional system involving solute
transport in a finitelength aquifer with flow of constant velocity. Boundary
conditons are thirdtype, although ELLAM approximates total flux bound
ary conditions using advective flux. The numerical results are compared to a
solution by Wexler [99]. Parameters for the model are specified in table 1.5.1.
A low dispersion and a high dispersion case are presented here. In
both cases, ELLAM results for CELDIS = \ (241 time steps), NSC = 4, NSR
= NSL = 2, NT = 128 are essentially identical to the analytical results, and so
are not plotted. Instead, the results for substantially fewer time increments are
shown. For the low dispersion case, Dxx = 0.01, or oil = 0.1cm. The analytic
solution and two ELLAM solutions are graphed in figure 1.7. For CELDIS
= 1, which requires 121 time steps, there is a very close match between the
analytical and numerical solutions, using NSC = 32, NSR = NSL = 2, NT =
128. Setting NSC = 4 results in a slightly low concentration at the first grid
node, but quite accurate values elsewhere. The NSC = 4 solution is computed
in about a quarter of the time used for the finer, NSC = 32, discretization
on a Data General Unix workstation [57] (4 min 20 sec for NSC = 32, and 1
min 0 sec for NSC = 4). For CELDIS = 10.1 (12 time increments, simulation
time, 15 sec [57]), NSC = 4, NSR = NSL = 2, and NT = 128, concentrations
at early times and short distances are somewhat low, but the results are good
elsewhere. Thus, ELLAM compares favorably with explicit MOC3D which
needed 2401 time steps to satisfy stability criteria, and with implicit MOC3D
which required 241 time increments. In the high dispersion case, Dxx = 0.1,
or = 1.0cm. The analytical solution, and concentrations for CELDIS = 1
(121 time steps) and CELDIS = 10.1 (12 time steps) are plotted in figure 1.8.
The CELDIS = 1 results are very close to the analytical solution, and without
the oscillations produced by MOC3D at short distances. The solution using
fewer time steps has concentrations which are somewhat low near the inflow
boundary, and high near the outflow boundary.
27
Parameter Value
T T xx 0.01 cm2/s
Â£ 0.1
OiL 0.1 cm
Ctj'H = CH'JV 0.1 cm
PERLEN (length of stress period) 120 sec
vx 0.1 cm/s
vy = vz 0.0 cm/s
Initial concentration (C) 0.0
Source concentration (C1) 1.0
Number of rows 1
Number of columns 122
Number of layers 1
DELR (Ax) 0.1 cm
DELC (Ay) 0.1 cm
Thickness (b) 1.0 cm
Table 1.1. Parameters used in ELLAM simulation of solute transport in a
onedimensional, steadystate flow system.
28
Figure 1.7. Plots of concentration as a function of cell node for onedimen
sional flow with a constant velocity field and low dispersion. Shown are the
analytical solution (lowest graph), ELL AM results using CELDIS = 1 (121
time steps), NSC = 32, NSR = NSL = 2, NT = 128, and ELLAM results using
CELDIS = 10.1 (12 time steps), NSC = 4, NSR = NSL = 2, NT = 128 (upper
graph). Results for CELDIS = 1 are virtually identical to the analytical.
29
Figure 1.8. Plots of concentration as a function of cell node for onedimen
sional flow with a constant velocity field and high dispersion. Shown are the
analytical solution (lowest graph), ELL AM results using CELDIS = 1 (121
time steps), NSC = 32, NSR = NSL = 2, NT = 128, and ELLAM results using
CELDIS = 10.1 (12 time steps), NSC = 4, NSR = NSL = 2, NT = 128 (upper
graph). Results using CELDIS = 1 are virtually identical to the analytical.
30
Various other tests were performed, including runs with a retardation
factor of other than one. Concentration profiles in space at various times were
also evaluated. Results comparable to the above were obtained. Even in the
extreme case shown in figure 1.9 of CELDIS = 61 (two time steps), NSC =
4, NSR = NSL = 2, NT = 128, qualitativity good results were calculated,
except near the outflow face. This demonstrates the apparent robustness of
the method.
The low dispersion, no sorption problem was also used to test the
way the method handles nonzero decay. Figure 1.10 exhibits results for A =
0.01 sec"1, CELDIS = 1, NSC = 32, NSR = NSL = 2, NT = 128 which
are in excellent agreement with the analytical solution. As in the case of no
decay, NSC = 4 (not shown) yields a slightly low concentration near the inflow
boundary.
In all cases described above, the mass balance error was less than
0.001 percent, in contrast to methodofcharacteristics code, which is not mass
conservative. Explicit and implicit MOC3D, for example, yielded mass balance
errors of up to a few percent in some cases reported above.
1.5.2 Uniform, ThreeDimensional Flow
Next, ELL AM results were compared with the analytical solution
developed by Wexler [99] for threedimensional solute transport from a contin
uous point source in a steady, uniform flow field in a homogeneous aquifer of
infinite extent. The problem and analytical solution are described in detail by
Konikow and others in [71], and specific parameter values for the test case are
shown in figure 1.5.2.
Since the flow velocity is aligned with the grid and dispersive flux
crossterms are zero, this problem provides a test of the accuracy of dispersive
flux calculations in three directions [71]. It also offers a test of the source mass
algorithm, and its use in representing the effects of a specified flux boundary
condition.
Analytical results are shown in figure 1.11, and ELL AM results are
plotted in figures 1.12, 1.13, and 1.14 for the xy plane passing through the
31
Figure 1.9. Concentration vs. cell node plot with CELDIS = 61 (two time
steps), NSC = 4, NSR = NSL = 2, NT = 128.
Figure 1.10. Plots of concentration as a function of cell node for decay
constant A = 0.01 sec1. Shown are the analytical solution (lower graph)
and ELLAM results using CELDIS = 1 (121 time steps), NSC = 32, NSR =
NSL = 2, NT = 128.
32
Parameter Value
T T Jxx ^yy 0.0125 m2/day
Â£ 0.25
OiL 0.6 m
OC'j'H 0.03 m
CX.'j'V 0.006 m
PERLEN (length of stress period) 400 days
vx 0.1 m/day
0.0 m/day
Initial concentration (Co) 0.0
Source concentration (C1) 2.5106 g/m3
Q (at well) 1.0106 m3/d
Source location row 8, column 1, layer 11
Number of rows 30
Number of columns 12
Number of layers 40
DELR (Ax) 0.5 m
DELC (Ay) 3 m
Thickness (b) 0.05 m
Table 1.2. Parameters used in ELLAM simulation of transport from a con
tinuous point source in a threedimensional, uniform, steadystate flow system
33
point source. For 1.12, CELDIS = 7 (two time steps), NSC = NSR = NSL
= 4, NT = 16; for 1.13, CELDIS = 1 (14 time steps), NSC = NSR = NSL
= 4, NT = 4; and for 1.14, CELDIS 0.1 (134 time steps), NSC = NSL =
4, NSR = 8, NT = 16. MOC3D reports, [71] and [70] discuss the greater
spreading in numerical models than the analytical solution, and note that this
is partially explained by the application of the source over an entire grid cell
in the models, while the analytical solution portrays a point source.
Figure 1.12 can be interpreted as demonstrating that more than two
time steps are needed to adequately resolve dispersion. As seen in figure 1.13,
ELLAM results with 14 time steps accurately characterize the dispersive flux,
without the spreading upstream that is produced by MOC3D. With 134 time
increments, ELLAM produces yet less upstream spreading; so little, however,
that the solution manifests the numerical oscillations typical where concentra
tion gradient is too steep relative to grid mesh density.
Graphs in vertical planes parallel and perpendicular to the flow di
rection are comparable to the above, in terms of the proximity of the ELLAM
results to the analytical.
1.5.3 TwoDimensional Radial Flow
The ELLAM solution was compared to the analytical solution give
by Hsieh [61] to a radial flow and dispersion problem with a finiteradius well
in an infinite aquifer of twodimensions. There is flow from a single injection
well, with velocities inversely related to distance from the well.
The parameters for the problem are given in table 1.5.3. More details
about the problem are presented in [71] and [70]. Due to symmetry, the problem
can be modelled on one quadrant of the radial flow field.
The analytical solution is plotted in 1.15, along with four ELLAM
solutions in figures 1.16 1.19. Each ELLAM solution captures the salient
qualitative features of the analytical solution, with the exception of the high
concentration contour of 1.16, which has an articulated rather than a smooth
shape.
The two time step run graphed in figure 1.16 uses CELDIS = 75,
34
Figure 1.11. Concentration contours of analytical solution in the horizon
tal plane containing the solute source (layer 1) for threedimensional solute
transport in a uniform steady flow problem.
Figure 1.12. Concentration contours in the horizontal plane containing the so
lute source (layer 1) for threedimensional solute transport in a uniform steady
flow problem with CELDIS = 7 (two time steps), NSC = NSR = NSL = 4,
NT = 16.
35
5 10 15 20 25 30
Figure 1.13. Concentration contours in the horizontal plane containing the so
lute source (layer 1) for threedimensional solute transport in a uniform steady
flow problem with CELDIS = 1 (14 time steps), NSC = NSR = NSL = 4, NT
= 4.
5 10 15 20 25 30
Figure 1.14. Concentration contours in the horizontal plane containing the so
lute source (layer 1) for threedimensional solute transport in a uniform steady
flow problem with CELDIS = 0.1 (134 time steps), NSC = NSL = 4, NSR 
8, NT = 16.
36
Parameter Value
T = T *xx yy 3.6 m2/hour
Â£ 0.2
OiL 10.0 m
OtrpH Oirpy 10.0 m
PERLEN (length of stress period) 1000 hours
Q (at well) 56.25 m3/hour
Source concentration (C) 1.0
Number of rows 30
Number of columns 30
Number of layers 1
DELR (Ax) = DELC (Ay) 10.0 m
Thickness (b) 10.0 m
Table 1.3. Parameters used in ELL AM simulation of
twodimensional, steadystate, radial flow case.
Figure 1.15. Contour plot of analytical solution for two dimensional radial
flow problem.
37
Figure 1.16. Contour plot of concentration for run with two time steps using
CELDIS = 75, NSC = NSR = 4, NSL = 2, NT = 16.
38
NSC = NSR = 4, NSL = 2, NT = 16. The high concentration contour has
nonphysical oscillations, which cannot be entirely eliminated even with higher
NS and NT values. The 29 time step run shown in figure 1.17 uses CELDIS
= 5, NSC = NSR = 4, NSL = 2, NT = 4. Although the high concentration
contour is somewhat flat, it is a close approximation to the analytical solution.
The advisability of using higher NS values when modeling with many time
steps is illustrated in figures 1.18 and 1.19, where 563 time steps are used.
Parameter values NSC = NSR = 4, NSL = 2, NT = 4 produce a rather flat set
of contours, whereas NSC = NSR = 8, NSL = 2, NT = 4 yield a much better
match to the analytical solution.
Both versions of MOC3D produce close to analytical results for this
problem. For comparison, explicit MOC3D needed 596 time steps and 12 min
30 s of cpu time; implicit MOC3D used 282 time steps and 7 min 25 s of cpu
time; and ELLAM runs discussed above needed 2 min 18 s, 33 min 10 s, and
96 min 20 s, respectively, all on a UNIX workstation [57].
1.5.4 Initial Condition in Uniform Flow
A problem of threedimensional solute transport from an instanta
neous point source in a flow field with uniform velocity was also considered.
Wexler [99] presents an analytical solution for a continuous point source in a
homogeneous aquifer of infinite extent, which was modified to yield quantita
tive analytical results for the instantaneous problem.
This initial condition is represented numerically by a nonzero concen
tration in a single finitedifference cell. For MOC3D, this signifies a constant
concentration distributed evenly across a cell. For ELLAM, with its assump
tion of concentrations varying piecewise trilinearly between cell nodes, this
signifies initial mass in 27 cells. Furthermore, ELLAM is known to be ill
equipped to handle advection of fronts discretized with fewer than four nodes
(see section 1.6.5), so that ELLAM is expected to have difficulty propagating
even this dispersed point source.
However, a nonzero initial concentration in one cell, with a uniform
39
5 10 15 20 25 30
Figure 1.17. Contour plot of concentration for run with 29 time steps using
CELDIS = 5, NSC = NSR = 4, NSL = 2, NT = 4.
5 10 15 20 25 30
Figure 1.18. Contour plot of concentration for run with 563 time steps using
CELDIS = 0.25, NSC = NSR = 4, NSL = 2, NT = 4. Concentration maximum
is 1.019.
40
5 10 15 20 25 30
Figure 1.19. Contour plot of concentration for run with 563 time steps using
CELDIS = 5, NSC = NSR = 8, NSL = 2, NT = 4. Concentration maximum
is 1.0056.
41
velocity field, has been tested, and results are compared to the analytical solu
tion of a problem with an instanteneous point source. Also, to better test the
ELLAM code, this analytical solution has been allowed to develop over time,
and the smoothed concentration distribution is used as an initial condition for
ELLAM. The uniform flow problem with this second initial condition, is also
solved with two different velocity fields to test for grid orientation effects on the
dispersion calculations. Parameters for the problem are given in table 1.5.4.
For the case of flow in the x direction only, Vx = 1.0275, and Vy =
Vz = 0. For flow at 45 degrees to the grid, Vx Vy = 1.0275, and Vz = 0.
So the distance the center of the plume moves in the x direction is the same
for both cases, for equal simulation times, but the magnitude of the velocity
is greater for the flow at 45 degrees, introducing more dispersion. For all runs
reported below, CELDIS = 5 (six time steps), NSC = NSR = NSL = 4, NT
= 2. All graphs are concentration contours in the plane of the initial source
of solute, and of movement. In figure 1.20 the analytical solution using an
initial point source is plotted, and in figure 1.21, the ELLAM solution. The
numerical results show some spreading relative to the analytical, in both the
transverse and longitudinal directions. Increasing the number of time steps
does not completely eliminate the spreading, and causes some loss of peak
concentration, even with increased NS values. In contrast to both versions of
MOC3D, ELLAM results exhibit the symmetry of the analytical solution. In
figure 1.22 the analytical solution using an initial dispersed source is plotted,
and in figure 1.23, the ELLAM solution. These solutions resemble each other
more closely, although the ELLAM solution is still somewhat dispersed.
The analogous results with a dispersed initial condition are plotted
for flow at 45 degrees to the grid. In figure 1.24 is the analytical solution,
and in figure 1.25 the ELLAM solution is plotted. The characteristic symme
try of the analytic solution is captured by ELLAM, but there is longitudinal
spreading, and some distortion of the shape of the plume. This narrowing is
characteristic of a gridorientaion effect caused primarily by offdiagonal terms
of the dispersion tensor.
42
Parameter Value
T T Jxx ^yy 10.0 m2/day
Â£ 0.1
aL 1.0 m
OLj'H = OL'J'V 0.1 m
PERLEN (length of stress period) 90 days
vx 1.0 m/day
Vy = Vz 0.0 m/day*
Initial concentration at source 1 x 106
Source location column 11,
in transport row 36,
subgrid layer 4
Number of rows 72
Number of columns 72
Number of layers 24
DELR (Ax) 3.33 m
DELC (Ay) 3.33 m
Thickness (b) 10.0 m
* For flow at 45 degress to x and y axes, Vy = 1.0275m/day.
Table 1.4. Parameters used in ELLAM simulation of threedimensional
transport from a point source with flow in the xdirection and flow at 45
degrees to x and yaxes.
43
Figure 1.20. Contour plot showing log of concentrations in analytical soluion
of Dirac problem at t = 90. Concentration maximum is 25195.
Figure 1.21. Contour plot showing log of concentrations for spike initial
condition, and CELDIS = 5, NSC = NSR = NSL = 4, NT = 2. Concentration
maximum is 15160.
44
Figure 1.22. Contour plot showing log of concentrations in analytical soluion
of Dirac problem at t = 130. Concentration maximum is 14539.
Figure 1.23. Contour plot showing log of concentrations for dispersed initial
condition, and CELDIS = 5, NSC = NSR = NSL = 4, NT = 2. Concentration
maximum is 16909.
45
10 20 30 40 50 60 70
I I i I I I I
Figure 1.24. Contour plot showing log of concentrations of analytical solution
to Dirac problem at t = 130. Concentration maximum is 8645.
10 20 30 40 50 60 70
i11111r~
Figure 1.25. Contour plot showing log of concentrations for dispersed initial
condition, and CELDIS = 5, NSC = NSR = NSL = 4, NT = 2 at t = 130.
Concentration maximum is 8167.
46
20
10
0
0
Figure 1.26. Contour plot showing concentrations of a twodimensional finite
element solution of the Burnett and Frind problem. Contours shown are 0.1
to 0.9.
1.5.5 Constant Source in Nonuniform Flow
A problem having a constant source of solute over a finite area at the
surface of a homogeneous aquifer with boundary conditons producing nonuni
form flow was introduced by Burnett and Frind [21]. They used an alternating
direction Galerkin finite element method to solve the flow and transport equa
tions in two and three dimensions. This Burnett and Frind problem has been
used as a test case to evaluate both versions of MOC3D, with the finite ele
ment results as the benchmark. It is a complex problem, including a variable
dispersion coefficient. Konikow and others [71] present a detailed desciption
of the problem geometry including a discussion of the differences between the
finite element grid and the cellcentered finitedifference grid. Parameters used
for the ELL AM simulation are given in table 1.5.5.
The ELL AM simulator was used on the two dimensional case, using
CELDIS = 30 (seven time steps), NSC = NSR = NSL = 4, NT = 32. The plume
calculated by Burnett and Frind is shown in figure 1.26, and the ELLAM results
in figure 1.27. The contours are very similar, although the high concentration
contour from ELLAM doesnt extend as far downgradient as that of Burnett
and Frind, while the lowest concentration curve extends further. Nevertheless,
the ELLAM solution provides a closer match to the Burnett and Frind contours
than do MOC3D results using 381, 1901, or 4218 time steps [70]. Increasing
the number of time increments yields a solution with a greater downgradient
extent, but still short of the Burnett and Frind results.
Fourteen rows were added to the twodimensional grid in order to
expand it to three dimensions. Transport results for a vertical plane at the
47
Parameter Value
K 1.0 m/day
Â£ 0.35
OIL 3.0 m
Oi'jpH 0.10 m
Ci'j'V 0.01 m
PERLEN (length of stress 12,000 days
period)
Q (at well) 56.25 m3/hour
Source concentration(C") 1.0
Number of rows 1
Number of columns1 141
Number of layers1 91
DELR (Ax) 1.425 m
DELC (Ay) 1.0 m
Thickness (b) 0.22220.2333 m
1 One row and layer were allocated to defining boundary conditions, so
concentrations calculated in only 140 columns and 90 layers were used for
comparision.
Table 1.5. Parameters used for ELLAM simulation of transport in a verti
cal plane from a continuous point source in a nonuniform, steadystate, two
dimensional flow system (described by Burnett and Frind, 1987).
48
Figure 1.27. Contour plot showing concentrations of ELLAM solution to
twodimensional Burnett and Frind problem using CELDIS = 30, NSC = NSR
= NSL = 4, NT = 32. Contours shown are 0.1 to 0.9.
49
20
10
0
r uu (a) 3D FiniteElement Model
 1 (1 Contours: 0.1 to 0.9
i i i i '  i .. i . > : i
0 40 80 120 160 200
Figure 1.28. Contour plot showing concentrations of a threedimensional
finite element solution of the Burnett and Frind problem. Contours shown are
0.1 to 0.9.
middle of the plume are plotted for runs with lower and higher dispersion.
Figure 1.28 shows Burnett and Frind results for the low dispersion case where
a.Tv = 0.01m and olth = 0.1m, while ELLAM are given in figure 1.29. The
ELLAM plume extends slightly further downstream with low concentrations
of solute than the Burnett and Frind. Again, the ELLAM solution is a closer
match to the Burnett and Frind than any MOC3D solution.
With vertical transverse dispersivity increased by a factor of ten so
that ctTv = olth = 0.1m, the results of simulation were again compared. The
ELLAM solution using CELDIS = 30 had noticably low concentrations near
the source. Plotted in figures 1.30 and 1.31 are the Burnett and Frind solution
for higher dispersion, and an ELLAM solution using CELDIS = 21 (10 time
steps), respectively. The contours are in close agreement.
1.6 OneDimensional ELLAM
Consider a onedimensional problem with e = Rf = 1, and D > 0,
and assume, for simplicity, a uniform grid and v = constant > 0. Also, forego
the use of approximate test functions, that is, let
where
wjik(x,y,z) = f{x)g(y)h{z)
f(x)
0
< 1 \
and g = h = 1.
50
Figure 1.29. Contour plot showing concentrations of ELLAM solution to low
dispersion Burnett and Frind problem using CELDIS = 30, NSC = NSR =
NSL = 4, NT = 32. Contours shown are 0.1 to 0.9.
20
.10
o
p (a) 3D FiniteElement Model ,
Contours: 0.1 to 0.9 
0 40 80 120 180 200
Figure 1.30. Contour plot showing concentrations of a threedimensional
finite element solution of the Burnett and Frind problem with high vertical
transverse dispersivity. Contours shown are 0.1 to 0.9.
51
Figure 1.31. Contour plot showing concentrations of ELL AM solution to
high dispersion Burnett and Frind problem using CELDIS = 21, NSC = NSR
= NSL = 4, NT = 32. Contours shown are 0.1 to 0.9.
52
Neglecting sources and inflow, which do not affect stability, the system
of ELL AM cell equations can be expressed
(A + M)cn+1 = Bcn (1.25)
where A is the coefflcent matrix for integration of concentration at the new
time level; M is the matrix of diffusion coefficients obtained from a backward
Euler in time and centered differences in space approximation to the diffusion
term of the continuous equation; B is the coefficient matrix for concentration
from the old time level, incorporating an integration to determine mass in a
cell, followed by a shift operation to represent advection; and c is the vector of
nodal concentration values at the time level denoted by the superscript. Due
to the assumption of a uniform grid, dependence of the numerical integrals on
cell length has been divided out.
Using the trapezoidal rule to integrate the unknown piecewise lin
ear concentration function, and using zero Dirichlet boundary conditions for
simplicity, the A matrix is tridiagonal with
X = [1,6,1] (126)
denoting the sub, main, and superdiagonal entries, respectively.
The matrix M of centered difference coefficients at the new time level
is
M = At k[l, 2, 1], (1.27)
where for diffusion coefficient D > 0, k = Note that A + M is symmetric
positive definate for all k > 0, so (A + M)_1 exists.
The matrix B depends on Courant number, defined by Cr =
The matrix considered here is for a constant Courant number, 0 < Cr < To
construct a B which aportions the correct mass to each cell at tn+1, integrate
exactly at the old time level using integration points at cell centers and pre
images at the old time level of cell boundaries at the new time, as shown in
figure 1.32. Denoting by x* the position at the old time level of a point x at
53
i1 i i+1
xI/2 x+1/2
Figure 1.32. Mass at old time level is integrated exactly using cell centers
and preimages of cell boundaries as integration points, and applying the trape
zoidal integration rule.
tn+1, and interpolating concentration values between nodes, the trapezoid on
the left has area,
(1r+CrAx)
'C(x*_ k) + Q
and the area on the right is,
/Ax
VT~ ~~
CrAx
Ci + C(
Since
C(x*_i) = (i + Cr) + (i Cr) Ct
and
C(x+i) = (l + Cr) ft + (i Cr) Ci+U
the mass in cell i at the new time level is
M(d + Cr)(J + Â¥)Ci
+ ((5+Cr)(ff) + aCr)(! +
+ ((5C'r)(Jf)Ci+1'
))Ci
Thus, disregarding boundary conditions, the B matrix is tridiagonal with
(1.28)
1 Cr Cr2
B = [1,6,1] + [1,0, 1] + [1, 2,1],
denoting the sub, main, and superdiagonal entries, respectively. In the case
of Cr = , B is singular. The more general case of arbitrary, still constant,
Courant number is considered analogously: the Courant number is decomposed
into the sum the nearest integer, and a fractional part with \ fractional part\ <
In (1.28), Cr in the second two terms is replaced by the fractional part, and
the entire matrix sum is diagonalshifted by minus the integer nearest to the
Courant number. For simplicity, only the case 0 < Cr < \ is treated here.
The integration points and rule used in the above construction differ
from those in the threedimensional ELL AM implementation. Application of
the threedimensional code to a onedimensional problem on a uniform grid
55
with Courant number a multiple of A., n = 1,2,3,... likewise tracks and inte
grates advected mass exactly, and yields a solution identical to that produced
by A~lB, up to treatment of boundary conditions.
Convergence of the method (1.25) will follow from its consistency and
stability. The exact solution, c(,t) satisfies
(A + M)c(, tn+1) = Bc(, tn) + Atec, (1.29)
where ec denotes the vector of truncation (consistency) errors at each node.
Letting en denote the error cn c(,tn), and subtracting (1.29) from (1.25),
we have the error propogation
(A + M)en+l = Ben + A tec.
Because A + M \I is nonnegative definate, this implies
en+1 = (A + M)~1Ben + A tec,
with a different constant incorporated into ec. This leads to
en+1 = Â£((A + M)~lB)jAtec + ((A + M)~1B)n+1e,
3=0
where e denotes any initial error. For stability we then require,
((A + M)15)n < K
in some norm, for 0 < nAt < T, where the constant K is independent of At,
and T is simulation time. Then en+1 < 0(iC)ec, and consistency yields
convergence as Ax, At 0.
Stability depends on the properties of the matrix (A + M)~lB which
operates explicitly upon a error vector at the old time level to yield error at
the new time.
1.6.1 Limiting Case of No Advection
In the limiting case of nonzero diffusion, but no advection,
A + M=[\K6 + 2k,1lk]. (1.30)
56
(1.31)
The righthand side matrix B becomes
S = i[l,6,l],
since with no advection, Cr = 0.
Lemma 1.1 In the limiting case of no advection, the method given by equa
tion (1.25), using (1.30), and (1.31) is consistent, with error of order,
A tce = 0(Ax2At + At2).
Proof: We use standard finite difference arguments.
First note that integrating the Taylor expansion of a function / gives,
l'xi+1/2
[ + f(x) dx
Jxil/2
= f(Xi) + f(Xi)(X ~ Xi) + y(xi)(x ~ xi)2 + 0((x Xif dx
= Axf(xi) + 0 + (a: Xi)2 dx + 0
= Axf(xi) + ^f"{xi) + 0(Ax3),
where it is useful to separate the and 0(Ax3) terms. Then, we have
fÂ£+1 fg+% ct{x,t) dxdt
= Ax J?+1 ct(xi, t) + ct{xi, t) + 0( Ax2)
= Ax C1 Ct(Xi, t) + Â£f fr*1 + 0(Ax2^ + 0(Ax2) dt
= Ax (c"+1 cf + l [c^1 cj*_! 2(c+1 c?) + dg? cf+1
+0(AÂ£Ax4) + 0(AtAx2))
= Ax (Kct+i1 cU) + Wrl ~ c?) +  c?+1) + 0(AtAx2)) .
Also,
rn+1
,{x, tn+1) J cxxt(x, s) ds
Ax2 d2
dt
Cxx  CX.
rtn+l
Cxi(Xi,tn+1) +Cxxx(Xi,fn+1)(:r Xi) + 0(Ax2) j cxxt(x, s) ds
<Â£!? 2 c+1 + eft1
(Ax)2
+ csa;a;(xi,Â£n+1)(x Xi) + 0(Ax2)
rt"+1
J CxxtijX, s)
ds.
57
Integrating then gives
and
L
"i+1/2
~(x, t) dx
i1/2
Ax
_ rx,
JXi
,n+l_
2c"+1+c^+_A
+0(Ax3)
*.+1/2 rtn+1
xi1/2
(Ax)2 J
ir cxxt(x, s) dsdx;
Dcxx(x,t) dtdx
= + 2c"+1 cm) + /t"+1 /f+1 c**t(z, s) dsctedf
 Az^cft1 + 2c+1 c^1) + 0(Ax3At) + 0(AxAt2).
Combining terms gives the estimate,
/t+1 ct(. t) dxdt //;+1 Dc^rr, i) dxdt
= Aa: (Kcft1
A^(c'!+1 + 2c?+1 c^1) + 0(A:r3Ai) + 0(AzAt2).
Dividing by Ax to reflect (1.25) completes the proof.
Lemma 1.2 In the case of no advection, the method given by (1.25) is uncon
ditionally stable in the l2 norm, independent of the meshsize.
Proof: Write (1.30) as
A + M = I+(k 1) [1,2, 1]
and (1.31) as
B = I ^[1,2,1]
where / is the identity matrix. The matrix (A + M)_1B has eigenvalues
1 f A
1 + (A: )A
58
where A is an eigenvalue of [1,2,1], Since 0 < A < 4, (A + M)~XB has
nonnegative eigenvalues less than or equal to one for any k > 0. Because
(A + M)~lB is normal, this gives \\{A + M)~l B\\% < 1. Therefore, the method
is stable.
Convergence follows directly from consistency and stability.
Theorem 1.3 In the case of no advection, the method given by (1.25) is con
vergent, independent of the meshsize.
It is interesting to note that for k > , A~lB satisfies a discrete
maximum principle. For 0 < k < , however, the operator need not satisfy a
maximum principle; consider, for example, 3x3 A~XB with k = ~. Here,
( 2721 14 1 \ ( 6 i
4 14 8 8
14 196 14 1 6
679 8 8
V 1 14 195 0 1 8
/ 16130 1559 1 \
A 112 112
91 574 56
679 4 4 4
56 578
V 1 4 4 /
0 N
i
8
The absolute row sum of the second row is . Since this is greater than 1,
the matrix does not satisfy a maximum principle.
1.6.2 Limiting Case of No Diffusion
Consistency analysis for the purely hyperbolic problem has been done
by Stynes and Russell [84]. Consistency error of 0((A:r)2 At) for the interior
of the domain, and 0(Ax + At) near an inflow boundary is established.
Theorem 1.4 In the case of no diffusion, the method given by (1.25) is stable
in the l2 norm, for any fixed meshsize.
Proof: Since for no diffusion, k = 0,
(.A + M)~lB =A~XB
= I + A~1 (f[l,0,l] + Â¥[l,2,l])
= / + f A1 ([1,0, 1] + Cr[ 1, 2.1]).
59
Then, since 0 < Cr <
\\(A + M)~lB\\2 =\\A~lB\\2
< 1 + %\\A% ([1,0, 1]2 + CV[1, 2,1]2)
< 1 + f A12(2 + 4Cr)
< 1 + 2CVA 112
<1 + 4^
= 1 + O(At) for fixed Ax.
We have
(A1B)^[2 <\\AlB\\f
< (l + 0(At))&
< eA/,
where T is the total simulation time. Thus the method is stable for fixed
meshsize.
Convergence again follows from consistency and stability.
Theorem 1.5 In the case of no diffusion, the method given by (1.25) is con
vergent, for fixed meshsize.
1.6.3 Infinite Spatial Domain
An infinite spatial domain is used in order to eliminate boundary
effects. Here, both the lefthand side (1.26) and righthand side (1.28) operators
are Laurent and Ax > 0. For an introduction to the theory of Laurent matrices,
see [14], which serves as a reference for the following discussion. We will show
that in the case of an infinite spatial domain, that A~XB is a Laurent matrix
with norm of one, and note properites of its spectrum.
Given a sequence {2;_}^i_0O of complex numbers, a doublyinfinite
60
matrix with zn along the nth diagonal is called a Laurent matrix:
... Zq ... Zi Z1 Zq Z2 z1 Z3 z2 z4 ... Z3
... %2 Zi Zq Z1 Z2
... z3 Z2 Zl ZQ Zi ...
... ZA Z3 z2 Zl Zq ... J
The following theorem is wellknown:
(1.32)
Theorem 1.6 (Bottcher and Silberman [14], Theorem 1.1) The Laurent ma
trix (1.32) generates a bounded linear operator on l2(Z) if and only if there is a
function z G L(T) such that {zn}^_^ is the sequence of Fourier coefficients
of z:
1 /2tt
= = / z{eme)e~m6d0 (n G Z),
V27T Jo
where Z means the integers, and T is the complex unit circle.
Should a Laurent matrix indeed represent a bounded linear operator,
denote the matrix and the operator by C(z). Let L := L(T) and L2 :=
L2(T).
Laurent matrices represent multiplication operators on L2 with re
spect to the orthonormal basis
AnO
nÂ£Z
Define
M(z) : L2 > L2, f^zf.
This multiplication operator is bounded if and only if z G L00, in which case,
ll^WII = Iklloo
(1.33)
61
Consider the operator which maps a function to the sequence of its
Fourier coefficients,
$ : L2 ) /2(Z), / {fn}nez
The operator $ is bijective and
(134)
for / G L2. Then
C{z) = {z)$\
(1.35)
and properties of C can be established by determining the analogous property
of M.
Theorem 1.7 For the case of no diffusion and an infinite spatial domain,
A~lB = C(a~lb),
(1.36)
with
a lb
3 + cos(0)
Cr ^ cos(0) + Cr eld + ^ Cr2
and Â£(a 1&) = 1 for any Courant number Cr G [0, ].
gL,
(1.37)
Proof: The continuous function a G L(T) with Fourier coefficients com
posing the diagonals of (1.26) is given by
W =57K(e_a, + 6 + e")
= 472^(3 + COS(^(('
The continuous function b G L(T) with Fourier coefficients composing the
diagonals of (1.28) is given by
 vfc '(f + T)^+(Cr2) + ( l+f+ Â£y)e;
_ 1 (Cr ^ cos(0) + Cr el6 +  Cr2
(1.38)
62
We thus have Â£(a) and Â£(&). Furthermore, since zero is not in the range of a,
the inverse of Â£(a) is Â£(a_1), (see [14], Theorem 1.2), where
4\/2tr
a~l =
3 + cos (9)'
From (1.35), it follows that
Â£(2i)Â£(.z2) = C{ziz2) for all Zi,Z2 G Lc
(1.39)
(1.40)
So Â£(a x)Â£(&) = Â£(a x&) Multiplying continuous functions (1.39) and (1.38),
we get the bounded linear Laurent operator on an infinite domain,
AXB = Â£(a_16),
with a~lb given by (1.37).
From (1.35), (1.34), and (1.33), it follows that
Â£(F) = IMloo for all zeL.
(1.41)
Thus, to determine the norm of Â£(a 16), we find the maximum mod
ulus of the continuous function a~lb. For
p = (4Cr2l)2,
q = 2(4CV2 + 1)(3 4Cr2),
r = 16CV4 8Cr2 + 9,
we have
Then
Ui2 pcos2(9) + qcos(9) + r
b\ =
(3 + cos(0))2
= (rfsfe Pr 3, + ( + 6p) 008(9)]
~w
JsSfe (04Cr2(2Cr2 1)) (cos(9) 1)
= 0
for 9 such that sin(0) = 0 or cos(0) = 1. Thus the critical points for a lb\2 are
0 and 7r when 9 G [0, 27t). By (1.37), for 0 = 0,
a 1b = 1,
63
and for 9 = 7r,
a~1b = 1 4Cr2 G [0,1].
Thus the maximum modulus occurs at 9 = 0 and is equal to 1. Then (1.41)
gives Â£(a_16) = 1 for any Courant number Cr G [0, ], where the operator
norm is that induced by the l2 norm.
For an element z in a Banach algebra, define the spectrum of z as
the set
spz := {A Â£ C : z Ae is not invertible in the algebra},
where e is the algebra identity element.
For z G L, further define the essential range of z as the set
7Z(z) := {A G C : meas(i G T : z(t) A < e}) > 0 for every e > 0},
where meas means the Lebesgue measure. The essential range of z is the
spectrum of z as an element of the Banach algebra L, and equivalently,
the spectrum of the multiplication operator M(z) on L2 as an element of the
Banach algebra of bounded linear operators on L2. It then follows from (1.35)
that spÂ£(z) =7Z(z), whenever z G L.
We can now consider features of the spectrum of C(a~lb) as an element
of the Banach algebra of bounded linear operators on the complex Banach space
l2. We note that this definition of spectrum coincides with the usual definition
for the spectrum of a linear operator defined on a domain in a complex, normed
space (see Kreyszig [72] pg 397).
Since a_16 = 1, the spectrum of C{a~lb) is contained in the closed,
complex unit disk. Since a~lb is not (globally or locally) a complex constant,
C(a~lb) has no eigenvalues. The entire spectrum is seen to be continuous.
Also, 0 G sp(o_15) only for a Courant number of
Lemma 1.8 For a~lb given by (1.37) and Cr G [0, ],
0 G sp(a_16) if and only if Cr =
64
(1.42)
Proof: We have a 1b = 0 if and only if
0 = cos(0) (Cr2 Cr + ) + Cr eiB Cr2 + Â§
= (cos(0) 1 )Cr2 + (el9 cos{9))Cr + \ cos(9) + f.
For 0 ^ 0, using the quadratic theorem,
^ _ cos(0)ez6y/e2ie 2cos(0)e,0+cos2(0) (cos(0) l)(cos(0)+3)
= 2(cos(0)l)
cos (8)e'e 2 sin ()
4sin2(f) ''
Since Courant number is real, this is satisfied if and only if,
zsin(0) = 0,
so sin(0) must be identically zero. Since 6 G [0, 27t), this means 0 = 0 or 0 = tv.
If 0 = 0, substitution into (1.42) yields 0 = 1, hence no solutions. For 9 = tt,
we have Cr
1.6.4 Stability Independent of Meshsize
In Section 1.6.2, stability of the method was demonstrated for a fixed
meshsize on a finite spatial domain with zero Dirichlet boundary conditions.
A bound of A_1jB[ < 1 using the l2, A, or other norm equivalent to 
2 independent of matrix size would guarantee stability independent of Ax.
Shown in table 1.6.4 are MATLAB calculations, for various Courant numbers
and problem sizes, of the Euclidean and A norms of A~lB. Here, we consider
both the tridiagonal A_1B, and A~XB constructed with A^i = An>n =  to
reflect Neumann boundary conditions. In all examples,  A~1B\\a is bounded by
one for the tridiagonal case corresponding to the Dirichlet boundary condition.
Numerical study of the eigenvalues of A_1B, with A as in (1.26) and
B given by (1.28), suggests the moduli of the eigenvalues of this matrix are
bounded by one for any order of the matrix. If this is indeed the case, and the
inequality is strict, then the existence of a matrix norm with a value less than
or equal to one is assured:
65
71 tridiagonal U^b\\ Alti An,n r U'bu tridiagonal \\A~lB\\a <4.1,1 <4n,n = s
5 10 50 100 200 300 0.99798782758721 0.99984912202045 0.99999971570515 0.99999998183770 0.99999999885218 0.99999999977242 1.02731655221192 1.02646524138853 1.02646524324671 1.02646524324671 1.02646524324671 1.02646524324671 0.99367278911431 0.99945791347173 0.99999884309114 0.99999992482548 0.99999999520802 0.99999999904715 1.04572629392131 1.04220454113547 1.04216137329937 1.04216137329937 1.04216137329937 1.04216137329937
5 10 50 100 200 300 0.99911499209546 0.99998979236483 1.00002196618037 1.00002196618549 1.00002196618549 1.00002196618549 0.99811477727258 0.99994432525461 1.00001117880143 1.00001117886019 1.00001117886019 1.00001117886019 0.99578523861235 0.99962924605506 0.99999920055886 0.99999994803407 0.99999999668714 0.99999999934125 0.99365010768407 0.99948444126033 0.99999912922811 0.99999994570318 0.99999999661270 0.99999999933138
5 10 50 100 200 300 1.00154201321822 1.00154201699366 1.00154201699366 1.00154201699366 1.00154201699366 1.00154201699366 1.00009454693413 1.00009495089198 1.00009495089972 1.00009495089972 1.00009495089972 1.00009495089972 0.99998843721148 0.99999896202120 0.99999999774353 0.99999999985328 0.99999999999065 0.99999999999814 0.99996021692627 0.99999809688996 0.99999999745662 0.99999999984424 0.99999999999036 0.99999999999810
5 10 50 100 200 300 1.00023639568756 1.00023639647505 1.00023639647505 1.00023639647505 1.00023639647505 1.00023639647505 1.00003687264714 1.00003687289972 1.00003687289972 1.00003687289972 1.00003687289972 1.00003687289973 0.99999981047593 0.99999998298590 0.99999999996301 0.99999999999760 0.99999999999985 0.99999999999997 0.99999933999554 0.99999996868398 0.99999999995829 0.99999999999745 0.99999999999984 0.99999999999997
5 10 50 100 200 300 1.00006026622643 1.00006026665609 1.00006026665609 1.00006026665608 1.00006026665608 1.00006026665609 1.00001031695580 1.00001031868571 1.00001031868571 1.00001031868571 1.00001031868571 1.00001031868571 0.99999998815467 0.99999999893661 0.99999999999769 0.99999999999985 0.99999999999999 1.00000000000000 0.99999995869556 0.99999999804191 0.99999999999739 0.99999999999984 0.99999999999999 1.00000000000000
5 10 50 100 200 300 0.99799236275199 0.99984936308447 0.99999971662918 0.99999998194544 0.99999999886528 0.99999999977629 1.02487950575327 1.02391532781825 1.02391531969142 1.02391531969142 1.02391531969141 1.02391531969141 0.99368383249979 0.99945833836339 0.99999884356986 0.99999992485559 0.99999999520993 0.99999999904753 1.04147443273890 1.03829256968434 1.03823646002098 1.03823646002098 1.03823646002098 1.03823646002098
Table 1.6. Sections represent Courant numbers, Cr = 1/2, 1/3, 1/64,
1/500, 1/2000, 99/200, respectively.
66
Lemma 1.9 (Horn and Johnson [60], Lemma 5.6.10) Let A be an n x n com
plex matrix and e > 0 be given. There is a matrix norm   such that
q(A) < A < g(A)+e.
The equivalence of the unknown norm to the l2 norm may be related to the
order of the matrix, however.
Various attempts to prove the bound on the moduli of the eigenvalues
have failed: Bounds on the numerical radius have been studied (see Gustafson
and Duggirala [49]). The generalized eigenvalue problem, Ax. = XBx, has been
solved to express an eigenvalue in terms of a component of its (scaled) eigen
vector. The inverse matrix of (1.26) has been expressed explicitly (see [73]),
reformulated in terms of exponentials, and applied to B of (1.28). In no case
have eigenvalue bounds of one been forthcoming.
Sharp bounds on
\\AxB\\2 = q((A~1B)t A~XB)
and
WA'BWa = q(A*BtA~1BA~2) = q((B A~l)T A~l B),
where g is the spectral radius, have not been analytically established.
It is easy to see that the matrix equation Acn+1 = Bcn, with A given
in (1.26), and B defined by (1.28), results from a forward difference in time,
centered difference in space discretization of the equation,
Cr2{ Ax)2
C l C Ct^XC
C't ng xxt 1
2A t
~ (XXX, 0,
where the coefficients are taken to be constant. With the opposite sign on the
spacetime derivative, this would be an equation of the Sobolev, or pseudo
parabolic, type. Disregarding the Cx term, this continuous operator has eigen
values greater than 1.
1.6.5 Numerical Dispersion and Oscillations
It is wellknown (cf. [83]) that ELLAM methods, in general, require
sufficient meshdensity for about four grid nodes on a front for accurate model
ing of advected concentration. (Compare with Eulerian methods where number
67
1
_0 511111111111
20 40 60 60 100 120 140 160 180 200
Figure 1.33. Initial and advected peaks with Cr = 10.
of nodes on a front must increase as the front sharpens, as discussed in Sec
tion 1.1.) Since an arbitrary concentration distribution is a linear combination
of chapeau functions, the behavior of an initial condition of nonzero concentra
tion at one node only under repeated application of the nodiffusion operator
was studied. Qualitative features of the solution are independent of the problem
size. The following represent results for the onedimensional problem currently
under discussion.
In each of figures 1.33, 1.34, 1.35, 1.36, and 1.37, an initial condition
of one at node four is plotted. Also plotted is the solution after time T = 180
under pure advection, with v = 1 and Courant number as indicated. With
time discretization to yield an integer Courant number, the peak is advected
without distortion as shown in figure 1.33. If Cr = , symmetric numerical
dispersion is evident, as shown in figure 1.34. An nonsymmetric, oscillatory
graph is produced for Cr 6 (, ). As Cr > 0, the maximum decreases and
the oscillations become more extensive. This behavior is shown is figures 1.35
and 1.36. The train of oscillations extends primarily to the right of the max
imum for Cr < 0 and to the left for positive Courant number. Figure 1.37
illustrates the advisability of using large time steps: fewer applications of the
68
Figure 1.34. Initial and advected peaks with Cr
Figure 1.35. Initial and advected peaks with Cr = \.
69
Figure 1.36. Initial and advected peaks with Cr = ^
Figure 1.37. Initial and advected peaks with Cr = 10^.
70
operator with noninteger Courant number results better preservation of the
initial contour.
This behaviour can be considered in light of the eigensystem of the re
lated operator S(A~1B)& in cases where Cr = intlger Here, S is the diagonal
shift by 1 operator,
s = [1,0,0].
Disregarding boundary effects, application of (S(A~1B)'^)T results in a con
centration distribution with the same form as the advected solution, but sta
tionary. The shifted matrix acts repeatedly on a linear combination of its eigen
vectors, and the solution evolves accordingly. Ideally, if the initial condition
were advected without distortion, this matrix would be the identity. Consider
the 10 x 10 problem with an initial condition of one at node seven, v = 1,
and Cr = . All eigenvectors of S(A~1B)<^ are plotted in 1.38. The place
ment of entries with largest modulus near the beginning of any eigenvector is a
feature that characterizes matrices (
seen for ((A~1B)'&:S, which first shifts, then advects; and also (S(i)(A~lB)'^:)
and ((A1B)s?
corresponding to the maximum eigenvalue is the smoothest vector, with its
minimum at node three. Eigenvalues, and coefficients of their respective eigen
vectors to produce the initial condition are shown in table 1.6.5. (The matrix
is singular, due to the shift.) The initial condition holds large contributions
from eigenvectors with small eigenvalues, so dimunition and change of shape
of the solution is to be expected under repeated matrix application. Indeed,
this is observed. The initial condition is plotted along with the solution after
different elapsed times (number of matrix applications). Shown in figure 1.39
are graphs for T = 1, T = 10, T = 100, and T = 1000, corresponding to
the curves with highest, second highest, next, and smallest peaks, respectively.
The eigenvector of the maximum eigenvalue is seen to emerge as predominate.
As the size of the problem increases, the system of eigenvectors be
comes increasingly illconditioned. The entries with largest modulus remain
near the beginning of any eigenvector. It is not possible to approximate the
71
Figure 1.38. Eigenvectors of 10 x 10 S(A lB)cr with Cr =
eigenvalues
0.02512380295106
0.14348830209936
0.34537563664597
0.56442393631698
0.74559736660795
0.87006283785412
0.99708068200960
0.94405857162778
0.98197582016563
0
coefficients
11.08455231527136
17.99738031282700
6.21949567802237
5.06404951515980
6.99680498309164
2.97173348456930
0.76361709637109
0.71286631865321
1.50653675213281
0
Table 1.7. Eigenvalue and coefficient of the respective eigenvector in eigen
vector decomposition of spike initial condition at node seven on 10 x 10 grid.
72
Figure 1.39. Initial condition and results using 10 x 10 (S{A 1B)cr)T with
T = 1,10,100,1000 and Cr =
73
initial condition very well without using many of the eigenvectors. The largest
eigenvalues approach one (from below) in modulus. The eigenvector corre
sponding to the maximum eigenvalue remains the smoothest eigenvector.
A system with 40 unknowns is still small enough to find reasonable
approximations to the coefficients for an eigenvector expansion of the initial
condition, in spite of illconditioning of the matrix of eigenvectors. The ini
tial condition and the results of 1,10,100,1000, and 100000 applications of
S{A~lB)'^: are shown is figure 1.40. Again, the emergence of the dominant
eigenvector can be traced. Analogous results are plotted in figure 1.41 for a
system of order 200. The matrix of eigenvectors is too illconditioned for MAT
LAB to find coefficients for an eigenvector expansion of the initial condition,
but the eigenvalues are numerically clearly distinct, so presumably a full set
of eigenvectors exists. Figure 1.41 would then be understood as showing the
emergence of eigenvector components of eigenvalues close to one. The single,
smooth mode visible for T = 10000 in the smaller systems still is not isolated
in the 200 x 200 system. The case of Cr = \ is anomolous: The nonzero
entries of the eigenvector are periodically distributed evenly across the domain.
Next shown are graphs of the initial condition, and the concentration
profile after T 1, and T = 100 for Cr =  for the 40 x 40 system with
two (figure 1.42), three (figure 1.43), four (figure 1.44), and five (figure 1.45)
nodes across a front or peak. Superposition of the oscillatory contours from
the spikes composing the smoothed initial condition yields an advected concen
tration profile much closer to the original for the single spike initial condition.
The peak maximum initially may increase slightly with respect to the initial
condition. Still, with enough time steps, the peak tends to smooth, oscillate,
and be displaced upstream, as would be expected considering the underlying
eigensystem. Thus, even with welldiscretized fronts, using few (large) time
steps may be advisable with respect to the advection algorithm.
Note that the forgoing discussion applies exclusively to the consid
eration of finer time discretization, without accompanying refinement of the
spatial grid. In figures 1.46 and 1.47 are shown superimposed plots illustrating
74
Figure 1.40. Initial condition and results using 40 x 40 (S(A 1B)cr)T with
T = 1,10,100,1000,10000 and Cr =
Figure 1.41. Initial condition and results using 200 x 200 {S(A 1B) cr )T with
T = 1,10,100,1000,10000 and Cr =
75
1.2 
)
40
Figure 1.42. Initial condition with 2 nodes on a front, and results using
40 x 40 (S(A~1B)'^)T with T = 1,100 and Cr = Â§.
Figure 1.43. Initial condition with 3 nodes on a front, and results using
40 x 40 (S(A~1B)^:)T with T = 1,100 and Cr =
76
1.2
Figure 1.44. Initial condition with 4 nodes on a front, and results using
40 x 40 with T = 1,100 and Cr = f.
Figure 1.45. Initial condition with 5 nodes on a front, and results using
40 x 40 (S,(A1Â£)^r)T with T = 1,100 and Cr = J.
77
the convergence of the method as both Ax and At approach zero. A spike
initial condition on a 5 x 5 grid is advected two cells under a uniform velocity
field. The grid is refined four times, each time by a factor of three. Simulation
time and Courant number are constant for each figure, meaning a finer time
discretization with each spatial refinement. Increasing height of peak height
corresponds to progressively more refined grids of 5 x 5, 15 x 15, 45 x 45,
135 x 135, and 405 x 405.
1.7 Conclusion
The ELLAM algorithm for solution of the advectiondiffusion equa
tion can simulate the transient, threedimensional transport of a solute subject
to decay and retardation. The accuracy of the ELLAM numerical results were
tested and evaluated by comparision to analytical and numerical solutions to a
set of test problems. These tests indicate that ELLAM can accurately simulate
threedimensional transport and dispersion of a solute in flowing ground water.
The method appears robust, with demonstrated stability in a variety of test
situations. It compares favorably to the method of characteristics codes used
as benchmarks, in some cases. To avoid nonphysical oscillations and loss of
peak concentrations, care must be taken to use a grid with sufficient mesh den
sity to adequately resolve sharp fronts. ELLAM is globally and locally mass
conservative, and can provide good solutions using large time steps.
1.8 Further Research
With regard to an analysis of the onedimensional ELLAM meth
od (1.25), a proof of stability independent of the problem size remains to be
accomplished. Incorporation of various boundary conditions into a stability
proof is another direction for further effort. Generalization of convergence re
sults to higher dimensional problems remains to be done. With respect to an
ELLAM implementation, the possibility of specifying different spatial subdis
cretization for mass tracking (NS values) in different parts of the computa
tional grid, might result in greater efficiency of computation. Incorporation of
a backtracking scheme to identify the preimage at the old time level of a finite
78
1.2 
Figure 1.46. Superposition of results showing initial condition and advected
peak using Cr = Grids are 5x5, and 4 refinements, each by a factor of
three. Initial condition is a spike on 5 x 5 grid. Advected peaks show increasing
height with decreasing Arc.
Figure 1.47. Superposition of results showing initial condition and advected
peak using Cr = . Grids are 5x5, and 4 refinements, each by a factor of
three. Initial condition is a spike on 5 x 5 grid. Advected peaks show increasing
height with decreasing Ax.
79
difference cell at the new time level, could lead to improvements in accuracy
and efficiency.
80
2. Smoothed Aggregations Algebraic Multigrid
2.1 Introduction
In engineering practice, problems are often encountered in which the
material properties change by orders of magnitude from one part of the model
to the next. Such changes usually adversely affect the conditioning of the
linear algebraic systems obtained from discretization of these problems, and
pose a challenge to the methods used for solving these systems. Similar effects
may result from discretizations using locallyrefined meshes. Adding to the
difficulty, the problems are often discretized using unstructured meshes. This
renders classical multilevel methods impractical. Algebraic multilevel schemes
have been designed specifically to deal with the problems discretized over un
structured grids [18, 80, 86].
Efficient methods have been sought to handle algebraic systems aris
ing from discretization of second order partial differential equations having
subregions with markedly larger coefficients than the rest of the domain. Sev
eral nonmultigrid methods are studied in [7, 4, 5, 6]. We present a variant
of algebraic multigrid (AMG) designed for this purpose. This work is a joint
effort with Petr Vanek and Marian Brezina.
Aggregationbased multilevel methods [45, 44] using plain aggregation
involve low computational complexity. In this chapter, we present an algebraic
multilevel method of smoothed aggregation class. Such methods have previ
ously been proposed in [86, 91, 87, 19, 92], and have proved effective in solving
a variety of problems ranging from linear elasticity to the Helmholtz problem.
The classical aggregation method is a variational multigrid [9], i.e.,
the restriction operator is the transpose of the prolongation in which the pro
longation operator P/+1 from the coarse level / +1 to fine level l is constructed
based on the decomposition of finelevel nodes into disjoint sets, called aggre
gates. The aggregates then determine the nonzero structure of the prolongation
81
operator from the coarse level to the fine level. Each finelevel aggregate cor
responds to a single node on the coarse grid. As the aggregates are disjoint,
so is the nonzero structure of the columns of the prolongation operator corre
sponding to distinct aggregates. Given the finest level matrix Ax> the coarse
problems are constructed by the recurrence Ai+X = (P/+1)TA;P/+1.
For illustration, we consider the simplest example of the tentative
prolongator for the onedimensional Laplace equation discretized on a mesh
consisting of nx = 2>L~lni, nodes:
f 1
1
i
Pl
M+l
73
\
v
i
i
i
Each column corresponds to disaggregation of one IRn,+1 variable into three
IR71' variables, ni = 3n/+1. So, P/+1 can be thought of as a discrete piecewise
constant interpolation.
Since the matrix Ax A is tridiagonal, the aggregated coarselevel
matrices At, l = 2,..., L are tridiagonal as well.
The construction of aggregates in 2D or 3D [19] can be carried out
by a straightforward generalization of the above example. A typical rate of
nodal coarsening in aggregation will be by a factor of 3d, where d, denotes the
dimension of the problem.
Since the aggregates can be very inexpensively constructed by purely
algebraic means, based solely on the connectivity in finelevel matrix, aggre
gation suggests itself as a natural way to construct AMG methods. Note that
82
for scalar problems, the prolongation operator from coarse level of dimension
m to finelevel of dimension n stores exacly n nonzero entries. This sparsity of
the transfer operators, together with a fast rate of coarsening, leads to solvers
with very low algebraic and geometric complexities, as defined in [80]. These
methods, as well as improvements based on them, can be easily generalized
to solving both high order [87] and nonscalar [90] problems. This is achieved
by associating more than one column with each aggregate, and the resulting
method is referred to as the generalized aggregation method.
Unfortunately, the price for simplicity of aggregation methods is that
their multilevel convergence is suboptimal, as observed in practical computa
tion.
It is now understood that the reason for this suboptimality is that
the disaggregated coarselevel basis functions possess large energy. This un
derstanding has lead to development of the smoothed aggregation methods,
which address the above deficiency. Characteristic for the smoothed aggrega
tion methods is that their prolongation operators are constructed in two steps.
First, a socalled tentative prolongator P}+1 is constructed as in the general
ized aggregation methods. As noted above, this operator has a simple sparsity
structure. A smoothing operator Si is then applied to obtain the final prolon
gator Ill+l = SiPl+1 used in the method. The hierarchy of coarse problems is
then constructed by the recurrence
= (/;+1)tm+1. (2,i)
The prolongation smoothing improves energetic properties of the coarse space
basis functions and allows for approximation property (2.2) with constants
uniform with respect to the level.
Recently, multilevel convergence result has been presented by Vanek,
Brezina and Mandel in [89] for a smoothed aggregation method applied to
i71equivalent problems, including nonscalar equations of linear elasticity dis
cretized on unstructured meshes. The key assumption of that theory is satis
faction on all levels of the socalled weak approximation property, which can
83
