Citation
Eulerian-Lagrangian localized adjoint method and smoothed aggregations algebraic multigrid

Material Information

Title:
Eulerian-Lagrangian localized adjoint method and smoothed aggregations algebraic multigrid
Creator:
Heberton, Caroline
Publication Date:
Language:
English
Physical Description:
130 leaves : ; 28 cm

Subjects

Subjects / Keywords:
Eulerian graph theory ( lcsh )
Lagrange equations ( lcsh )
Multigrid methods (Numerical analysis) ( lcsh )
Eulerian graph theory ( fast )
Lagrange equations ( fast )
Multigrid methods (Numerical analysis) ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 119-130).
Statement of Responsibility:
by Caroline Heberton.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
45117599 ( OCLC )
ocm45117599
Classification:
LD1190.L622 2000d .H33 ( lcc )

Downloads

This item has the following downloads:


Full Text
EULERIAN-LAGRANGIAN 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
Jan Mandel
Thomas Manteuffel
^-dJo-QC)
Date


Heberton, Caroline (Ph.D., Applied Mathematics)
Eulerian-Lagrangian Localized Adjoint Method and Smoothed Aggrega-
tions Algebraic Multigrid
Thesis directed by Professor Thomas Russell
ABSTRACT
A three-dimensional implementation of an Eulerian-Lagrangian Lo-
calized Adjoint Method (ELLAM) for the solution of an advection-diffusion
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 one-dimensional 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
second-order 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


high-coefficient 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
IV


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. DMS-9706866 and Army Research Office Grant No. 37119-GS-AAS.
Chapter Two represents joint work with Petr Vanek and Marian Brezina.


CONTENTS
Figures ............................................................. ix
Tables.............................................................. xvi
Chapter
1. Eulerian-Lagrangian Localized Adjoint Method ...................... 1
1.1 Introduction..................................................... 1
1.2 Three-Dimensional ELLAM.......................................... 4
1.2.1 Governing Equation for Solute Transport........................ 5
1.3 Formulation of ELLAM 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 One-Dimensional Flow........................................... 26
1.5.2 Uniform, Three-Dimensional Flow ............................... 29
1.5.3 Two-Dimensional Radial Flow.................................... 33
1.5.4 Initial Condition in Uniform Flow.............................. 40
1.5.5 Constant Source in Nonuniform Flow............................. 42
1.6 One-Dimensional ELLAM........................................... 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 high-coefficient............................ 102
2.8 Multiple subdomains with high-coefficients........................ 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+1.............................................. 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 one-dim-
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 one-dim-
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 sec-1. 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 three-dimensional
solute transport in a uniform steady flow problem......... 35
1.12 Concentration contours in the horizontal plane containing the
solute source (layer 1) for three-dimensional 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 three-dimensional 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 three-dimensional 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 t = 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 two-dimensional 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 three-dimensional 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 three-dimensional 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
pre-images 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 with 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~1B)'S7)T with
T = 1,10,100,1000,10000 and Cr = §........................ 75
1.42 Initial condition with 2 nodes on a front, and results using 40x40
(S(A~1B)c7)t with T = 1,100 and Cr = |.................... 76
1.43 Initial condition with 3 nodes on a front, and results using 40x40
(S(A-1B)&)T with T = 1,100 and Cr = f..................... 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~1B)^)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 high-coefficient subregions touching at
a single node
118
xv


TABLES
Table
1.1 Parameters used in ELLAM simulation of solute transport in a
one-dimensional, steady-state flow system.................... 28
1.2 Parameters used in ELLAM simulation of transport from a con-
tinuous point source in a three-dimensional, uniform, steady-
state flow system............................................ 33
1.3 Parameters used in ELLAM simulation of....................... 37
1.4 Parameters used in ELLAM simulation of three-dimensional
transport from a point source with flow in the x-direction and
flow at 45
degrees to x- and y-axes................................... . 43
1.5 Parameters used for ELLAM simulation of transport in a vertical
plane from a continuous point source in a nonuniform, steady-
state, two-dimensional 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 lO*7, 10_t7................... 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 (1CT<7, 1017)..................... 118
xvii


1. ELL AM
1.1 Introduction
Solute transport in flowing ground-water has been intensively studied
in recent years in an effort to predict long-term 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 advection-diffusion 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 cost-effective for application to large systems over extensive
time periods, admitting an advective component to the equation may require
prohibitively fine discretizations. Consider a one-dimensional prototypical so-
lute transport equation,
dc dc d2c
dt ~>rV dx 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
Peg
vAx
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
L Pe nx ~i < Ax ~~ 2 (1.2)
A 2 L Ax~Te (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
, \ ( x vt\
c{x't)=eTflvm)
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 y/Pe
on a front to satisfy Ax proportional to A_. 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 dispersion-dominated, often varying in
different parts of the space-time 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 advection-diffusion transport equations. Eulerian
methods are characterized by the use of a fixed spatial grid. Petrov-Galerkin
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], transport-diffusion
method [79], characteristic-mixed finite-element method [2, 101, 3], operator
splitting method of [40, 100, 30], the Lagrangian-Galerkin method [76], and
Eulerian-Lagrangian 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 Eulerian-Lagrangian 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 constant-coefficient one-dimensional equations, then
extended to one-dimensional variable-coefficient 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 one-dimensional
finite-volume setting. One-dimensional nonlinear advection-diffusion equations
have been treated by Ewing [33] and Dahle, Ewing, and Russell [31]. Addi-
tional work has been done on ELLAM schemes for one-dimensional 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 one-dimensional case [82], Various
two-dimensional 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 three-dimensional ELLAM in a framework within which
all characteristic methods can be viewed.
In this thesis, a three-dimensional 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 one-dimensional context.
1.2 Three-Dimensional ELLAM
A three-dimensional finite-volume Eulerian-Lagrangian 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


MODFLOW-96, the USGS three-dimensional, finite-difference, ground-water
flow model [52], [54], [53]. The code is written in FORTRAN-77.
ELLAM [26] solves an integral form of the solute-transport equation.
The code uses an implicit method for dispersion calculations, which allows
for large time steps without stability constraints. The Eulerian-Lagrangian
approach involves tracking mass through time, then solving a dispersion equa-
tion on a fixed-in-space grid. This is particularly advantageous for advection-
dominated problems, as are typical of many field problems involving ground-
water contamination, since Eulerian-Lagrangian approaches tend to generate
less numerical dispersion than standard Eulerian finite-difference 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 (sD§i) E CW (1.4)
+X(eC + 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,
e is the effective porosity (dimensionless), V is a vector of interstitial fluid ve-
locity components (LT_1), D is a second-rank tensor of dispersion coefficients
(L2T-1), W is a volumetric fluid sink (W < 0) or fluid source (W > 0) rate
per unit volume of aquifer (T_1), C' is the volumetric concentration in the
sink/source fluid (ML~3), A is the reactive decay rate (T-1), 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 phase-equilibrium relationship in which C is proportional to C. The
retardation factor is defined as:
Rf 1 +
QbC
eC
An integral form of the solute-transport equation, which is a state-
ment of conservation of mass over the domain of integration, is the governing
equation for this finite-volume ELLAM approach. Equation (1.4) is integrated
against a space-time 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 space-time 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:
If
T V (eCV sBVC) - V C'W + u\eC) dtdx = 0,
otRf Rf J
(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
dt
~meC
and
Qj 1 1
V (eCV eBVC) = V (ueCV ueDVC) Vu (eCV eDVC)
Rf Rf Rf
6


to yield the global equation,
If
d^Cl + -Ly (ueCV ueDVC) + -EVu eDVC
dt Rf Rf
~TT £ C'w eC(lT + Vu Au)^ dtdx = 0.
Rf dt Rj J
(1.6)
The Eulerian-Lagrangian aspects of the method derive from the re-
quirement that the test function satisfy the adjoint equation,
Xu = 0 Thus, for the time step from tn to tn+1, we choose u of the form
u(x,t) = /(x, where ^ V/ = 0, so that / is constant
along characteristics of the retarded interstitial velocity field. Note that with
u = e-Mtn+1-t) (that is, / = 1) in the following, we arrive at a statement of
global conservation of mass over the time step:
- Ja(ueC)ndx
+ In +1 W} V {ueCV ueDVC) + ^ Vu eD VC
-jj-Z C'W dtdx
= 0.
To obtain a system of equations, each representing mass conservation
on a cell, we use u = X)z wj where l is the index for the finite-difference cells
Qi in the transport subdomain, and local space-time test functions are defined
with fi(x,tn+1) = 1 on Oi and //(x, tn+1) = 0 elsewhere:
ui
e_A(t"+1-t) on characteristics from any (x, t) G O x (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,
Jo, (ecr+1*c e"iA (£C) dx
+ JT \f{;v eDVC) n dtds
suppuiarn+1 *
- JT V, 2li^(eDVC) n dtds (1.8)
S ^ S CW dtdx
supp U[ ftsupp W 1
= 0
where 3- 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; Cl* is the pre-image at time tn under advection of cell Clt at
time tn+l; Tn+1 = <9f2 x (tn, tn+1) is the space-time 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 space-time 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 right-hand side
of the system of equations.
1.3 Formulation of ELLAM equations
ELLAM approaches the hyperbolic-parabolic 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 finite-difference
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. User-input inflow con-
centrations, together with the outflow solutions, then appear on the right-hand
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-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,
/n,(sC)+1 dx Atj9ni i(eDVC)+1 -11* = e~^ /n;(eCy dx
S e-si-tn*'-t)eCiHi,cw^r ndtds
suppuiC\Tn+l
+ If
supp U(flsupp W f
where f2* means the pre-image in the spatial domain at tn of at tn+1, Tn+1 =
dQ x (tn,tn+1) is the space-time 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.
(1.9)
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:
fQ j£+1 V (usCV ueDVC) dtdx =
C1 fdn jqiueCV ueDVC) n dsdt.
Considering just the outflow portion of the boundary, this becomes,
j£+1 w- (ueCV ueDVC) n dsdt
J(on)outfiow Rj \ / (1.10)
f(dn)out/iow hn (VZCoutflow-^) n dtds,
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:
' e_A(r+i_t)
uu =
0
on characteristics from fi at time level n
into boundary area (dd)u at any time during time step
otherwise.
(l.n)
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,
S{dn)u J?+1 e X{tn+1 t]£Coutfiowij ndtds = (1.12)
10


123456789 10
source
Figure 1.2. Outflow boundary receives advected mass from inflow, source in
cell 1, and storage in all cells.
e~XAt f(dn)h(£C}n d*
If e-W^-QeCmfiovjt ndtds
suppuuC\do.infiow X (tn ,tn+1)
+ if e-xr*'->zc%dtdx,
suppuur\suppWx(tn,tn+1) 1
where (dQ,)u is a boundary face, and (dQ)^ is the preimage at time tn of
(d£l)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 steady-state flow condition. The specific discharge, or flow
11


rate per unit area, across each face of every finite-difference 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 finite-difference grid, the integrals on the
right-hand 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 piecewise-constant 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_Vz_
dt Rf
(1.15)
This is accomplished by introducing a set of moving points that can be
traced within the stationary coordinates of a finite-difference 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 pre-image 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 right-hand 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 finite-difference grid in-
dexing, with this index order indicating an x,y,z sequencing. A right-handed
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+1.
14


are calculated using centered differences as in M0C3D, generalized for varying
grid dimensions. (See [71], page 64, for the form of the ebDm §£- expansions.)
OXm
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 lAr A'ihhn+1 Y'8 (1
porosity 64-^corner=l ^corner fl 17l
= ^corner=l ^2nbr=l(We^9^)nbrCnbr,
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 right-hand 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 Xj^i + Ax.
-) + (
Axi
Axj
j+1
+ Axj
)-4)
15


x (() + () 4)
{{+ Ay/^ {Ayi+1 + Ayi} j
x r-Tk) + (i
_ ^
)(

AXji -|- AXj Ayj_i I- Ayj bj,i,k 4" bj,i,k1
b

_l_ / \ /- \ /
Axj-i + Axj A^_i + Ayj bj^ + fy.i,,
fc+i
)c:
j-l,i-l,A;+l
+ (
Arc,-
3 -Ktt.--1? x..")(i bjAA-
A yi
Aa;j_i 4- AxjJ v Ai/j+i 4- AyjJ v + bj^k+i J +
+ (

Arrj_i + AXj 'v Ayi+i + Ayj' v bjiiik + bj^k-1
^ ^_________________V_____ bjtj,k
)Cj
^Vi1 ^Vi AXj+i + bj}i,k I
bj,i,k
v Jhhk \n
JU i------------JW+l,i-l,fc-i

\/ bjj tk ^
'Al------~--------J^+i,i-i,fc+i
^Vi1 I ^2/z ^*^3+1 t- &Xj I ^j,z,A:+l
-)(-
A yi
^ ^ bjtitk
+ (------AXj
Azj+i 4~ A Xj A^/j_j_i 4- A yi bjjtk 4" bjtiki
_l_ / ^xj w ^Vi w bjtitk
A.Xj+1 4- A Xj A 2/j+i 4- Ayj bj^^k 4- bjtij~+1
,, Aa;,- , , Az,- .
+ (fc---------4rr) + fc----------4rr)-4)
)Q+l,i+l,fc+l)
3 Axj_i + A Xj Azj+i 4- Axj
X ((
+ (
A Vi bjtitk
Ayi_i 4" A yi bjtitk 4- bj^^ki
)Cj,i-l,k-l
^Vi \f bjtiik ^
Hi z. )'^j,i-l,k+1
^Vii 4~ Ayj bjtitk 4" bj^^+i
, ( AVi \f bj^k ^
+ r^rrAi n--------Jw,i+i,fc-i
+ (
^2/i+l 1 ^Vi bjtitk 4- bjtitkl
^Ui ^ ^ bjtj}k___________________
^Vi+i H Ayj bjtitk 4- bj^^+i
)Cj ,i+l,A:+l)

X ((
' Ayj_i 4- Ayj Ayi+l 4- Ay,
A Xj s
)(l
A:Cj_i 4- Asj bjtitk 4-
)Q-i4,
i,A:+l
, ASj 6j,j,fc
Axj+1 + 4- b^kJ^1^-1
16


+ (
+ (
^-------)(l----------)Ci+iA
/S.Xji -f- AXj bj,i,k 4 bj>i,k+l
k+1
Ax,-
)(l
bi.i
j,i,k
A.Xj1 4 AXj bj,i,k 4

+ ((
bj,i,k -) + (-
bj,i,k + bj,i,kl bjjfi 4 1
)-4)

' A yi_! + A yi Axj+i + Ax

' Axj-i 4 Axj A^_i + A y{
il,k

- --------------------)((
Axj
bj,i,k 4" i A Xji 4- A Xj
x ((-----------------) + (-----------------) 4)
{{Ay^ + AVi] ^ [ Ayi+1 + AVi} }
AXi ) 4)
+ (*r^hrr)((
Aj/j_i + Ain Ax j+i + Ax
+ (
A Xj
)((-.-. -VA-.) ~4)
Axj+i + A Xj Ayi+i + Ay,
- 4((
- (
Aj/i
Aj/i+i 4- Ay,
bj,i,k

)((-
Azj
bj,i,k 4 1 Axji -I- Ax
' Ayj_! + Ay, Ayi+1 + A yt J
+ (A.. Ay\j((
Ax,-
A^-x 4- Ayi Axj+1 + A xj
)-4)
, Ax,-
+ (- J
)T^C7)-4)
Axy+1 + A Xj Ayi+1 + Ay,
- 4((
- ((
A yi
A yi+i 4 A y{
) ~ 4))Cj>i)fc+i
bj,i,k
-) + (i
b<
j,i,k
bj,i,k 4 bj i^fci bjtitk 4 bjtitk+l
)-4)
17


X(((
+ (
X ((
Axj
Axj-i + Axj
Ay,
A yi+i + A yi
Axj
Axj-i + Ax,
) + (-
Axj
A Xj+i + Axj
XW)(( Ar
)-4)((
A yi-i + Ay,
) + (
A yi
Ayi_i + A yi
A yi
Axj
Axj+i + Ax,
A yi+i + Ay,,
)Cj+l,i,k))-
)"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 right-hand 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 fli.
This yields the formulation,
AAt
/n* &
AxjAyib
Hy; = p-AAt V V
k k (NSC)(NSR)(NSL)
subcell
center
Wi
Cp/)C(p)
(1.18)
where summation runs through all subcells of each cell in the transport sub-
domain, and p-f 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


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 one-dimensional approximate test function in each direction:
wjik(x,y,z) = f(x)g{y)h{z).
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
X Xj
X = Axj
and similarly for y and z. For an interior cell on a uniform grid with all
surrounding cells active,
/(£) =
/
0
NSCx + \{NSC + 1)
< 1
-NSCx + l(NSC + 1)
0
\
2 2 NSC
_1___1_ <£<_!+ 1
2 2NSC ^ x ^ 2 2NSC
_I -|_1 < x < I_L_
2 ^ 2NSC x 2 2NSC
1 __1_ 2 2NSC ^ ^ ^ 2 T 2NSC
2 ^ 2NSC x
and similarly for g and h. These functions are shown graphically in figures 1.4
and 1.5.
19


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
test functions corresponding to / defined above are illustrated graphically in
figure 1.6.
For a point in fi/ with reference coordinates x,y,z G [ \,\], f is
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=)) giy) h(z),
Wj,i+i,k = f(x) (1 g(y)) h{z),
Wji,k+i = f(x) g(y) (1 h(z)),
Wj+i,i+i,k = (1 fix)) (1 g{y)) h{z),
Wjj+iMi = fi*) i1 ~ 9{y)) (1 /i(2)),
wj+i,i,k+i = (1 f{x)) giy) (1 h{z)),
wj+1>i+1>k+i = (1 fix)) (1 giy)) (1 h(z)).
Figure 1.6. Approximate test functions in one direction on a nonuniform grid,
with NS=2.


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 [5,5], 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
sub-time 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 sub-time 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 right-hand side of the local
conservation of mass equation for cell 11; with the following integration. The
21


domain of integration is all sources which intersect the space-time test function
for cell Q.I, so that all source mass that flows into Qi during a time step is
integrated. Multiple sources within a source cell are summed.
If
supp uiHsupp W
dfdx

NT
E EE1
all P= m= 0
sources source
subcell
center

A T
NSC NSR NSL
Wi
(p) E c.%- (1'20)
Kf
source
in
cell
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 = Al,
or if m = 0 or NT; and tm = tn + At 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 0/, 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 right-hand sides of the equation
for sink cell Ot with coordinates j, i, k. This gives
22


SUppUlPiSUpp W
If
e_A(t+i_t)i ^C'f- dtdx
itf
Atfnsinlc e {Qa,ve E + Cetjff) dx
( \
\
(1.21)
At (C"+1+Cn) C Qet V
ej.i.k 2 ^ Eu R >
5= *7
sink
\ in I
\ cell /
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
composite trapezoidal rule is applied in time. At each sub-time 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 two-dimensional boundary face is discretized,
while for a source, the entire cell is discretized. For a cell Qi, the integration
is performed over the intersection of the space-time test function for that cell
and the transport subdomain boundary; that is, all mass entering through the
boundary and advected to Q,i during the time step is accumulated to the right-
hand side of local equation l. Mass advected into S7; during a time step from
the inflow portion of the transport subdomain boundary is give by,
discretized into a number of sub-time steps determined by parameter NT. The
ff e~*(tn+1-QC
suppuill rn+1
NT
E E E e-A(+'")
ref area
wi(pf)Cinfi0WQ (1-22)
all P= 771=1
inflow fo.ce
faces subarea
center
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 + represents the time during time step at which discretized
source mass enters the system; and AT = or 7^ if to = 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 fi;.
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,
f(dn)u f£+1 e_A(f"+1 _t) Coutflow jij n dtds
~ ^outflow 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 right-hand 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 right-hand side of (1.12) is accomplished using a
trapezoidal in time, midpoint on subcell in space approximation:
24


ff e~^tn+1^C^--ndtds
suppuucrn+1 *
0XAt
?fe 5 (NSC){NSR)(NSL)Wu^C^
5 ... L_11
NT
+ E E E
all P= 771=0
inflow fo-ce
faces subarea
center
e-M*n+1-tm)
AT
ref area
IVll ip^ ) in flow Q
(1.24)
NT
ej,i,k
E EE1
p= m=0

A T
all
sources source
3ubcell
center
NSC NSR NSL
Mp) E
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 fi; has been replaced by the approximate test function for boundary
face U. 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 half-life 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 ground-water 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 three-dimensional ELLAM code was tested by running the same
suite of test cases as was applied to the USGS three-dimensional 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 One-Dimensional Flow
The first test case is a simple one-dimensional system involving solute
transport in a finite-length aquifer with flow of constant velocity. Boundary
conditons are third-type, 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 -* yy 0.01 cm2/s
£ 0.1
OiL 0.1 cm
CXj'H = OiJ'V 0.1 cm
PERLEN (length of stress period) 120 sec
vx 0.1 cm/s
Vy = Va 0.0 cm/s
Initial concentration (C0) 0.0
Source concentration (C1) 1.0
Number of rows 1
Number of columns 122
Number of layers 1
DELR (Aar) 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
one-dimensional, steady-state flow system.
28


Figure 1.7. Plots of concentration as a function of cell node for one-dimen-
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 one-dimen-
sional flow with a constant velocity field and high dispersion. Shown are the
analytical solution (lowest graph), ELLAM 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 method-of-characteristics 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, Three-Dimensional Flow
Next, ELL AM results were compared with the analytical solution
developed by Wexler [99] for three-dimensional 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
cross-terms 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


1
0.4 -
0.2 -
01-------1------1------1-------1-------1------1-------1-------'------1-------1-------1------1
10 20 30 40 50 60 70 80 90 100 110 120
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 sec-1. 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 J-xx ^yy 0.0125 m2/day
£ 0.25
OIL 0.6 m
aTn 0.03 m
Oijpv 0.006 m
PERLEN (length of stress period) 400 days
vx 0.1 m/day
K = VZ 0.0 m/day
Initial concentration (Co) 0.0
Source concentration (C') 2.5106 g/m3
Q (at well) 1.010-6 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 three-dimensional, uniform, steady-state 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 Two-Dimensional Radial Flow
The ELLAM solution was compared to the analytical solution give
by Hsieh [61] to a radial flow and dispersion problem with a finite-radius well
in an infinite aquifer of two-dimensions. 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 three-dimensional 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 three-dimensional solute transport in a uniform steady
flow problem with CELDIS = 7 (two time steps), NSC = NSR = NSL = 4,
NT = 16.
35


Figure 1.13. Concentration contours in the horizontal plane containing the so-
lute source (layer 1) for three-dimensional 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 three-dimensional 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 -i-yy 3.6 m2/hour
£ 0.2
OIL 10.0 m
dj'H = Oij'V 10.0 m
PERLEN (length of stress period) 1000 hours
Q (at well) 56.25 m3/hour
Source concentration ((7) 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 ELLAM simulation of
two-dimensional, steady-state, 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
non-physical 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 three-dimensional 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 finite-difference 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


Figure 1.17. Contour plot of concentration for run with 29 time steps using
CELDIS = 5, NSC = NSR = 4, NSL = 2, NT = 4.
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 grid-orientaion effect caused primarily by off-diagonal terms
of the dispersion tensor.
42


Parameter Value
T T -* xx yy 10.0 m2/day
£ 0.1
aL 1.0 m
CXj'H OLj'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 (As) 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 three-dimensional
transport from a point source with flow in the x-direction and flow at 45
degrees to x- and y-axes.
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 r~
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
-|------------1------------1--------------1--------------1------------1--------------r~
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 two-dimensional 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 cell-centered finite-difference 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 two-dimensional 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
OLL 3.0 m
Oij OL'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.2222-0.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, steady-state, two-
dimensional flow system (described by Burnett and Frind, 1987).
48


Figure 1.27. Contour plot showing concentrations of ELLAM 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


20
10
0
r uu (a) 3D Finite-Element Model
- 1 1 L Contours: 0.1 to 0.9
0 40 80 120 160 200
Figure 1.28. Contour plot showing concentrations of a three-dimensional
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
&tv 0.01m and ath = 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 cxtv = 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 One-Dimensional ELLAM
Consider a one-dimensional problem with e = Rj = 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 a;<
< 1 < x < 1 ,
. 0 k 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
o
Figure 1.30. Contour plot showing concentrations of a three-dimensional
finite element solution of the Burnett and Frind problem with high vertical
transverse dispersivity. Contours shown are 0.1 to 0.9.
(a) 3D Finite-Element Model ,
200
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+l = 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] (1-26)
denoting the sub-, main, and super-diagonal 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


i-1 i i+1
x-I/2 x+1/2
Figure 1.32. Mass at old time level is integrated exactly using cell centers
and pre-images of cell boundaries as integration points, and applying the trape-
zoidal integration rule.
54


tn+1, and interpolating concentration values between nodes, the trapezoid on
the left has area,
and the area on the right is,
- CrAx
Ci + C(
Since
c_i) = (i + cr) a-1 + (5 o) a
and
C+i) = (5 + Cr) ft + (5 Cr) Ci+1,
the mass in cell i at the new time level is
Ax
((2 + ^r)(l + Ci-1
+ ((2 + CV)(| T1) + (2 ~ Cr)(l + T"))
+ ((5 Cr)(| f) ci+1;.
Thus, disregarding boundary conditions, the B matrix is tridiagonal with
1 Cr Cr2
S = -[1,6,1] + [1,0,-1] + -2-[i, -2,1],
(1.28)
denoting the sub-, main, and super-diagonal 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 diagonal-shifted 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 three-dimensional ELL AM implementation. Application of
the three-dimensional code to a one-dimensional 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+l) = Bc{-, tn) + A iec, (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 is nonnegative definate, this implies
en+1 (A + M)~1Ben + A tec>
with a different constant incorporated into ec. This leads to
en+1 = f^((A + M)~lB)jAtec + ((A + M^B^+'e0,
i=o
where e denotes any initial error. For stability we then require,
\\{(A + M)~lB)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(iL)||ec||, and consistency yields
convergence as Ax, At > 0.
Stability depends on the properties of the matrix (^4 + 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 = [\-k^ + 2k,\-k], (1.30)
56


(1.31)
The right-hand side matrix B becomes
B = |[1,6,1],
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,
[ + f(x) dx
= /^i/a f(xi) + f'(xi)(x ~ xi) + \f"{xi){x ~ xif + ((x ~ xif dx
= Axf(xi) + 0 + yfa) /**+;y22 (x Xi)2 dx + 0
= Axf(xi) + ^f-f"(xi) + 0{ Ax3),
where it is useful to separate the and 0(Ax3) terms. Then, we have
l/n+1 Jxf'f ct(x,t) dxdt
= Ax fn+1 Ct(xi, t) + ^~-g^ct(xi, t) + 0(Ax2) dt
= Ax C1 ct(xi} t) + ^f i^lct(x^)+ct(xi+1,t),f) + 0(Aa.2^ + 0(Ax2) dt
= Ax (crx c?+1 [cr1 (?_! 2(Cr +0(A£Ax4) + 0(AtAx2))
= Aa; - i~ i)+K^r1 cD+i(Cr ci+1)+o{AtAx2)).
Also,
rn+1
Ca,2,(x,£ ) I ^xxt(^j &) ds
rtn+l
~ ^'xxis^i}^' ) "t" Cxxxi^iit ) (x X^) T O(Ax ) J ^xxt(^j &) ds
r"+l o^-K _i_
= (A^ +c:Ea;a;(^,r+1)(a:-a:i) + 0(A3;2)
- / cMt(a:, s) ds.
57


Integrating then gives
L
i+l/2
cxx(x,t) dx
i1/2
= Ax
n+l _9_n+l i _n+l
Ci-1 ZCi +ci-
(Ax)2
rtn+1
+ 0(Ax3)
ft cXxt(x, s) dscfo;
and
>xi-l/2
,,72+1 , n_Jl+l
[xi+1/2 rtn+1
/ / Dcxx(x,t)dtdx
JXi_ ,/> Jtn
= Aa;^(-C"-11 + 2cil+1 cm) + /t-+I /x3;^2 Jf"* cXIt(x, 5) dsdxdt
= Axf^(-C+11 + 2c+1 c^1) + 0(Ax3At) + 0(AxAt2).
Combining terms gives the estimate,
rtn+1
/t"+1 fltl',1 C*(X dxdt ~ f
c-1/2
tn+1
tn
fx£$ Dcxx(x, t) dxdi
= Ax (Kc?/-,1 Ax-gjl(-c+1 + 2c+1 £"/,') + 0(Ax3A() + 0(AxAt2).
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
^ + M = / + (t-l)[-l,2,-l]
and (1.31) as
B = /-i[-l,2,-l]
where I is the identity matrix. The matrix (A + M)~lB has eigenvalues
1 |A
l + (*-!)A
58


where A is an eigenvalue of [1,2, 1], Since 0 < A < 4, (A + M)~1B has
nonnegative eigenvalues less than or equal to one for any k > 0. Because
(A + M)~lB is normal, this gives \\(A + M)~lB\\2 < 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 A; > |, A~lB satisfies a discrete
maximum principle. For 0 < k < |, however, the operator need not satisfy a
maximum principle; consider, for example, 3x3 A~lB with k = T. Here,
A B
4
679
/ 2721
14
-14 1 \
-14 196 -14
1 -14 195
679
V
/ 16130
112
91
4
1559
112
574
4
56
4
/
-l\
__^ 5b 578
56
4
578
4
/ 6 I 0 \
8 8 u
8 8
\ U 8 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((Ax)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)_1B = A~lB
= / + .4-1(f[l,0,-l] + ^[l,-2,l])
= / + f A-' ([1,0, -1] + Cr[ 1, -2,1]).
59


Then, since 0 < Cr < ~,
|\(A + M)-lB\\2 =\\a~1b\\2
-l]||2 + Cr||[l>-2,l]||2)
< 1 + f |]A-1||2(2 + 4C'r)
< l + 2Cr||A-1||2
= 1 + O(At) for fixed Ax.
We have
||(A-1B)^||2 <\\A-'B\\?
< (1 + 0(At))S
< 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 left-hand side (1.26) and right-hand 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~lB is a Laurent matrix
with norm of one, and note properites of its spectrum.
Given a sequence {^-n}^T_oo of complex numbers, a doubly-infinite
60


matrix with 2_n along the nth diagonal is called a Laurent matrix:
... z0 ... Zi 21 zo Z-2 2-1 Z-3 2-2 2_4 ... 2-3 ...
... ^2 Zi zo Z-1 2-2 ...
... 23 Z2 Zl Zo 2_i ...
... 24 z3 Z2 Z1 Zo ... J
The following theorem is well-known:
(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 {^n}^L-oo is the sequence of Fourier coefficients
of z:
1 f2?r
2n = -p= / z{eme)e-medd (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 £(z). Let L := L(T) and L2 :=
L2(T).
Laurent matrices represent multiplication operators on L2 with re-
spect to the orthonormal basis
lx/2T
AnO
n£Z
Define
M(z) : L2 -> L2, f^zf.
This multiplication operator is bounded if and only if z G L^, in which case,
11X0011 = Nice,
(1.33)
61


Consider the operator which maps a function to the sequence of its
Fourier coefficients,
$ : L2 l2(Z), / i ^ {fn}nez-
The operator $ is bijective and
WII2HI/II2
for / G L2. Then
C(z) = $M{z)$-\
(1.34)
(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~XB = C(a~lb),
(1.36)
with
a~lb ------ ( (cr ^ cos(0) + Cr elB + 7 Cr2
3 + cos(0) \\ 2J y J 4
and ||£(a_1J>)|| = 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
(*) =i75T(e-" + 6 + e")
= + COS(0)).
The continuous function b G L(T) with Fourier coefficients composing the
diagonals of (1.28) is given by
K^) jkz
1
1
V27
(I f + ¥) + (I CC + (1 + f + SI) e
(Cr cos(0) + Cr e10 + | Cr2
(1.38)
62


We thus have £(a) and £(&). Furthermore, since zero is not in the the inverse of £(a) is £(a_1), (see [14], Theorem 1.2), where range of a,
_x 4x/2tt 3 + cos(0) (1.39)
From (1.35), it follows that
C(zi)C(z2) = £(2422) for all Zi,Z2 G L. (1.40)
So £(a l)C{b) = C(a 1b). Multiplying continuous functions (1.39) and (1.38), we get the bounded linear Laurent operator on an infinite domain,
A-XB = £(a_1&),
with a~lb given by (1.37). From (1.35), (1.34), and (1.33), it follows that
||£(*)|| = IMloo for all z G L. (1.41)
Thus, to determine the norm of £(a 1b), we find the maximum mod-
ulus of the continuous function a-1 b. For
p = (4Cr2-l)2,
q = 2(4Cr2 + 1)(3 4Cr2),
r = 16CV4 8CV2 + 9,
we have
I ui2 Pcs2(0) + qcos(6) + r
|a 1 (3 + cos(0))2
Then
^ =frisfe[2'3, + (5 + 6P)coS(9)l
= (34C'r2(2Cr2 1)) (coS(0) 1)
= 0
for 0 such that sin(0) = 0 or cos(0) = 1. Thus the critical points for |a_16|2 are
0 and 7r when 0 G [0, 27t). By (1.37), for 6 = 0,
a~lb = 1,
63


and for 9 = tt,
a~lb = 1 4CV2 G [0,1].
Thus the maximum modulus occurs at 9 = 0 and is equal to 1. Then (1.41)
gives ||£(a_1&)|| = 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 :={AeC:z Ae is not invertible in the algebra},
where e is the algebra identity element.
For z 6 L, further define the essential range of 2 as the set
TZ(z) := {A e 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 £ 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 £(a_16) 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_15| = 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,
£(a_16) has no eigenvalues. The entire spectrum is seen to be continuous.
Also, 0 G sp(a-15) only for a Courant number of
Lemma 1.8 For a~lb given by (1.37) and Cr G [0, |],
0 G sp(o_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 eie Cr2 + §
= (cos(0) 1 )Cr2 + (el9 cos{9))Cr + \ cos(9) + |.
For 9 ^ 0, using the quadratic theorem,
__ cos(0) e,9y/e2lt> 2cos(d)es +cos'2(0) (cos($) l)(cos(d)+3)
L,r = 2(cos(0)-l)
cos (8)ee 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 9 G [0, 27t), this means 9 = 0 or 9 = tt.
If 9 = 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_1H, 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 \\A-^B\\ tridiagonal ^1,1 An)n U-'bu tridiagonal \\A~lB\\a Ai}i A.nin =
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
g{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 = ABx, 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
\\A-xB\\2 = q{(A~1B)t A~lB)
and
p-^IU = e(A*BTA-1BA-2) = q((BA~1)tA~1B),
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,
Cr\Axf
c, + AALClIt CrACt
2A t
-a
o,
where the coefficients are taken to be constant. With the opposite sign on the
space-time 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 well-known (cf. [83]) that ELLAM methods, in general, require
sufficient mesh-density for about four grid nodes on a front for accurate model-
ing of advected concentration. (Compare with Eulerian methods where number
67


1
_Q5I___I___I____I___I____I___I____I____I_____I___I__
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 no-diffusion operator
was studied. Qualitative features of the solution are independent of the problem
size. The following represent results for the one-dimensional 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 E (|, |). 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


1
0.8
0.6
0.4
0
-0.2
_J________I________l_
20 40 60
80
100
_l________I________I
120 140 160
200
Figure 1.34. Initial and advected peaks with Cr
Figure 1.35. Initial and advected peaks with Cr =
69


1
0.2
_0 2^_______I_______1______1_______1_______1_______I_______1_______1_______1_______1_
20 40 60 80 100 120 140 160 180 200
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~1Bin cases where Cr = intlger Here, S is the diagonal-
shift by 1 operator,
5 = [1,0,0],
Disregarding boundary effects, application of 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 (S'(A-1H)^) of various sizes. (This is also
seen for ((A~lB)^S, which first shifts, then advects; and also (S(i)(A~1B)^)
and ((A~1B)'&S(i)), with advection and shifting of i cells.) The eigenvector
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 ill-conditioned. 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 5(A 1B)oT 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 (i>(A lB)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 ill-conditioning 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 ill-conditioned 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 well-discretized 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) cT )t with
T = 1,10,100,1000,10000 and Cr =
75


1.2
)
35
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-1.B)£)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(A~1B)'^)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 advection-diffusion equa-
tion can simulate the transient, three-dimensional 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
three-dimensional 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 non-physical 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 one-dimensional 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 Ax.
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 locally-refined 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 non-multigrid 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.
Aggregation-based 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 l + 1 to fine level l is constructed
based on the decomposition of fine-level 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 fine-level 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 A\, the coarse
problems are constructed by the recurrence Ai+1 = (P/+1)TA;P/+1.
For illustration, we consider the simplest example of the tentative
prolongator for the one-dimensional Laplace equation discretized on a mesh
consisting of n\ = 3L-1n£, nodes:
(l
1
i
\
1
1
1
V
1
1
1

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 coarse-level
matrices Ai, 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 fine-level 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 fine-level 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 coarse-level 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 so-called 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+1 = SiP[+1 used in the method. The hierarchy of coarse problems is
then constructed by the recurrence
= (IW)TM+1. (2.1)
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
/^-equivalent 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 so-called weak approximation property, which can
83


Full Text

PAGE 1

EULERIAN-LAGRANGIAN 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

PAGE 2

This thesis for the Doctor of Philosophy degree by Caroline Heberton has been approved by Andrew Knyazev Jan Mandel Lj.-Jb-OO Date

PAGE 3

Heberton, Caroline (Ph.D., Applied Mathematics) Eulerian-Lagrangian Localized Adjoint Method and Smoothed Aggrega tions Algebraic Multigrid Thesis directed by Professor Thomas Russell ABSTRACT A three-dimensional implementation of an Eulerian-Lagrangian Lo-calized Adjoint Method (ELLAM) for the solution of an advection-diffusion 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 one-dimensional problem is developed. Stability is demonstrated in the purely diffusive and purely advective limits, under assumptions on boundary conditions, and in the case of an infinite 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 second-order scalar problem with jumps in coefficients is presented. An algorithm is developed with components shown to satisfy abstract convergence criteria. The algorithm uses a coarsening strategy which respects boundaries of iii

PAGE 4

high-coefficient 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 candidate's thesis. I recommend its publication. Signed __ IV

PAGE 5

ACKNOWLEDG MENTS Portions of Chapter One of this thesis are being published as a part of a U.S. Geological Survey Water Resources Report by C. 1. 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. DMS-9706866 and Army Research Office Grant No. 37119-GS-AAS Chapter Two represents joint work with Petr Vanek and Marian Brezina.

PAGE 6

CONTENTS Figures Tables Chapter 1. Eulerian-Lagrangian Localized Adjoint Method 1.1 Introduction......... 1.2 Three-Dimensional ELLAM 1.2.1 Governing Equation for Solute Transport. 1.3 Formulation of ELLAM equations 1.3.1 Cell Integral Equations .... 1.3 2 Outflow Boundary Equations 1.3.3 Assumptions ... 1.4 Numerical Methods. 1.4.1 Mass Tracking .. 1.4.2 Numerical Integration 1.4.3 Dispersion . 1.4.4 Storage at New Time Level 1.4.5 Mass Storage at Old Time Level 1.4.6 Approximate Test Functions 1.4.7 Source Integral 1.4.8 Sink Integral VI IX xvi 1 1 4 5 8 9 9 11 11 12 14 14 14 17 18 21 22

PAGE 7

1.4.9 Inflow Boundary Integral 1.4.10 Outflow Integrals 1.4.11 Decay .... 1.4.12 Assumptions 1.5 Test Problems 1.5.1 One-Dimensional Flow 1.5.2 Uniform, Three-Dimensional Flow 1.5 3 Two-Dimensional Radial Flow ... 1.5.4 Initial Condition in Uniform Flow. 1.5.5 Constant Source in Nonuniform Flow. 1. 6 One-Dimensional ELLAM . 1.6.1 1.6.2 1.6.3 1.6.4 1.6.5 1.7 Limiting Case of No Advection Limiting Case of No Diffusion Infinite Spatial Domain ... Further Considerations for Finite Spatial Domain Numerical Dispersion and Oscillations Conclusion .... 22 23 25 25 26 26 29 33 40 42 51 54 57 58 63 65 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 Vll

PAGE 8

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 sub domain with high-coefficient 102 2 8 Multiple subdomains with high-coefficients 106 2 9 Computational Experiments 112 2.10 Conclusion. 114 References 116 VIll

PAGE 9

FIGURES Figure 1.1 Cell 7 receives advected mass from inflow, source in cell I, and storage in cell 1.. . . . . . . .. 10 1 2 Outflow boundary receives advected mass from inflow, source in cell I, and storage in all cells. . . . . .. 11 1.3 Preimage of cell may be irregularly shaped and not easily d e limited by backtracking. Instead, the known mass distribution at time tn is tracked forward along streamlines of the advective flow to time tn+1. . . . . . . .. 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

PAGE 10

1.7 Plots of concentration as a function of cell node for one-dim 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 one-dim ensional flow with a constant velocity field and high dispersion. Shown are the analytical solution (lowest graph), ELLAM 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 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 constant A = 0.01 sec-i. 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

PAGE 11

1.11 Concentration contours of analytical solution in the horizontal plane containing the solute source (layer 1) for three-dimensional solute transport in a uniform steady flow problem. . .. 35 1.12 Concentration contours in the horizontal plane containing the solute source (layer 1) for three-dimensional 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 three-dimensional 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 three-dimensional 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. Concentration maximum is 1.019. . . . . .. 40 Xl

PAGE 12

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 1.20 Contour plot showing log of concentrations in analytical soluion of Dirac problem at t = 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 two-dimensional finite element solution of the Burnett and Frind problem. Contours shown are 0.1 to 0.9. . . . . . . 47 Xll

PAGE 13

1.27 Contour plot showing concentrations of ELLAM 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 three-dimensional finite element solution of the Burnett and Frind problem. Contours 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 three-dimensional fi-nite element solution of the Burnett and Frind problem with high vertical transverse dispersivity. Contours shown are 0.1 to 0.9 .................. 1.31 Contour plot showing concentrations of ELLAM solution to high dispersion Burnett and Frind problem using CELDIS = 21, NSC 51 = 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 pre-images of cell boundaries as integration points, and applying the trapezoidal integration rule. . 1.33 Initial and advected peaks with Cr = 10. 1.34 Initial and advected peaks with Cr = 1.35 Initial and advected peaks with Cr = %. 1.36 Initial and advected peaks with Cr = Xlll 54 68 69 69 70

PAGE 14

1.37 Initial and advected peaks with Cr = 10 3 12' 1.38 Eigenvectors of 10 x 10 S(A-1B)dr with Cr = 1.39 Initial condition and results using 10 x 10 (S(A -1 B) dr Y with 70 72 T = 1,10,100,1000 and Cr = . . . . .. 73 1.40 Initial condition and results using 40 x 40 (S(A-1 B) dr Y with T = 1,10,100,1000,10000 and Cr . . . .. 75 1.41 Initial condition and results using 200 x 200 (S(A-l B) dr)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(A-l B) dr Y with T = 1,100 and Cr . . .. 76 1.43 Initial condition with 3 nodes on a front, and results using 40x40 (S(A-IB)dr)T 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-IB)dr)T with T = 1,100 and Cr . . .. 77 1.45 Initial condition with 5 nodes on a front, and results using 40 x 40 (S(A-l B) dr Y with T = 1,100 and Cr . . .. 77 1.46 Superposition of results showing initial condition and advected peak using Cr = Grids are 5 x 5, 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 !:J.x. .. 79 1.47 Superposition of results showing initial condition and advected peak using Cr = Grids are 5 x 5, 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 !:J.x. .. 79 XIV

PAGE 15

2.1 Configuration with two high-coefficient subregions touching at a single node. . . . . . . . 118 xv

PAGE 16

TABLES Table 1.1 Parameters used in ELLAM simulation of solute transport in a one-dimensional, steady-state flow system. . . .. 28 1.2 Parameters used in ELLAM simulation of transport from a con tinuous point source in a three-dimensional, uniform, steady-state flow system. .......... ... 33 1.3 Parameters used in ELLAM simulation of 37 1.4 Parameters used in ELLAM simulation of three-dimensional transport from a point source with flow in the x-direction and flow at 45 degrees to xand y-axes 1.5 Parameters used for ELLAM simulation of transport in a vertical plane from a continuous point source in a nonuniform, steady state, two-dimensional flow system (described by Burnett and 43 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 eigenvector decomposition of spike initial condition at node seven on 10 x 10 grid . . . . . . . .. 72 xvi

PAGE 17

2.1 Checkerboard pattern, coefficients 10<7, lO-CT. . . .. 116 2.2 Two cubes touching at a node Coefficient=10r in dark subregions (see Figure 2.1). ....... . 117 2 3 Element coefficients random in (lO-CT, 10CT). 118 xvii

PAGE 18

1. ELLAM 1.1 Introduction Solute transport in flowing ground-water has been intensively studied in recent years in an effort to predict long-term 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. S/olution of an advection-diffusion 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 cost-effective for application to large systems over extensive time periods, admitting an advective component to the equation may require prohibitively fine discretizations Consider a one-dimensional prototypical so lute transport equation, 8c 8c 82c -+v--D-=O 8t 8x 8xx (1.1) where c is concentration, v is velocity, and D is diffusion. Define the dimen sionless Peclet number by vL Pe=D' where L is a characteristic length of the system, and the grid Peclet number by vL1x Peg=D' wher e L1x is the mesh size of the grid For mesh density such that L1x satisfies 1

PAGE 19

standard Eulerian numerical methods, such as centered finite differences or Galerkin finite elements, perform well. These methods produce nonphysical oscillations 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 L Pe n --.-< x -!::.x -2' 2L !::.x < -Pe' (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 ( X -vt) c(x, t) = er f v'4Dt where er f is the error function. The width of this front is proportional to VD, hence to ffe. Eulerian methods need a number of cells proportional to VISe on a front to satisfy !::.x proportional to Thus the number of cells on a moving front must increase as a problem becomes more advection dominated A method requiring !::.x proportional to ffe 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 !::.t multiplied by a higher time derivative of the analytic 2

PAGE 20

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 stepping 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 advectionor dispersion-dominated, often varying in different parts of the space-time 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 advection-diffusion transport equations. Eulerian methods are characterized by the use of a fixed spatial grid Petrov-Galerkin 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], transport-diffusion method [79], characteristic-mixed finite-element method [2, 101, 3], operator splitting method of [40, 100, 30], the Lagrangian-Galerkin method [76], and Eulerian-Lagrangian 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

PAGE 21

The Eulerian-Lagrangian 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 constant-coefficient one-dimensional equations, then extended to one-dimensional variable-coefficient 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 one-dimensional finite-volume setting. One-dimensional nonlinear advection-diffusion equations have been treated by Ewing [33] and Dahle, Ewing, and Russell [31]. Addi tional work has been done on ELLAM schemes for one-dimensional 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 one-dimensional case [82]. Various two-dimensional 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 three-dimensional ELLAM in a framework within which all characteristic methods can be viewed. In this thesis, a three-dimensional 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 method's apparent robustness is discussed in a one-dimensional context. 1.2 Three-Dimensional ELLAM A three-dimensional finite-volume Eulerian-Lagrangian 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

PAGE 22

MODFLOW-96, the USGS three-dimensional, finite-difference, ground-water flow model [52], [54], [53]. The code is written in FORTRAN-77 ELLAM [26] solves an integral form of the solute-transport equation. The code uses an implicit method for dispersion calculations, which allows for large time steps without stability constraints. The Eulerian-Lagrangian approach involves tracking mass through time, then solving a dispersion equa tion on a fixed-in-space grid. This is particularly advantageous for advection dominated problems, as are typical of many field problems involving ground water contamination, since Eulerian-Lagrangian approaches tend to generate less numerical dispersion than standard Eulerian finite-difference 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): + aeae) + a (cOVi) t t a;; L: O'W (1.4) +A(cO + g/J) =0 (summation over repeated indices is understood), where 0 is volumetric con centration (mass of solute per unit volume of fluid, ML-3), (2b 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, M M-1), c is the effective porosity (dimensionless), V is a vector of interstitial fluid velocity components (LT-1), D is a second-rank tensor of dispersion coefficients (L2T-l), W is a volumetric fluid sink (W < 0) or fluid source (W > 0) rate per unit volume of aquifer (T-1), 0' is the volumetric concentration in the sink/source fluid (ML-3), A is the reactive decay rate (T-1), t is time (T), and Xi are the Cartesian coordinates (L). Porosity is the ratio of pore volume to 5

PAGE 23

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 phase-equilibrium relationship in which C is proportional to C. The retardation factor is defined as: (]bC Rf = 1 + cC' An integral form of the solute-transport equation, which is a state ment of conservation of mass over the domain of integration, is the governing equation for this finite-volume ELLAM approach. Equation (1.4) is integrated against a space-time test function to provide the formulation, including treat ment of (cell or sub domain) boundary conditions and solute decay. The test function effectively specifies the domain of integration for the transport equation by the portion of the space-time 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: In loT (u + \7 (cCV cD\7C) L C'W + UACC) dtdx = 0, (1.5) where n 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 U a(cC) = a(ucC) au cC at at at and u 1 1 -\7. (cCV cD\7C) = -\7. (ucCV -ucD\7C) --\7u (cCV cDVC) R f R f Rf 6

PAGE 24

to yield the global equation, { {T (a(ucc) 1 1 in io at + Rf v (ucCV ucDVC) + Rf VU cDVC (1.6) u au v ) -I: C W cC(+ -. Vu AU) dtdx = O Rf at Rf The Eulerian-Lagrangian aspects of the m ethod derive from the re quirement that the test function satisfy the adjoint equation, fIJ: + Vu AU = 0 Thus, for the time step from tn to tn+1, we choose u of the form u(x, t) = f(x, t)e->'(tn+1-t), where ?fit + v f = 0, so that f is constant along characteristics of the retarded interstitial velocity field. Note that with u = e->.(tn+1-t) (that is, f = 1) in the following, we arrive at a statement of global conservation of mass over the time step: In(ucC)n+1dx In (ucc)ndx + In Itt;:+l if V (ucCV ucDVC) + if VU cDVC _..1!... EC'W dtdx Rf = O. To obtain a system of equations, each representing mass conservation on a cell, we use u = El Ul where l is the index for the finite-difference cells Ol in the transport sub domain, and local space-time test functions are defined with !l(X, tn+1) = 1 on Ol and !l(X, tn+1 ) = 0 elsewhere: { e->.(tn+l_t) on characteristics from any (x, t) E 0 x (tn, t n+1 ) Ul = into Ol at time level n + 1 o otherwise. (1.7) 7

PAGE 25

We thus arrive at a system of local ELLAM equations, (1.8) =0 where a 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; nj is the pre-image at time tn under advection of cell nl at time t n+\ rn+1 an x (tn, t n+1) is the space-time 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 space-time 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 right-hand side of the system of equations. 1.3 Formulation of ELLAM equations ELLAM approaches the hyperbolic-parabolic 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 finite-difference code, MOD FLOW. 8

PAGE 26

This ELLAM method approximates total solute flux across the trans port sub domain 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. User-input inflow con centrations, together with the outflow solutions, then appear on the right-hand 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, In/ (cc)n+l dx f:lt Ian/ if (cDV'c)n+l n ds = e-ALlt Inj (cc)n dx jjrr _A(tn+l_t) C v d d J e c inflow Ii: n t s supp u/nrn+1 f (1.9) + II e-A(tn+l-t) 2:= c': dtdx, suppu/nsupp W f where ni means the pre-image in the spatial domain at t n of nl at tn+l, rn+1 -an x (tn, tn+l) is the space-time boundary at time step n + 1, n is the unit outward normal vector, and supp f = {xlf(x)::f O}. 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

PAGE 27

234 5 678 9 10 sOurce 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 if Itt;:+l V (ucCV -ucDVC) dtdx = rtn+1 J. 1 ( D C) Jtn on Rf UcCV -Uc V n dsdt. Considering just the outflow portion of the boundary, this becomes, r 1 ( C) Jtn -R U cCV UcDV n dsdt ::::; U" outflo w f (1.10) Itt;:+! (UCCoutflOWRV) n dtds, u" outflow f where total flux across the boundary is now approximated by advective flux We index the outflow boundary faces with II and define the following test functions: Ull = o on characteristics from n at time level n into boundary area (On)ll at any time during time step otherwise { e-).(tn+l_t) (1.11) The mass across outflow boundary face II 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 fa c e II during the time step, as illustrated in figure 1.2 is given by J. rtn+l -).(tn+!-t) C V (On)ll Jtn e c outflow Rf n dtds 10 (1.12)

PAGE 28

1 2 3 4 5 6 7 8 9 10 source Figure 1.2. Outflow boundary receives advected mass from inflow, source in cell 1, and storage in all cells jjrr e->.(tn+l_t) cO-II Y. n dtds m ow Rf supp ulln8ninflow X (t n ,tn+l) + II e->.(tn+1-t) L: C': dtdx SUppUllnSUppW x (tn,t n + 1 ) f where (8n)1l is a boundary face, and (8n)l/ is the preimage at time tn of (8n)1l x [tn, tn+l]. This system of equations will be solved for Coutllow' 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 steady-state flow condition. The specific discharge, or flow 11

PAGE 29

rate per unit area, across each face of every finite-difference 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 finite-difference grid, the integrals on the right-hand 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 piecewise-constant 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)j Z = z(t)] along which the fluid is advected: dx Vx (1.13) --dt Rj dy -Vy (1.14) dt Rj dz 11. (1.15) --dt Rj This is accomplished by introducing a set of moving points that can be traced within the stationary coordinates of a finite-difference 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 pre-image 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

PAGE 30

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 doesn't 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 right-hand 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 rll will denote the spatial finite-difference grid in dexing, with this index order indicating an x, y, z sequencing. A right-handed 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: 6.t R1 (cD\7Ct+l n ds = Jonl C f (1.16) where m = 1,2,3 is the summation index for the dispersion term. Finite difference approximations to the space derivatives in the dispersion integral 13

PAGE 31

t"n+l n Advection of known I1l3SS d istributio n to new tirru= level ( C oUTMt=I). Figure 1.3. Preimage of cell may be irregularly shaped and not easily delimited by backtracking. Instead, the known mass distribution at time tn is tracked forward along streamlines of the advective flow to time tn+l. 14

PAGE 32

are calculated using centered differences as in MOC3D, generalized for varying grid dimensions. (See [71], page 64, for the form of the cbDm expansions.) 1.4.4 Storage at New Time Level The quantity mass/porosity in a cell at the new time level tn+l 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 porosity (1.17) 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 right-hand 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: 15

PAGE 33

x ( ( .6.Yi ) + ( .6.Yi ) 4) .6.Yi-l + .6.Yi .6.Yi+l + .6.Yi X (( bj,i,k ) + ( bj,i,k ) 4)Cj ,i,k bj,i,k + bj,i,k-l bj,i,k + b j ,i,k+1 (( .6.Xj )( .6.Yi )(" bj,i,k )C.. A A A A b W.Xj_l + w.Xj W.Yi-l + W.Yi bj,i,k + j,i,k-l ( .6.Xj )( .6.Yi )( bj,i,k )C.. + A A A A b b W.Xj_l + w.Xj W.Yi-l + W.Yi j,i,k + j,i,k+1 ( .6.Xj )( .6.Yi )( bj,i,k )C .. + A A A A b b W.Xj-l + w.Xj W.Yi+l + W.Yi j,i,k + j,i,k+1 ( .6.Xj )( .6.Yi )( bj,i,k )C.. + A A A A b W.Xj-l + w.Xj W.Yi+l + W.Yi bj,i,k + j,i,k-l ( .6.Yi )( .6.Xj )( bj,i,k )C .. + A A A A J+l,2-1,k-l W.Yi-l + W.Yi W.Xj+l + w.Xj bj,i,k + bj,i,k-l ( .6.Yi )( .6.Xj )( bj,i,k )C.. + A A A A b b W.Yi-l + W.Yi W.Xj+l + w.Xj j,i,k + j,i,k+l ( .6.Xj )( .6.Yi )( bj,i,k )C.. + A A A A b b J+l,2+1,k-l w.Xj+1 + w.Xj W.Yi+l + W.Yi j,i,k + j,i,k-l ( .6.Xj ) ( .6.Yi ) ( bj,i,k ) C.. ) + A A A A b J+l,2+1,k+l W.Xj+1 + w.Xj W.Yi+1 + W.Yi bj,i,k + j,i,k+l + (( .6.Xj ) + ( .6.Xj ) 4) .6.Xj_l + .6.Xj .6.Xj+1 + .6.Xj .6.y. b k X (( 2)( J,t, )Cj,i-l,k-l .6.Yi-l + .6.Yi bj,i,k + bj,i,k-l + ( .6.Yi ) ( bj,i,k )C A + A b + b J,2-1,k+l W.Yi-l W.Yi j,i,k j,i,k+l + ( .6.Yi ) ( bj,i,k )C A A b + b J,2+1,k-l W.Yi+1 + W.Yi j,i,k j,i,k-l .6.y. b k + ( .6. (b ) C j ,i+1,k+1) Yi+1 Yi j,i,k j,i,k+1 + (( .6.Yi ) + ( .6.Yi ) 4) .6.Yi-l + .6.Yi .6.Yi+l + .6.Yi .6.x b k X (( J)( J,2, )Cj-1,i,k+l .6.Xj_l + .6.Xj bj,i,k + b j ,i,k+1 + ( .6.Xj )( bj,i,k )C A + A b + b J+l,t,k-l W.Xj+l W.Xj j,i,k j,i,k-l 16

PAGE 34

17

PAGE 35

x((( ) + ( ) 4)(( )Cj,i-l,k + + + + ( )Cj,i+!,k) (( ) + ( ) -4) + + + x )Cj+l,i,k)). Xj-l Xj Xj+! Xj 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 right-hand 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 sub cells determined by parameters NSC, NSR, and NSL, specifying the number of sub cells 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 mayor may not also be found in that destination cell. In order to mitigate the effects of unwarranted mass lumping, sub cell mass is distributed among cells neighboring the destination cell using the "approximate test functions", WI, described below. The value of WI at the sub cell center destination point is the fraction of sub cell mass to be distributed to cell [ll. This yields the formulation, e-,\.lt fni cn dx = e-,\.lt L L WI (pl)C(p) j,i,k p= ... beell center (1.18) where summation runs through all sub cells of each cell in the transport sub domain, and pI 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

PAGE 36

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 sub domain 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 one-dimensional approximate test function in each direction: Wjik(X, y, i) = f(x)g(Y)h(i). On a uniform grid, we define reference coordinates X, y, i with respect to a cell 01 with node indices j, i, k by x-x X J J and similarly for y and z. For an interior cell on a uniform grid with all surrounding cells active, f(x) = o NSCx + HNSC + 1) 1 -NSCx + HNSC + 1) o and similarly for g and h. These functions are shown graphically in figures 1.4 and 1.5. 19

PAGE 37

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 test functions corresponding to f defined above are illustrated graphically in figure 1.6. For a point in r2l with reference coordinates x, y, 2 E [f is defined generally, on a potentially nonuniform grid by, J-J J-J J-J -t:Sx<-t+ 2N1SC f(x)= 1 2NSC(( )-l)x+( )+NSC(l-( )) 1 __ 1_ 0 their values given by, Wj+l,ij = (1 f(x)) g(y) h(2), Wj,i+l,k = f(x) (1 g(y)) h(2), Wji,k+l = f(x) g(y) (1 h(2)), WH1,i+l,k = (1 f(x)) (1 g(Y)) h(2), Wj,i+l,k+l = f(x) (1 g(y)) (1 h(2)) Wj+l,i,k+l = (1 f(x)) g(Y) (1 h(2)), Wj+l,i+l,k+l = (1 f(x)) (1 g(y)) (1 h(2)). Figure 1.6. Approximate test functions in one direction on a nonuniform grid, with NS=2. 20

PAGE 38

f For X, i), z in any other octant ofthe cell, there are also eight potentially nonzero test functions evaluated similarly. Thus equations of the form (1.19) for j, 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 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 sub-time 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 transportl subgrid depending on when in the transport time step the mass enters the system. At each sub-time 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 right-hand side of the local conservation of mass equation for cell 01 with the following integration. The 21

PAGE 39

domain of integration is all sources which intersect the space-time test function for cell nl so that all source mass that flows into nl during a time step is integrated. Multiple sources within a source cell are summed. II e->.(tn+1-t} l: C': dtdx suppu/nsuppw f _1_ '" '" ->.(tn+1_tm } t1Tm ( f) '" c' Qs (1.20) Cj,i,k L.i L.i L.i e NSC NSR NSL WI P L.i S R all p= m=O s= I sou.rces source source 3ubcell in center cell where summation runs through all sources in all subcells of every source cell in the transport sub domain; pI is the image of p under forward tracking to the new time level; Q s is the flow rate of source s in the source cell; t1T = or if m = 0 or NT; and tm = tn + ;Tt1t 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 spacetime test function for a cell nl intersected with any sink cells. To approximate, this term is formulated only if nl 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 approximatation contributing to both the left-and right-hand sides of the equation for sink cell nl with coordinates j, i, k. This gives 22

PAGE 40

II e->-(tn+l-t): L C'][ dtdx supp Ul nSUpp W f C 9. + '" Qs etRf R s= f sink in cell (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 sub-time steps determined by parameter NT. The composite trapezoidal rule is applied in time. At each sub-time 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 two-dimensional boundary face is discretized, while for a source, the entire cell is discretized. For a cell ill, the integration is performed over the intersection of the space-time test function for that cell and the transport subdomain boundary; that is, all mass entering through the boundary and advected to ill during the time step is accumulated to the right hand side of local equation l Mass advected into ill during a time step from the inflow portion of the transport sub domain boundary is give by, (1.22) 23

PAGE 41

where ref area = NBC NBR, NBC NBL, ... depending on plane of face; pf is the image of p under forward tracking to the new time level; Q is inflow rate; tm = t n + :FTtlt represents the time during time step at which discretized source mass enters the system; and tlT = 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 nl 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, (1.23) :::::: tltQIl Coutflow where II is the index for boundary faces; and QIl 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 II is the un known in the boundary equation. The right-hand 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 right-hand side of (1.12) is accomplished using a trapezoidal in time, midpoint on sub cell in space approximation: 24

PAGE 42

e->..6..t tlXjtlYibj,i,k ( /)G( ) (NSG) (NSR) (NSL) Wll P P j,i,k p= subcell center +_1_ ->..(tn+1_tm ) tlTm ( f) G' Qs ej,i,k e NSG NSR NSLWI P sR all p= m=O s= / BOU1'ces SOUTce source subcell in center cell (1.24) 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 nl has been replaced by WIl, 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 Gout/low. 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->"t, 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 half-life is on the order of or smaller than the transport time step, however, some accuracy will be lost. 25

PAGE 43

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 ground-water 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 n e w tim e 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 sub domain 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 three-dimensional ELLAM code was tested by running the same suite of test cases as was applied to the USGS three-dimensional 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

PAGE 44

applications. 1.5.1 One-Dimensional Flow The first test case is a simple one-dimensional system involving solute transport in a finite-length aquifer with flow of constant velocity. Boundary conditons are third-type, 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 0'. = O.lcm. 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 0'. = 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

PAGE 45

Parameter Txx = Tyy C CY.L CY.TH = CY.Tv PERLEN (length of stress period) Vx Vy = Vz Initial concentration (GO) Source concentration (G') Number of rows Number of columns Number of layers DELR (6.x) DELe (6.y) Thickness (b) Value 0 .01 cm2/s 0.1 0 1 cm 0 1 cm 120 sec 0.1 cmls 0.0 cmls 0.0 1.0 1 122 1 0.1 cm 0 1 cm 1.0cm Table 1.1. Parameters used in ELLAM simulation of solute transport in a one-dimensional, steady-state flow system. 28

PAGE 46

0 8 ...... ........ . ................... . . .... ...... ..... 0 .6 ...... .... ....... ... ..... ...... ........ 0.4 ....... ..... ............ ..................... . ................ ........ 0 2 10 20 30 40 50 60 70 80 90 100 110 120 Figure 1.7. Plots of concentration as a function of cell node for one-dimen sional flow with a constant velocity field and low dispersion. Shown are the analytical solution (lowest graph), ELLAM 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

PAGE 47

. 0.8 ...................................... 0 6 .... .... .... ..... .. :.. ... .. .. ... .................... .. .. ...... 0.4 0 2 .............................. ....................... 10 20 30 40 50 60 70 80 90 100 110 120 Figure 1.8. Plots of concentration as a function of cell node for one-dimen sional flow with a constant velocity field and high dispersion. Shown are the analytical solution (lowest graph), ELLAM 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

PAGE 48

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 see1 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 method-of-characteristics 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, Three-Dimensional Flow Next, ELLAM results were compared with the analytical solution developed by Wexler [99] for three-dimensional 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 cross-terms 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 ELLAM results are plotted in figures 1.12, 1.13, and 1.14 for the xy plane passing through the 31

PAGE 49

O B ..... ................... ......... 0 6 ......... .. ......................... 0 4 .... ...... ....... ..... ...... ....... .. ... ; ..... 0 2 10 20 30 40 50 60 70 80 90 100 110 120 Figure 1.9. Concentration vs. cell node plot with CELDIS = 61 (two time steps), NSC = 4, NSR = NSL = 2, NT = 128. 1 ...... :.. ......... ..... ................. .. ..... : ........ : . : 0.8 ---.---. .......... ......... ........ 0 6 0 4 ... ... : ... ... ..... : .... : .... : . : .. 0 2 10 20 30 40 50 60 70 BO 90 100 110 120 Figure 1.10. Plots of concentration as a function of cell node for decay constant A = 0.01 sec-I Shown are the analytical solution (lower graph) and ELLAM results using CELDIS = 1 (121 time steps), NSC = 32, NSR = NSL = 2, NT = 128. 32

PAGE 50

Parameter Txx = Tyy C CiL CiTH CiTv PERLEN (length of stress period) Vx Vy=Vz Initial concentration (Co) Source concentration (Cf ) Q (at well) Source location Number of rows Number of columns Number of layers DELR (.6.x) DELe (.6.y) Thickness (b) Value 0.0125 m2/day 0.25 0.6 m 0.03 m 0.006 m 400 days 0.1 m/day 0.0 m/day 0.0 2.5106 g/m3 1.010-6 m3/d row 8, column 1, layer 11 30 12 40 0.5 m 3m 0.05 m Table 1.2. Parameters used in ELLAM simulation of transport from a con tinuous point source in a three-dimensional, uniform, steady-state flow system 33

PAGE 51

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 Two-Dimensional Radial Flow The ELLAM solution was compared to the analytical solution give by Hsieh [61] to a radial flow and dispersion problem with a finite-radius well in an infinite aquifer of two-dimensions. 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.161.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

PAGE 52

15 20 .... . ....... ; .... . 30 . ...... ...................... Figure 1.11. Concentration contours of analytical solution in the horizontal plane containing the solute source (layer 1) for three-dimensional solute transport in a uniform steady flow problem. 5 10 15 2 0 2 5 30 .. i .. . ; ....... ... ... ....... .. ...... . ........ .... .. ..... 25 ..... ; ..... . ..... . 30 .......... ........... .... .... ............................ .... .... . . Figure 1.12. Concentration contours in the horizontal plane containing the solute source (layer 1) for three-dimensional solute transport in a uniform steady flow problem with CELDIS = 7 (two time steps), NSC = NSR = NSL = 4, NT = 16. 35

PAGE 53

15 20 25 30 ...... ... .... . ... ..... . . .. ; ....... ......... ; ............. ; ...... ...... ; . . ... .... ... ...... ; ........... ; . 30 .......... ........................... .......... ........... ... .... Figure 1.13. Concentration contours in the horizontal plane containing the so lute source (layer 1) for three-dimensional 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 3 0 ... .. .... ... ; ...... ... ... ... ; ............ ... ........ ; .... ....... ; ........ ....... ........ 20 ... ..... ......... . 30 ........... .............................. .. ....................................... .. Figure 1.14. Concentration contours in the horizontal plane containing the so lute source (layer 1) for three-dimensional solute transport in a uniform steady flow problem with CELDIS = 0.1 (134 time steps), NSC = NSL = 4 NSR = 8, NT = 16. 36

PAGE 54

Parameter Txx = Tyy C CiL CiTH = CiTv PERLEN (length of stress period) Q (at well) Source concentration (CO) Number of rows Number of columns N umber of layers DELR (.6.x) = DELC (.6.y) Thickness (b) Value 3.6 m2/hour 0.2 10.0 m 10.0 m 1000 h o urs 56.25 m3/hour 1.0 30 30 1 10. 0 m 10.0 m Table 1.3. Parameters used in ELLAM simulation of two-dimensional, steady-state, radial flow case. 10 2 5 1 51:, ...... .. 20 .... ... ....... .. : .... ...... ...... : ...... . -. 25 ........ 30 .......... ........... .. ............ .. Figure 1.15. Contour plot of analytical solution for two dim e nsional radial flow problem 37

PAGE 55

30 Figure 1.16. Contour plot of concentration for run with two time steps using CELDIS = 75, NSC = NSR = 4, NSL = 2, NT = 16. 38

PAGE 56

NSC = NSR = 4, NSL = 2, NT = 16. The high concentration contour has non-physical 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 fiat, 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 fiat 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 three-dimensional solute transport from an instanta neous point source in a fiow 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 finite-difference 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 illequipped 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

PAGE 57

10 25 . 30 .................. . . . . . .................................................... Figure 1.17. Contour plot of concentration for run with 29 time ste ps using CELDIS = 5 NSC = NSR = 4, NSL = 2 NT = 4. 1 0 5 ............................. 25 ............. 3 0 ........................ .. Figure 1.18. Contour plot of concentr a tion for run with 563 time steps using CELDIS = 0 25, NSC = NSR = 4 NSL = 2, NT = 4. Concentration m a ximum is 1.019. 4 0

PAGE 58

25 30 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

PAGE 59

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 = O. For flow at 45 degrees to the grid, Vx = Vy = 1.0275, and Vz = o. 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 symmetry 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 grid-orientaion effect caused primarily by off-diagonal terms of the dispersion tensor. 42

PAGE 60

Parameter Txx = Tyy E: aL aTH = aTv PERLEN (length of stress period) Vx Vy=Vz Initial concentration at source Value 10.0 m2/day 0.1 1.0 m 0.1 m 90 days 1.0 m/day 0.0 m/day* 1 x 106 Source location column 11, in transport row 36, subgrid layer 4 N umber of rows 72 Number of columns 72 N umber of layers 24 DELR (6.x) 3.33 m DELe (6.y) 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 three-dimensional transport from a point source with flow in the x-direction and flow at 45 degrees to xand y-axes. 43

PAGE 61

10 20 30 40 50 6 0 70 10 ....... ................. : .. .. .......... : .................... 20 ................................ ... ..... ................ . . .... . ...... ...... 30 .............. ............... .. ........... : .. .............. 40 .. .. .. : ..... ....... I .. ........ I 50 .......... : ............. : .. .... .......... .. ..... ....... ........ M.. .. .. ............ ............. .. 70 Figure 1.20. Contour plot showing log of concentrations in analytical soluion of Dirac problem at t = 90. Concentration maximum is 25195 10 20 30 40 50 60 70 .. 10 ........... .............. ..... ....... .. : ... ... ...... ... .. 20 ................ .. ....................... .............................. .. 30 ............................................ ................. ............ 50 ........ : .............. : ..................... .... : .... .................... ...................... .. . 60 ......... : .............. ; .............. : ............................... ............ ................. 70 ....................... .. .............. .... ......... .......... .... ................... .. 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

PAGE 62

10 2 0 30 40 50 60 70 10 ..... ...... . ........ ... 20 ... ............... ... ........ ... : .. ......... : ....... ... : ...... . . .. '.' . . ... .......... 30 ........................... 40 ...... .... ........... 50 .......... ; . ......... : . . ...... : ............. : .............. : 60 ................. .: . ..... .. .... Figure 1.22. Contour plot showing log of concentrations in analytical soluion of Dirac problem at t = 130. Concentration maximum is 14539 10 20 30 40 50 60 70 10 ........... : .............. : .... .. 20 ........... : .......... ... ... .......... ....... ... .. ........ ... ... .......... 50 60 .. ............... ............ ... ................ .......... ............. 70 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

PAGE 63

10 20 30 40 50 10 20 ....... .... ....... 40 .... . ....... ........ ... .. 70 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 10 ........................................... 20 .... 30 40 ... ...... ..... .. ....... ... ....... 50 60 ...... ........ 70 ........... : ... .......... : ............... ................... ........ ...... .... ....... Figure 1.25. Contour plot showing log of concentrations for dispersed initial condition, and CELDIS = 5, NSC = NSR = NSL = 4 NT = 2 a t t = 130. Concentration maximum is 8167. 46

PAGE 64

(a) 20 Finite-Element Model 201---mP Conlours : 0 1 to 0 9 40 80 t20 160 200 Figure 1.26. Contour plot showing concentrations of a two-dimensional 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 element 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 cell-centered finite-difference grid. Parameters used for the ELLAM simulation are given in table 1.5.5. The ELLAM 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 doesn't 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 two-dimensional grid in order to expand it to three dimensions Transport results for a vertical plane at the 47

PAGE 65

Parameter Value K 1.0 m/day c 0.35 O'.L 3.0 m O'.TH 0.10 m O'.TV 0.01 m PERLEN (length of stress 12,000 days period) Q (at well) 56.25 m3/hour Source concentration( G') 1.0 N umber of rows 1 Number of columnsl 141 Number of layersl 91 DELR (.6.x) 1.425 m DELe (.6.y) 1.0 m Thickness (b) 0.2222-0.2333 m lOne 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, steady-state, two dimensional flow system (described by Burnett and Frind, 1987). 48

PAGE 66

00 20 40 60 80 100 120 140 160 180 200 100 200 300 400 500 600 700 800 900 Figure 1.27. Contour plot showing concentrations of ELLAM 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

PAGE 67

(a) 30 Finite-Element Model 20 Contours: 0.1 to 0.9 40 .. 80 120 160 200 Figure 1.28. Contour plot showing concentrations of a three-dimensional 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 aTy = O.Olm and aTH = O.lm, 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 aTy = aTH = O.lm, 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 tim e steps), respectively. The contours are in close agreement. 1.6 One-Dimensional ELLAM Consider a one-dimensional problem with c = Rf = 1, and D 2: 0 and assume, for simplicity, a uniform grid and v = constant> O. Also, forego the use of" approximate test functions," that is, let where and 9 = h = 1. Wjik(X, y, z) = J(x)g(y)h(z) f(x) = U x <_1 -2 _1. < x < 1. 2 -2 l
PAGE 68

100 ....... ; ......................... ........ ........................... ;. ........ : 200 .... ".; ... ... : ........ ; ..... .............. ................... .. 300 ................ .. 400 500 600 ..... ......... .... ..... ;.. .. ...... ; ... ... : .. ... .. : .......... : .... 700 ........ ; ....... ,,; ......... : ... ..... : .. 800 .................. .. ... .: .. .... ... ; ... ... ... ;. ... . 900 .................. 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. (a) 3D Finite-Element Model 20 Contours: O."to 0 9 10 .80 120 Figure 1.30. Contour plot showing concentrations of a three dimensional finite element solution of the Burnett and Frind problem with high vertical transverse dispersivity. Contours shown are 0.1 to 0.9. 51

PAGE 69

Figure 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

PAGE 70

Neglecting sources and inflow, which do not affect stability the system of ELLAM cell equations can be expressed (1.25) where A is the coefficent 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 1 A=S[1,6,1] (1.26) denoting the sub-, main, and super-diagonal entries, respectively The matrix M of centered difference coefficients at the new time level IS M = b..t 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 = vf:. 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+l, 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

PAGE 71

n+I t n t i-I X x-Ill X x+112 i+I Figure 1.32. Mass at old time level is integrated exactly using cell centers and pre-images of cell boundaries as integration points, and applying the trape zoidal integration rule. 54

PAGE 72

tn+I and interpolating concentration values between nodes, the trapezoid on the left has area, and the area on the right is, Since and C+t) = + Cr) Ci + -Cr) CHI, the mass in cell i at the new time level is 6.x [( Cr)(i + Ci -I + + Cr)(i + Cr)(i + Ci + (( -C r ) (i -Ci+ 1] Thus, disregarding boundary conditions, the B matrix is tridiagonal with 1 Cr Cr2 B = S[I, 6,1] + 2[1,0, -1] + -2-[1, -2, 1], (1.28) denoting the sub-, main, and super-diagonal 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 I fractional partl ::; In (1.28), Cr in the second two terms is replaced by the fractional part, and the entire matrix sum is diagonal-shifted 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 three-dimensional ELLAM implementation. Application of the three-dimensional code to a one-dimensional problem on a uniform grid 55

PAGE 73

with Courant number a multiple of 2;' n = 1,2,3, ... likewise tracks and inte grates advected mass exactly, and yields a solution identical to that produced by A-I B, 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 (1.29) where ec denotes the vector of truncation (consistency) errors at each node Letting en denote the error cn c(, t n), and subtracting (1.29) from (1.25), we have the error propogation Because A + M is nonnegative definate, this implies with a different constant incorporated into ec This leads to n en+! = 2:((A + M)-I B)j .6.tec + ((A + M)-I B)n+IeO, j=O where eO denotes any initial error. For stability we then require, in some norm, for 0 n.6.t T, where the constant K is independent of .6.t, and T is simulation time. Then Ilen+111 O(K)llecll, and consistency yields convergence as .6.x, .6.t -t O. Stability depends on the properties of the matrix (A + M)-l B 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, 161 A + M = [8 -k, 8 + 2k, 8 -k]. 56 (1.30)

PAGE 74

The right-hand side matrix B becomes 1 B = 8[1,6,1], since with no advection, Cr = O. (1.31) 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, !::ltce = o (!::lx2!::lt + !::lt2). Proof: We use standard finite difference arguments. First note that integrating the Taylor expansion of a function f gives, lXi+1 / 2 f(x) dx X i-l/2 f(Xi) + !'(Xi)(x -Xi) + -Xi)2 + O((x -Xi)3 dx !::lXf(Xi) + 0 + -Xi)2 dx + 0 !::lxf(Xi) + .6.:3 !"(Xi) + o (!::lx3), where it is useful to separate the .6.:3 and O(!::lx3) terms. Then, we have Also, rn+1 Cxx cxx(x, t n+1 ) -it Cxxt(x, s) ds tn+1 Cxx(Xi, tn+1) + Cxxx(Xi, tn+1) (x -Xi) + O(!::lx2) -1 Cxxt(X, s) ds + + (. tn+1)( ) + O( /\ 2) (!::lx)2 Cxxx Xi, X -Xi uX tn+1 1 Cxxt(X, s) ds. 57

PAGE 75

Integrating then gives and Combining terms gives the estimate, O(ci!l-Ci-l) + Hci+1 ci) + Hcitl-ci+1) + + 2ci+1 citl) + + Dividing by 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 1 A + M = 1+ (k -8)[-1,2, -1] and (1.31) as 1 B = I --[-1 2 -1] 8 where I is the identity matrix. The matrix (A + M)-l B has eigenvalues 58

PAGE 76

where A is an eigenvalue of [-1,2, -1]. Since 0 :::; A :::; 4, (A + M)-l B has nonnegative eigenvalues less than or equal to one for any k O. Because (A+M)-lB is normal, this gives II(A+M)-lBlb :::; 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-I B satisfies a discrete maximum principle. For 0 < k < however, the operator need not satisfy a maximum principle; consider, for example, 3 x 3 A I B with k = 116 Here, (: i -14 195 0 ( 16130 4 112 = 91 679 4 -1 1559 ll2 574 '""4 56 "4 ) 5!S The absolute row sum of the second row is the matrix does not satisfy a maximum principle. 1.6.2 Limiting Case of No Diffusion Since this is greater than 1 Consistency analysis for the purely hyperbolic problem has been done by Stynes and Russell [84]. Consistency error of O( (.6.X)2 + .6.t) for the interior of the domain, and O(.6.x + .6.t) 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 [2 norm, for any fixed meshsize. Proof: Since for no diffusion, k = 0, (A + M)-lB = A-1B = I + A-I 0, -1] + C;2 [1, -2, 1]) = 1+ A-I ([1,0, -1] + Cr[l, -2, 1]). 59

PAGE 77

Then, since 0 :::; Cr :::; II(A + M)-1 BI12 = IIA1 BI12 We have :::; 1 + IIA-1112 (11[1,0, -1]112 + Crll[l, -2, 1]112) :::; 1+ :::; 1 + 2Cr1lA1112 < 1 + 4vLlt Llx = 1 + o (b.t) for fixed b.x. 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 left-hand side (1.26) and right-hand side (1.28) operators are Laurent and b.x > O. 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 -1 B is a Laurent matrix with norm of one, and note properites of its spectrum. Given a sequence of complex numbers, a doubly-infinite 60

PAGE 78

matrix with Z-n along the nth diagonal is called a Laurent matrix: Zo Z-l Z-2 Z-3 Z-4 Zl Zo Z-l Z-2 Z-3 Z2 Zl Zo Z-l Z-2 (1.32) Z3 Z2 Zl Zo Z-l Z4 Z3 Z2 Zl Zo The following theorem is well-known: 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 E L 00 (T) such that {zn} is the sequence of Fourier coefficients of z: Zn = _1_ r27r z(einO)e-inOd() y'2;IT io (n E Z), 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 (z). Let Loo := Loo(T) and L2 := L2(T). Laurent matrices represent multiplication operators on L2 with re spect to the orthonormal basis { I ino} y'2;ITe nEZ Define This multiplication operator is bounded if and only if Z E Loo, in which ca se, IIM(z)11 = Ilzlloo. (1.33) 61

PAGE 79

Consider the operator which maps a function to the sequence of its Fourier coefficients, if,. : L2 -t l2(Z), f H {f } ':I:' n nEZ The operator is bijective and (1.34) for f E L2. Then (1.35) and properties of C can be established by determining the analogous property ofM. Theorem 1.7 For the case of no diffusion and an infinite spatial domain, (1.36) with 4 (( 1)2 3 ) a-1b = 3 + cos(O) Or -"2 cos(O) + Or ezB + 4 Or2 (1.37) and iiC(a-1b)II = 1 for any Courant number Or E [0, Proof: The continuous function a E LOO(T) with Fourier coefficients com posing the diagonals of (1.26) is given by a(O) = + 6 + eiB) = + cos(O)). The continuous function b E LOO(T) with Fourier coefficients composing the diagonals of (1.28) is given by b(O) = vk [0 + C;2) e-iB + Or2) + + + C;2) eiB] 2 (1.38) -_1_ (Or -1) cos(O) + Or eiB + Q Or2 -y'2; 2 4' 62

PAGE 80

We thus have (a) and (b). Furthermore, since zero is not in the range of a, the inverse of (a) is (a-1), (see [14], Theorem 1.2), where -1 4V2Jf a = 3 + cos(O)' (1.39) From (1.35), it follows that (1.40) So (a-1)(b) = (a -1b). Multiplying continuous functions (1.39) and (1.38), we get the bounded linear Laurent operator on an infinite domain, with a-1b given by (1.37). From (1.35), (1.34), and (1.33), it follows that 11(z) II = Ilzlloo for all z E Loo. (1.41) Thus, to determine the norm of (a-1b), we find the maximum modulus of the continuous function a-lb. For p (4Cr2 1)2, q 2(4Cr2 + 1)(3 -4Cr 2), r 16Cr4 -8Cr 2 + 9, we have la-1bl2 = pcos2(O) +qcos(O) +r (3 + cOS(O))2 Then OlaaIJ1W (IJ) = [2r 3q + (q + 6p) cos( 0) 1 = (64Cr2(2Cr2 1)) (cos(O) 1) =0 for 0 such that sin(O) = 0 or cos(O) = 1. Thus the critical points for la-1bl2 are o and 1f when 0 E [0, 21f). By (1.37), for 0 = 0, 63

PAGE 81

and for f) = 1r, Thus the maximum modulus occurs at f) = 0 and is equal to 1. Then (1.41) gives IIC(a-1b)11 = 1 for any Courant number Cr E 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 E C : z -Ae is not invertible in the algebra}, where e is the algebra identity element. For z E Loo, further define the essential range of z as the set R(z) := {A E C : meas(t E T: Iz(t) -AI < c:}) > 0 for every c: > O}, where meas means the Lebesgue measure. The essential range of z is the spectrum of z as an element of the Banach algebra Loo, 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 spC(z) = R(z), whenever z E Loo. We can now consider features of the spectrum of C( a-1b) 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 la-1bl = 1, the spectrum of C(a -1b) is contained in the closed complex unit disk. Since a-1b is not (globally or locally) a complex constant, C(a-1b) has no eigenvalues. The entire spectrum is seen to be continuous. Also, 0 E sp(a-1b) only for a Courant number of Lemma 1.8 For a-1b given by (1.37) and Cr E [0, o E sp(a-1b) if and only if Cr = 64

PAGE 82

Proof: We have a-1b = 0 if and only if o = cos(O) (Cr 2 Cr + i) + CreiO Cr2 + = (cos(O) 1)Cr2 + (eiO cos(O))Cr + i cos(O) + For 0 =I 0, using the quadratic theorem, Cr cos(O)-eilJ ye2iIJ -2 cos(O)eilJ +cos2(O)-( cos( 0)-1)( cos( 0)+3) 2(cos(O)-I) Since Courant number is real, this is satisfied if and only if, i sin(O) = 0, (1.42) so sin(O) must be identically zero. Since 0 E [0,271"), this means 0 = 0 or 0 = 71". If 0 = 0, substitution into (1.42) yields 0 = 1, hence no solutions. For 0 = 71", 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 IIA-IBII :S 1 using the l2, A, or other norm equivalent to II 112 independent of matrix size would guarantee stability independent of 6.x. 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-I B, and A-I B constructed with AI,1 = An,n = to reflect Neumann boundary conditions. In all examples, IIA-1 BIIA is bounded by one for the tridiagonal case corresponding to the Dirichlet boundary condition. Numerical study of the eigenvalues of A-I B, 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

PAGE 83

n IIA BII IIA BIIA IIA BIIA tridiagonal tridiagonal A -A 7 1,1 -n,n -5 0.99798782758721 1.02731655221192 0.99367278911431 1 .04572629392131 10 0.99984912202045 1.02646524138853 0 .99945791347173 1 .04220454113547 50 0 .99999971570515 1.02646524324671 0.99999884309114 1 .04216137329937 100 0.99999998183770 1.02646524324671 0.99999992482548 1 .04216137329937 200 0.99999999885218 1.02646524324671 0.99999999520802 1 .04216137329937 300 0.99999999977242 1.02646524324671 0.99999999904715 1 .04216137329937 5 0.99911499209546 0.99811477727258 0.99578523861235 0.99365010768407 10 0.99998979236483 0.99994432525461 0.99962924605506 0.99948444126033 50 1.00002196618037 1.00001117880143 0.99999920055886 0.99999912922811 100 1.00002196618549 1.00001117886019 0.99999994803407 0 .99999994570318 200 1.00002196618549 1.00001117886019 0.99999999668714 0 .99999999661270 300 1.00002196618549 1.00001117886019 0 .99999999934125 0 .99999999933138 5 1.00154201321822 1.00009454693413 0.99998843721148 0.99996021692627 10 1.00154201699366 1.00009495089198 0.99999896202120 0.99999809688996 50 1 .00154201699366 1.00009495089972 0.99999999774353 0.99999999745662 100 1.00154201699366 1.00009495089972 0.99999999985328 0.99999999984424 200 1.00154201699366 1.00009495089972 0.99999999999065 0.99999999999036 300 1.00154201699366 1.00009495089972 0.99999999999814 0.99999999999810 5 1.00023639568756 1. 0000368 7264 714 0.99999981047593 0.99999933999554 10 1.00023639647505 1.00003687289972 0 .99999998298590 0 .99999996868398 50 1.00023639647505 1.00003687289972 0.99999999996301 0 .99999999995829 100 1.00023639647505 1.00003687289972 0.99999999999760 0.99999999999745 200 1.00023639647505 1.00003687289972 0.99999999999985 0.99999999999984 300 1.00023639647505 1.00003687289973 0.99999999999997 0.99999999999997 5 1.00006026622643 1.00001031695580 0.99999998815467 0.99999995869556 10 1.00006026665609 1.00001031868571 0.99999999893661 0.99999999804191 50 1.00006026665609 1.00001031868571 0.99999999999769 0.99999999999739 100 1.00006026665608 1.00001031868571 0.99999999999985 0 .99999999999984 200 1.00006026665608 1.00001031868571 0.99999999999999 0 .99999999999999 300 1.00006026665609 1.00001031868571 1.00000000000000 1.00000000000000 5 0.99799236275199 1.02487950575327 0.99368383249979 1.04147443273890 10 0.99984936308447 1.02391532781825 0.99945833836339 1 .03829256968434 50 0.99999971662918 1.02391531969142 0.99999884356986 1.03823646002098 100 0.99999998194544 1.02391531969142 0.99999992485559 1.03823646002098 200 0.99999999886528 1.02391531969141 0.99999999520993 1.03823646002098 300 0.99999999977629 1.02391531969141 0.99999999904753 1 .03823646002098 Table 1.6. Sections represent Courant numbers, Cr = 1/2, 1/3, 1/64, 1/500, 1/2000 99/200, respectively. 66

PAGE 84

Lemma 1.9 (Horn and Johnson [60], Lemma 5 .6.10) Let A b e an n x n com plex matrix and c > 0 be given. There is a matrix norm such that g(A) ::; IIAII ::; g(A) + c. 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 = >'Bx, 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.2S). In no case have eigenvalue bounds of one been forthcoming. Sharp bounds on and where g is the spectral radius, have not been analytically established. It is easy to see that the matrix equation Acn+l = Bcn with A given in (1.26), and B defined by (1.2S), results from a forward difference in time, centered difference in space discretization of the equation, Ct + -S-Cxxt Cr Cx Cxx = 0, where the coefficients are taken to be constant. With the opposite sign on the space-time 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 well-known (cf. [83]) that ELLAM methods, in general, requir e sufficient mesh-density for about four grid nodes on a front for accurate model ing of advected concentration. (Compare with Eulerian methods where number 67

PAGE 85

1 ......... . . .... ..... . ... ... .... . . .... ........ ....... 0 5 . ....... ............. . ................. . . ..... -0. 5 '----'----'-----'------''-----'---+---'-----'-----'---'--20 40 60 80 140 160 180 2 00 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 c ombination of chapeau functions, the behavior of an initial condition of nonzero concentra tion at one node only under repeated application of the no-diffusion operator was studied. Qualitative features of the solution are independent of the problem size. The following represent results for the one-dimensional 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 fig ure 1.34 An nons ymme tric, oscillatory graph is produced for C r E As C r -+ 0 the maximum decreases and the oscillations becom e more extensive. This behavior is shown is figures 1.3 5 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: fewe r applications of the 68

PAGE 86

0 8 0 6 0.4 0.2 _0.2'---,,'--L-_:'-_..L-_-'-_--'-_-'---_-'-_--'-_-'-_ 20 40 60 80 100 120 140 150 180 200 Figure 1.34. Initial and advected peaks with Cr = 1 0.8 . ..... .......... : . ... ; ........ : ... ...... : ..... : ......... : ...... : . . 0.6 ; .. .; . : ... ;. ; .. > 0.4 .... ..... .. .... ... .. ... : ......... : ........ -: .. ....... : .. 0.2 ...... ; .. ....... : .... ..... : ......... : ......... ; .. ....... ; ........ ,............. .. .... .. ........ -0.2 '----"2':-0 --4'c-0--S-!-:0--!:80--,:-L00 =------',2-0 Figure 1.35. Initial and advected peaks with Cr = i. 69

PAGE 87

0.8 0.4 0 2 --'4':-0 --,6'-:-0--8::':-0--' ... 00:---':-:-20:-----:-" '4-::-0 --:-"6::-0 Figure 1.36. Initial and advected peaks with Or = 3 1 2 .... ......... ... . .... .. 0 8 ... .. ; .. .. .. : ...... : ... . 0 6 0.4 ...... .. J- ... ',-" -............. ; 0 2 . . ........ . .... ... ...... .. ... 0 -0. 2 20 40 60 80 '00 1 2 0 '40 '60 '80 2 0 0 Figure 1.37. Initial and advected peaks with Or = lOA. 70

PAGE 88

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-1 B) Jr in cases where Cr = -. -t1-. Here, S is the diagonaleger shift by -1 operator, S = [l,O,OJ. Disregarding boundary effects, application of (S(A -1 B) Jr Y 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(A1 B)Jr 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 (S(A -1 B) Jr ) of various sizes. (This is also seen for ((A1 B)JrS, which first shifts, then advects; and also (S(i)(A-IB)Jr) and ((A-1B)JrS(i)), with advection and shifting of i cells.) The eigenvector 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 ill-conditioned. The entries with largest modulus remain near the beginning of any eigenvector. It is not possible to approximate the 71

PAGE 89

0.8 .......... . . ............. ...... ....... . ; . .... ; ... .... -... ..... .. ... \ . ..... . . ..... .. ...... .... ... ; .. .... ...... .; .... .. : . . Figure 1.38. Eigenvectors of 10 x 10 S(A-1 B) Jr with CT = eigenvalues 0.02512380295106 0.14348830209936 0.34537563664597 0.56442393631698 0.74559736660795 0.87006283785412 0.99708068200960 0.94405857162778 0.98197582016563 o coefficients 11.08455231527136 17.99738031282700 -6 21949567802237 -5.06404951515980 6.99680498309164 -2.97173348456930 -0.76361709637109 0.71286631865321 -1.50653675213281 o Table 1.7. Eigenvalue and coefficient of the resp e ctive eigenvector in eigen vector decomposition of spike initial condition at node seven on 10 x 10 grid. 72

PAGE 90

Figure 1.39. Initial condition and results using 10 x 10 (S(A -1 B) Jr f with T = 1,10,100,1000 and Cr = 73

PAGE 91

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 ill-conditioning of the matrix of eigenvectors. The ini tial condition and the results of 1,10,100,1000, and 100000 applications of S(A-l B) Jr 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 ill-conditioned for MATLAB 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 well-discretized 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

PAGE 92

1 ..... .... ................. 0 8 ...... ........ 0 4 .............. 0 2 -0. 2 -0. 4 0 L ----1.--.L.---L-----'-----'----'------'-------'40 Figure 1.40. Initial condition and results using 40 x 40 (S(A -1 B) J r )T with T = 1,10,100,1000,10000 and Cr = k. 0.8 O s 0.4 0 2 -0.2 _0.4L----'----'------'--_-'--_.L.----'-_-"-_---'-_-"------' o 100 1 2 0 140 180 200 Figure 1.41. Initial condition and results using 200 x 200 (S(A-l B) Jr f with T = 1,10,100,1000,10000 and Cr = k. 75

PAGE 93

0.4 0.2 _0.2L----'---"--'--'---'--------'-----'----'-----'------' o 10 15 20 25 30 35 40 Figure 1.42. Initial condition with 2 nodes on a front, and results using 40 x 40 (S(A-1 B) Jr f 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-1 B) Jr f with T = 1,100 and Cr = k 76

PAGE 94

1.2 0.8 0 6 ................... 0.4 """""'."'" ...... ; .. 0.2 -0.20'----'---:'::10--"----',----':25:---3"::-0----:3'=-5 ---'40 Figure 1.44. Initial condition with 4 nodes on a front, and results using 40 x 40 (S(A-l B) Jr f with T = 1,100 and CT = 1.2 ............ "" ................. "_ .. .... ,,.. .." ........... .. ........ .. :. ...... .; ........... 0.8 .......... .............. .. 0.6 .. 0.4 ......... . ... .......... 0.2 0 -0.2 0 5 10 25 30 35 40 Figure 1.45. Initial condition with 5 nodes on a front, and results using 40 x 40 (S(A-1B)Jr f with T = 1,100 and CT = 77

PAGE 95

the convergence of the method as both and 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 advection-diffusion equa tion can simulate the transient, three-dimensional 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 three-dimensional 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 non-physical 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 one-dimensional 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 sub dis 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

PAGE 96

1.2 1 0 6 0.4 0 2 400 50 100 15 0 200 250 3 0 0 350 Figure 1.46. Superposition of results showing initial condition and advected peak using Cr = Grids are 5 x 5, 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 1 2 0 8 0.6 0 4 0 2 _0.2L--_'___---'-_--'__ '--_-'-_ Figure 1.47. Superposition of results showing initial condition and advected peak using Cr = Grids are 5 x 5, 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 79

PAGE 97

difference cell at the new time level, could lead to improvements in accuracy and efficiency. 80

PAGE 98

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 locally-refined 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 non-multigrid 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. Aggregation-based 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 Pi+l from the coarse levell + 1 to fine levell is constructed based on the decomposition of fine-level nodes into disjoint sets, called aggre gates. The aggregates then determine the nonzero structure of the prolongation 81

PAGE 99

operator from the coarse level to the fine level. Each fine-level 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 A1 the coarse problems are constructed by the recurrence AIH = (PiHf ALP/H' For illustration, we consider the simplest example of the tentative prolongator for the one-dimensional Laplace equation discretized on a mesh consisting of nl = 3L 1nL nodes: 1 1 1 1 1 1 1 1 PIH = J3 1 1 1 Each column corresponds to disaggregation of one lRn1+1 variable into three lRn1 variables, nl = 3nl+l' So, P/+1 can be thought of as a discrete piecewise constant interpolation. Since the matrix Al = A is tridiagonal, the aggregated coarse-level matrices AI, 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 fine-level matrix, aggre gation suggests itself as a natural way to construct AMG methods. Note that 82

PAGE 100

for scalar problems, the prolongation operator from coarse level of dimension m to fine-level 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 sub optimality is that the disaggregated coarse-level 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 so-called tentative prolongator P/+l 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 If+l = SIP/+l used in the method. The hierarchy of coarse problems is then constructed by the recurrence (2.1) 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 HI-equivalent 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 so-called weak approximation property, which can 83

PAGE 101

be formulated as (2.2) Although the method greatly benefits from utilization of the smoothed transfer operators If+1' note that the approximation property under which convergence is proved is formulated in terms of the properties of the tentative prolongators, and the smoothing Sl involved in the definition of If+1 is present only in right hand side in the form of the spectral radius e(At). This makes verification of (2.2) easy in practice [89]. In this work we take advantage of the local nature of aggregation basis functions in treatment of problems with coefficient discontinuities. Our algo rithm uses a coarsening strategy that respects boundaries of high-coefficient subdomains. As a result, coarsening in the high-coefficient regions will be per formed independently of the rest of the domain. Eventually, a high-coefficient region may be represented by a single node on one of the coarse levels. Such node will subsequently be eliminated from the system and no longer consid ered in the coarsening. In order to accommodate this elimination process, the prolongation smoother Sl has to be carefully designed. We introduce changes to the prolongation smoother considered in [89], which will allow us to retain good computational complexity, and allow us to prove for the problems with coefficient discontinuities the same asymptotic rate of convergence as proved in [89] for problems with uniformly HI-equivalent forms. In our analysis, we will focus mainly on the case of a single high coefficient region. The method is first presented in an abstract setting, and analyzed for the case of a scalar elliptic problem. Theory applicable to prob lems with multiple high-coefficient regions is then developed in section 2.8. The key assumption of our theory in sections 2.7 and 2.8 is a modified weak approx imation property, which includes prolongator smoother Sl in its left-hand side. We show at the end of section 2 8 that our definition of Sl allows easy verifica tion, in spite of the presence of the smoother Sl in the left-hand side Finally, section 2.9 concludes with numerical experiments demonstrating the efficacy 84

PAGE 102

of our smoothed aggregation method applied to several model problems. 2.2 Abstract Convergence Theory In this section we present modification of the abstract theory used in [89] suitable for later application in the analysis of the problems with coef ficient discontinuities. Given a finest level nl x nl symmetric and positive definite matrix Al and recalling the definition of the prolongators of the form If+l = SIP/+ll one iteration of the smoothed aggregation multigrid x MG(x, b), solving (2.3) is described in abstract terms as the following variational multigrid algorithm Algorithm 2.1 Let Rl : IRnl -t IRnl, l = 1, ... L -1 be given smoothers and 1/, 'Y > 0 be a given smoothing and cycle parameter, respectively. Set MG = MGl where MG1(, .), l = 1, ... L -1 is defined by: Pre-smoothing: Perform 1/ iterations of RIAI)X I + R1bl Coarse grid correction: If l + 1 = L solve AI+!X I +! = b l +! by a direct method, otherwise set Xl+! = 0 and perform 'Y iterations ofxl+l MGI+!(xl+I, bl+l), Post-smoothing : Perform 1/ iterations of RIAI)X I + R1bl As noted in section 2.1, columns of the tentative prolongaror Pi+! will have disjoint nonzero structure. Here and in the following sections we assume 85

PAGE 103

that the columns of P/+l have been normalized, so that P/+l : lR,n/+ 1 -----+ lR,nl, nl = ord(A) > n2 > ... > nL is a full rank orthogonal matrix, (2.4) Given a prolongator smoother Sl : lR,n1 -----+ lR,nl, the hierarchy of coarse-level problems are constructed by the recurrence (2.1), so we have (2.5) Following the proof by Vanek, et al. in [89), we state an abstract convergence result for an AMG method based on smoothed aggregation with prolongator smoothers and tentative prolongators satisfying the requisite prop erties. The earlier result stipulated a weak approximation property to be sat isfied only by the (unsmoothed) tentative prolongators. Here, we incorporate the prolongator smoother into the weak approximation property, noting that the smoother may indeed have a nontrivial kernel. In fact, we will take advan tage of this in our analysis. The prolongation smoother will be constructed in a way which renders irrelevant for approximation selected entries of a coarse vector. We now introduce some notation. Define the composite tentative prolongator, P/ : lR,n 1 -----+ lR,nl, by and the smoothed composite prolongator Il : lR,n1 -----+ lR,n 1 by The transfer operators Il define a hierarchy of coarse spaces UL UL 1 ... U1 by Ul = Range Il, with the norm on Ul induced by the lR,nl-norm IIXlllR"l = (xT X)1/2, lIulll = min{llxlllRnl : u = Ilx}, 86

PAGE 104

and the associated inner product (u, V)l = (x, Y)IRn/. Note that by the con struction of coarse problems (2.1), (2.6) Our estimates are based on an abstract regularity-free convergence result proved in [17]. It can be written in our notation as follows: Lemma 2.2 (Bramble, Pasciak, Wang, Xu [17], Theorem 1). Assume there are linear mappings Ql : Ul H U1 Ql = I and constants Cl, C2 > 0 such that for all u E Ul and every levell = 1, ... ,L (2.7) for all u E Ul and every level l = 1, ... L -1 (2.8) Further assume that Rl are symmetric positive definite matrices satisfying with a constant CR > 0 independent of the level. Then, Algorithm 2.1 satisfies IIx -MG(x, b)IIA (1 Ilx -xiiA \Ix E Ul (2.9) where x is the solution of (2.3), and co(L) = (1 + Cl + C2CR)2(L -1). Moreover, the preconditioner P defined by the action of MG(O,) is symmetric with respect to (., )ntnl and cond(A, P) co(L). In the following lemma, Assumptions (2.7) and (2.8) of Lemma 2.2 are verified from the properties of Sl and Pl+1. 87

PAGE 105

Lemma 2.3 Let for every l = 1, ... L -1, .xl 2: e(AI) and be given linear operators. Assume that for some C1 C2 > 0 and all l = 1, ... ,L -1, I 2 < \lu E IRnl, (2.10) IISI(QI Pl+1QI+1)ulllRnl Al = 1, (2.11) IIStilAl < 1, (2.12) II (I SI)xllitnl \Ix E IRnl. (2.13) Then, for every u E U1 the mappings QI = flQI satisfy (2.14) with e1(l) = 1 + C1(l-1), and Proof: Let u E U1 From the definitions of Ql+1 and 1/+1' and (2.6), Then, using (2.12), (2.10), and (2.6), we have -1-< IISI(QI P I+1QI+1)uIlAl + IISIQIUllAl < e1 /2(AI)IISI(QI p/+1QI)ulllRnl + IIQIUllAl < e1 / 2(AI)C1.x;-1/21IuIIA + IIQIUIIA < C 111ullA + IIQIUIIA Estimate (2.14) follows by induction with Q1 = I. 88

PAGE 106

To prove (2.15), use the definition of the l-norm, (2.14), (2. 13), and (2.6) to get -1< II(QI SIPI+1QI+lulllRn / -1--IISI(QI Pl+lQI+l)U + (I SI)QIUlllRn / -1-< IISl(QI P1+lQl+l)ulllR n / + 11(1 Sl)QlUlllRn / C1 C2 < >-i/21IuIIA + e1 / 2(AI) IIQIUIIA/ < (Cl + C21IQtlIA) II II e1 / 2(AI) u A. Using (2.14) to see IIQtllA cl(l), we get (2.15). The following convergence result is immediate from above Lemmas 2.2 and 2.3. Theorem 2.4 Let for every l = 1, ... ,L -1, >-1 e(AI) and be given linear operators. Assume that for some C1 C2 > 0 and all l = 1, ... ,L -1, I 2 < \lu E IRnl, IISI(QI P 1 +lQl+l)ulllRn/ Al (p/+lfp/+l = 1, IISdIA/ < 1, II (1 SI)xllin/ < Ci II 1 2 e(A1 ) x lA/ \Ix E JRn/. Then, with Rl satisfying (2.9), Algorithm 2.1 satisfies Ilx -MG(x, b)IIA (1IIx XIiA \Ix E U1 where x is the solution of(2.3), with Cl (l) = 1 + C1 (l-1), c2(l) = C1 + C2C1 (l), and co(L) = (1 + Cl + C2CR)2(L -1). 89

PAGE 107

Moreover, the preconditioner P defined by the action of MG(O,') is symmetric with respect to (', and cond(A, P) ::; co(L). 2.3 Model Problem For clarity, we first present and analyze the AMG algorithm in the context of a scalar elliptic model problem with a single high-coefficient region. Consider the second order scalar elliptic problem, Find u E Vh such that a(u, v) = f(v) for every v E Vh (2.16) Here 0 C IRd d = 2,3 is a bounded domain; f E H-1(0); and a(,) is a coercive and bounded bilinear form on Hl(O) where a high coefficient applies on a single simply connected sub domain, Of:: a(u, v) = r \7u \7v dx + 12 r \7u \7v dx, in' c; ine with 0' = 0 \ Of: and c; 1. For concreteness, we assume a finite element discretization, with Th a quasiuniform finite element mesh on 0, and Vh a PI or Ql finite element space associated with Th. At some of the boundary nodes, a zero Dirichlet boundary condition is imposed for functions in Vh We assume the standard scaling of the finite element basis, II'PiliulO = 1. This yields a linear system, Ax = b. We will then consider the case where high coefficients apply over a number of simply connected sub domains Ok indexed by k. 2.4 Algorithm The algorithm depends upon the smoothed aggregations concept. Let L designate the coarsest level. On each level l, l = 1, ... L -1, nodes are organized in small, disjoint clusters called aggregates, where nlH is the number of aggregates on level l, or equivalently, the number of nodes on level l + 1. On the finest level, these clusters have to be specified, e.g., as the sets of degrees of freedom associated with the finite element vertices, while 90

PAGE 108

the coarse-level aggregates are created by our aggregation algorithm A simple greedy algorithm for generating aggregates based on the structure of the coarse level matrix is given in [87]. We introduce the composite aggregate which is the aggregate At understood as the corresponding set of degrees of freedom on the finest level. Formally, is defined by where = U Atl. Define the discrete l2-(semi)norm of the vector x E IRnl given by ( )1/2 Ilxllz2(A\) = I: dofs k of Ai (2.17) Note that since aggregates on each level make a disjoint covering of the nodes, nl = I: i=1 (2.18) Let the jth node on levell be denoted Wl,j = A;-I, j = 1, ... ,nl. The algorithm requires as input the kernel of the stiffness matrix obtained from the finite element model with no essential boundary conditions. Assume that we are given the system to solve, Ax = b, and the representation B of the zero-energy mode(s) of the bilinear form with respect to the finite element basis. For second order scalar problems discretized by standard linear Lagrange finite elements, B = 1, a vector of ones. In order to facilitate the analysis for problems with jumps in coefficients, instead of the original system Ax = b, we consider the equiva lent problem with diagonally-scaled finest-level matrix A1x = D-1 / 2b where D = diag(A) and Al = D-1 /2 AD-l/2. We note that as a consequence of this transformation the representation of the zero-energy modes changes to Bl = Dl/2B. Our coarsening algorithm starts by aggregating separately in the high coefficient region. The generalized aggregation coarsening method can be de scribed as follows. 91

PAGE 109

Coarse-level representations of BI and of the tentative prolongators P/+I can be constructed simultaneously by the following recursive procedure, starting from AI, BI : Based on the current level matrix AI, the nodes are decomposed into disjoint aggregates The aggregates determine the nonzero structure of the tentative prolongator P/+1. Using the knowledge of Bl and the disjoint nonzero structure of P/+1 the values of P/+1 are computed simultaneously with the coarse-grid representation of the kernel, B/+I by the recurrence [89] This construction is unique due to the normalization (2.4). As a result of separate coarsening of a high-coefficent region, this region will eventually be represented by a single node on one of the coarse levels. At finer levels, because of the diagonal scaling and separate treatment of high coefficient subregions, the tentative prolongators P/+1 obtained are identical to those in the uniform coefficient case. After the level where the high-coefficient area is first represented by a single node, the separate coarsening forces us to consider the node as a separate one-node aggregate. This is undesirable for two reasons. Applying prolongator smoother to a tentative prolongator with a single-node aggregate leads to extensive overlapping of coarse grid basis functions which mean coarse problem fill-in. Also, the presence of many high-coefficient regions might lead to a proportionately large coarsest grid size. To avoid this situation, single node aggregates will be eliminated from the system. These issues also force us to consider a modified prolongator smoother operator. While [89] considers prolongator smoothers of the form (2.19) we generalize here to a prolongator smoother of the form (2.20) 92

PAGE 110

where PA/,ul denotes a projection onto Ul orthogonal with respect to the Al -inner product, and Pv denote projection onto Vi orthogonal with respect to the Euclidean inner product. In order to describe this modification, we introduce for each level subspaces Ui of nodes to be eliminated and Vi of nodes to be smoothed over. Thus, if there is a node j* that represents a single-node aggregate, we set Vi = {x E IRn / : Xj> = O}. This choice makes PVz serve as a filter, and will guarantee that the value of Xj> will not be changed by prolongator smoothing This averts extensive overlapping of coarse-level basis functions which would otherwise result. If j* is a node to be eliminated, we set Ul = span {aej*, a E IR}. With this choice, the prolongator smoother performs a simple harmonic extension into a node being eliminated from its neighbors. This is necessary to guarantee that we can interpolate into all fine-level nodes after the node elimination. Application of this coarsening-smoothing strategy also requires an other concern be taken into account: we must guarantee that the energetic projection is stable in the Euclidean norm. The diagonal scaling implies that the coarse-level basis function corresponding to the single-node aggregate has very little energy. Because of this low energy, the corresponding diagonal entry in the coarse level matrix is small, and the Euclidean norm of the energetic projection may be large. To assure the stability of the energetic projection, provision is made to move a single-node aggregate to coarser levels, indepen dent of smoothing over the rest of the domain. The projection PVz serves to preserve the diagonal value corresponding to the one-node aggregate. Mean while, the spectral radius of the coarse-level matrices are decreasing. When the spectral radius of the coarse-level matrix is comparable to the small diagonal entry associated with the single-node aggregate, the high-coefficient node is eliminated. 93

PAGE 111

2.4.1 Stages of Processing The processing of a node associated with a high-coefficient sub domain is at one of four stages at each levell of coarsening. The creation of aggregates, the space VI which is smoothed, and the space Ul of nodes to be eliminated at the current level, all depend on the stage of processing of each high-coefficient region. Spaces VI and U1 are used to formalize the description of carrying, then eliminating, nodes, and generalize easily to the case of multiple high-coefficient subdomains. Stage 1. Separate coarsening in n E and in n \ nE; standard smoothing over the entire space. Aggregates A; are created so that if A; contains a node of ne all nodes of A; are in ne. Ul = 0, VI = lR nl Ends when there is a j such that A; contains all nodes in ne the closure of nE in n. Assume j = niH for convenience. Stage 2. Single node representing nE carried to coarser levels by means of single-node aggregate; extensive fill-in of coarse-level matrices prevented by means of "filter" PVr. Create aggregates covering the nodes 1 ... nl -1 (i.e skip the node which corresponds to ne.) Set niH = niH + 1 = {nl } (the" last aggregate" for convenience) Ul = 0; VI = {x E lRnl : xnl = O} to disable smoothing over {nl}. Ends when (Al+1)nl+l,nl+l 2: CL).Vr+ll where CL is a specified parameter. Stage 3. Coarsening eliminates single node aggregate representing ne ; smooth over entire space but with smoother kernel equal to the high-coefficient subspace on level l. Create aggregates covering the nodes 1, ... nl -1 (i.e. skip the node which corresponds to ne.) 94

PAGE 112

Uz = {aenr a E IR}; Vi = IRnr Stage at single level only Stage 4. Standard coarsening; standard smoothing. Partition nodes into non-overlapping aggregates {An Ul = 0; Vi = IRnr Ends at coarsest level, L Algorithm 2.5 Set Al = D-I/2 AD-I/2 BI = D I / 2 B. Set l = 1 Repeat: (1) Partition the active nodes into non-overlapping aggregates {An (2) Define an nz x nZH matrix P/H (3) Scale P/H columnwise and create BZH: For j = 1, ... nZH Set aj = Ileol j (P/+I)II Update col j (P/H ) r ; col j (P/H ) ) Set (Bl+I)j = aj End for ( 4) Define spaces Ul Vi c IRnr and set if i E A; otherwise. (2.21) where PAr,ur is the AI-orthogonal projector onto Ul P1IL is the projector onto Vi which is orthogonal in the Euclidian inner product, and >-1IL is an upper bound on the spectral radius of P1ILAlP1IL (see Remark 2.8). (5) Create AIH r (SIP/Hf AISZPl+I (6) Set l r l + l. Repeat until nl is small enough for Az to be treated by direct solver. 95

PAGE 113

We next show that the prolongator smoothers (2.22) and tentative prolongators (2.21) are such that the assumptions of Lemma 2.3 are satisfied. 2.5 Smoother properties Prolongator smoother (2.22) is simply of form (2.19) in case all pro cessing is in Stage lor Stage 4. The form (2.22) accommodates Stage 2 carrying a high-coefficient node to a coarser level without smoothing and the Stage 3 elimination of a high-coefficient node. Stage 2 reduces e(A/) while maintaining any diagonal entry associated with a single-node high-coefficient region. The following lemmas regarding prolongator smoothers (2.22) are stated and proved in sufficient generality for the case of multiple high-coefficient subdomains. Lemma 2.6 Let A be an n x n symmetric positive semidefinite matrix and U, V subs paces of IRn. Set w = 4/3 and, S = (1 PUA) (1 Pv A ) (I PUA), 'AV (2.23) where PU,A is the projection onto U which is orthogonal in the Hilbert space {IRn (A, .)}, Pv is the projection onto V orthogonal in {IRn (.,.)} and .\v 2: e(PvAPv). Then IISIIA 1, T 1e(PvS ASPv ) < gAv. If, in addition, U is nonempty and (Ax, x) > 2-mm // //2 C Av, xEU X for some C > 0, then 1 [ 1 4 e(A) VI) 1 II(I-S)xll < --+-----+IIxliA If; C 3 Av e(A) C 96 (2.24) (2.25) (2.26) (2.27)

PAGE 114

for all x E IRn and, IISII ::; (1 + (1 + 3 AV C AV (2.28) If U is empty, then (2.27) and (2.28) hold with b = O. Proof: Set SA = I -w 5.y 1 P v A and PO;A = I -PU,A' Since PO;A is an A -orthogonal projection, norm submultiplicativity gives (2.29) Let x E IRn. From the properties of the projection P v we have (Ax, P v Ax) = IIPv Axll 2 and ::; 5.vIIPvAxI12 Hence 2;v (Ax, PvAx) + 2 < 2;)IPvAxI12 + IIPv Axll 2 ;v (2 -w) IIPvAxl12 ::; The statement (2.24) now follows by (2.29). Next we prove (2.25). Let x E V. The well-known properties of orthogonal projections APO;A = (PO;Af APO;A' PyX = x, = P v and = Pv yield (ST ASx,x) IIP#,A (I w5.y1 P v A) P#,A II P#,A (I w5.y1 pv(p#,Af AP#,A) II P#,A P v (I w5.y1 pv(p#,Af AP#,APV) x 5.v (p(A) x, x), (2.30) where p is a polynomial 4 p(t) = (1 "3 t)2 t (2.31) and A is a projected and scaled matrix 1 .iT.i A = =-PV(P u A) APU APV' AV (2.32) 97

PAGE 115

Further, Denoting by o-(A.) the spectrum of A., and using e(A.) :::; 1, the spectral mapping theorem gives ( 4)2 (4)2 1 e (p(A.)) = rna]} 1 --t t:::; max 1 --3 t t = -. tEO"(A) 3 tE[O,l) 9 (2.33) Substituting (2.33) into (2.30) one gets (P S TASP ) (ASx, Sx) 1, e v v = T:JC IIxll2 :::; gAV' (2.34) proving (2.25) To prove (2.27) for nonempty U, we use the identity (1 PU,A)2 = 1 PU,A, rewrite S as w S = 1 PU,A 5.v (1 PU,A) P v A (1 PU,A) and estimate for x E lRn 11(1 S)xll < IIPU,AXII + ;v (1 + IIPu,AII) IIPvll IIAl / 211 11(1 PU,A)xIIA e(A) 1 < IIPU,A xII + w 5.v J e(A) (1 + IIPU,AII) IlxiIA. (2.35) By (2.26), it follows that for x E lRn IIPu,A xll2 :::; C2\V (A PU,A x, PU,A x) < C2\V and, by :::; e(A)llxIl2 one further gets liP, II < U,A C AV Substituting two estimates above into (2.35) completes the proof of (2.27) 98

PAGE 116

The statement (2.28) follows from (2.27) using inequalities If U is empty, then I PU,A = I and both (2.27) and (2.28) again hold with l-O c -. Remark 2.7 We will see from the bounds (2.27) and (2.28) that the parameter CL in the Stage 2 stopping condition must be chosen to balance conflicting requirements. If this parameter is extremely small, the constant C in (2.26) (2.28) will be close to zero; but a very large parameter value allows a larger value of the ratio Elimination of a single-node aggregate representing a high coefficient region may have to be postponed in order to stabilize the energetic projection which depends on The following lemma asserts how to calculate '\Vc used in the Stage 2 stopping criterion based on the value of '\Vl' In order to show that prolongator smoothers and tentative prolongators satisfy convergence assumptions, Al must be related to '\Vc' Remark 2.9 discusses how properties of the coarse level matrices, together with the Stage 2 stopping condition provide the necessary control over submatrix eigenvalues. Bounds on Al are then formulated in a lemma. Lemma 2.8 For ,\ Vc chosen as (1) 1-1 AVc:= 9" Av1 (2.36) where '\Vl is a given upper bound of e(PVt A1PVJ, satisfies (2.37) Proof: By construction, P/+l of Algorithm 2.5 maps VI+l into VI and (P/+lfP/+l = I. We have P/+lPVc+lx = PViP/+lPVc+1x and therefore, (PVc+! AI+lPVc+ 1 X, x) = (Sf AISIP/+l PVi+ 1 x, p/+l PVi+ 1 x) 99

PAGE 117

(pv,.sT AIS IPv,.P/+1Pv,.+l x, p/+1 Pv,.+l x) < (}(pv,.sT AISIPv,.)IIP/+1Pv,.HXI12, (2.38) where IIP/+1 Pv,.+l xl12 = ((P/+1)T P/+1 Pv,.+l X, Pv,.+l x) = IIPv,.H XIl2 ::s Ilxll2 Substituting this estimate into (2.38) gives, (2.39) The inequality (2.37) follows from (2.39) and (2.25) by induction. Remark 2.9 Let Av.-L and AUI be the nonempty blocks of matrix Al defined I by (1 -Pv,.)AI(I -Pv,.) and PUIAIPUp respectively. The purpose of the stopping condition in Stage 2 is to assure that on every level, (2.40) and at the same time, (2.41) In case of a single high-coefficient subdomain, Av.-L is a 1 x 1 matrix if \Il.l =I 0, I and AUI is a 1 x 1 matrix if Ul =I 0. Hence (2.40) and (2.41) are satisfied trivially. In the case of multiple high-coefficient subdomains, a straightforward gen eralization of Algorithm 2.5 applies the stopping criterion to each diagonal entry in Av/ independently, resulting in a AUIH which is a submatrix of AV:-L, with all of its diagonal entries approximately equal. Since AV:-L and I I AUI correspond to mutual interactions of high-coefficient subdomains, they are typically diagonal or sparse, and AUI is well-conditioned. Since AV[-L is 100

PAGE 118

Gram, the Gershgorin Theorem and the Cauchy-Schwarz Inequality together with the Stage 2 stopping condition guarantee (2.40) with G' = GCL. For ).min(Aul) Gmini(Aul)ii GcL.:xlt/, we have (2.41). Note that the stopping condition will be satisfied eventually for any CL > 0 since .:xlt/ decreases by a factor of 9 each level. To enforce (2.40) and (2.41) without indirect assump tions would require modification of Algorithm 2.5 and theory for Au/ to always be diagonal. Lemma 2.10 The largest eigenvalue of a coarse level matrix is related to the bound on the spectral radius of submatrix Alt/ by (2.42) and 1 2-G' ).1 :s; 91 1 (1 + G) ).Vl :s; 91 1 (2.43) for some positive constant G' dependent on A, the Stage 2 stopping condition, and finally, in (2.43), .:xVt. Proof: For any x E IRn/, so ).1/2 < .:x1 / 2 + .:x1/2 I It/ VIol' The bound (2.42) follows using (2.40). Substituting first (2.36), then then .:xVI into (2.42), one gets (2.43). Assuming no single node aggregates on the finest Lemma 2.11 Prolongator smoother (2.22) of Algorithm 2.5 given by W Sl = (1 PAz,U)(1 .:xl Plt/Al ) (I -PA/,u) for 1 = 1, ... ,L -1 satisfies convergence Lemma 2.3 conditions (2.12)-(2. 13) on the prolongator smoother. 101

PAGE 119

Proof: On each levell = 1, ... ,L -1, (2.22) is a smoother to which Lemma 2.6 applies Then, for every levell = 1, ,L -1 (2.24) assures (2.12). Also, Algorithm 2.5 produces spaces Ul and Vi with the following spectral properties (see Remark 2 9, (2.41), and (2.42)): either Ul = 0, or 2-. Amin(UI ) CAL, either Vi = lRnl and ::; 1, or Vi i= lRnl and ::; (1 + C')2. So if Ul i= 0, (2.26) applies with b = in (2.27). If Ul = 0, b = 0 in (2. 27). Then for any Ul and Vi, (2.27) and (2.42) give (2.13) with C2 = (1 + C') + (1 + C')2 (1: C' + ] Here, C' and C are the positive constants from (2.40) and (2.41), respectively. 2.6 Tentative Prolongator Condition (2.11) of Lemma 2.3 on the tentative prolongator is satisfied by construction. The remaining condition (2.10) is also satisfied merely by the properties of the aggregates, or equivalently, the tentative prolongator. Because of the separate coarsening in the highand low-coefficient regions, all finest level nodes of Wl,i are either in the high-coefficient sub domain or in a low-coefficient region if processing is at Stage 1 or Stage 2 on level l, in which case the coefficient applied to coarse-level basis functions to yield the matrix entry is e1 2 or 1. If at Stage 3 or Stage 4 on level l, then the high coefficient region is being eliminated or has already been eliminated, so the coefficient associated with all nodes involved in coarsening is 1 Algorithm 2.5 yields aggregates satisfying the following assumptions : 2.6.1 Assumptions on aggregates. For every level l and aggregate Ai, there exists a ball such that (1) All nodes of are located within (2) Every x E n belongs to at most K, balls 102

PAGE 120

(3) For each ball diam(BD Ch3l I (2.44) where h is the characteristic diameter of the quasiuniform finite element mesh Th. To discuss further properties of the aggregates, the following definition is required. Definition 2.12 We say that D c IRd is a shape-regular domain of characteristic diameter H if D can be mapped onto unit ball B C IRd using a WI,oo-diffeomorphism F such that almost everywhere on D. Here, of(x) denotes the Jacobian of F at Remark 2.13 The shape-regularity condition (2.45) is equivalent to Lipschitz conditions C2 IIF(x) -F(y) II H IIx -yll \Ix, y E D and IIF-I(x) F-I(y)11 IIx yll \Ix, y E B, or CI C2 H IIx -yll IIF(x) F(y) II H IIx -yll \lx,y E D. Proof: Follows by the mean value theorem for a CI-diffeomorphism, and its extension for a WI,oo-diffeomorphism using the density of C I in WI,oo. Algorithm 2.5 assures aggregates are one of the following two kinds: 103

PAGE 121

2.6.2 Kinds of aggregates (1) All nodes of belong to shape-regular nnE. (2) All nodes of are contained in n'. The domain n n' is shape regular. The shape-regularity assumption on these subdomains provides for use of the scaled Poincare inequality, stated as follows: Lemma 2.14 Let n be a shape-regular domain of characteristic diameter H. Then, for every U E HI (n) there is a constant q E JR such that (2.46) where Cp depends only on shape-regularity constants in (2.45). Proof: Since n is shape-regular, every U E H1(n) can be written as u(x) = u(F(x)), where F is a W1,oo-diffeomorphism satisfying (2.45). Hence, for every U E HI(n), (2.47) Letting x be an eigenvector of 8F, it follows from (2.45) that all eigenvalues A E a(8F) satisfy IAI E [%%], hence (deWF)-l = det(8F)-1 E ']. Substituting the last estimate into (2.47) and using (2.45) again gives (2.48) with constants c, C > 0 that depend exclusively on C1 and C2 in (2.45). By the same argument, one also gets CHdilu qlli2(B) < lIu qlli2(n) h(u -q)2det(8F-1 ) dB < CHdilu qlli2(B) V q E JR. (2.49) 104

PAGE 122

The proof is now completed using the Poincare inequality for a unit ball B[1]: where C is a Poincare constant for the unit ball. Indeed, setting q = q and substituting (2.48) and (2.49) into (2.50) gives (2.46). Remark 2.15 The expression on the left-hand side of (2.46) attains its mini mum when q is an L2(n)-orthogonal projection of u onto the space of constant functions. By well-known arguments this projection is given by (2.51) We recall several more well-known concepts used in the following proofs. For the set of finite element basis functions {'IIi}, define the finite element in terpolation operator II : lRn -7 Vh as IIv = ViCPi. Then, for stiffness matrix Ai,j = a(cpj,
PAGE 123

Since B = 1 for the model problem, (2.55) Also by construction, since the tentative prolongator P/+I formed when there is a high-coefficient area represented by a single-node aggregate at Stage 3 at levell has no column corresponding to that node, PBk = NiBi for every k > i. This means that PI Bk = Pl PBk = P/ NiBi whenever k > I, or (pI Bk)i = { (BI)i for i in low-coefficient region (2.56) o otherwIse. In the case of a single high-coefficient subdomain, we will now show that the tentative prolongators and prolongator smoothers constructed using Algorithm 2.5 on the model problem together satisfy the weak approxima tion property (2.10). The result is given in two lemmas: one showing the approximation properties of disaggregated coarse-level vectors over regions not being eliminated, satisfied using only assumptions on the aggregation; a sec ond showing that the smoothed disaggregated vectors indeed satisfy the weak approximation property. Lemma 2.16 For the single region of high-coefficient case, under Assump tions 2 6 1 and 2.6.2 on aggregates, for every u E IRnl, l = 1, ... L -1 we have, l 2 CA 2 IINl(Ql Pl+lQl+l)UlllRnl lIuIiD-I/ZAD-l/Z. with Nl defined by (2.54), and Ql = (Pl)T Proof: Since Pl is orthogonal, the lh column of Pl associates the lh degree of freedom on levell with the nodes of A;-l; and Nl is the identity except on span { enc} if it is the kernel of Nt, we have 106

PAGE 124

h M {nl if Nl = I were = nl -1 otherwise. Since PlQl is the orthogonal projector onto Range P/, and Range C Range Pl, Using the minimization property of orthogonal projection and (2.55), and let ting v = D-1/2U M IINl(QI QIH)ullinl < L Ilu j=l M Ln;in IIDl/2(v j = l The quasiuniformity of Th implies that Dii = Aii S; Caihd-2 where h is the charactenstlc mes lameter, an ai = e h d' d { 12 if WI i is in high-coefficient region 1 otherwise. Then we have M -12 '" d2 2 IINI(QI PIH QIH)ulllRnl S; C L.. aih IIV TilII12(wl,j)' j=l If II is the finite element projector, using (2.53), (2. 46), and (2.44), we have M IIN1(QI P/+1Q1H)ullinl < C L aih 2 l1IIv Tilli2(B I .nn/ ) j=l J M < L j = l J M < L j=l J 107

PAGE 125

Finally, the assumption of a bounded number of intersections 2.6.1-2, (2.52), and (2.43) give P/+1Ql+l)ullinl < < < Al concluding the proof. Lemma 2.17 Assumption (2.10) of convergence Lemma 2.3 is satisfied using tentative prolongators Pi+1 and prolongator smoothers (2.22) constructed using Algorithm 2.5 on the model problem, with QI = (Plf. Proof: Inequality (2. 28) holds for the prolongator smoothers (2.22) with the constant dependent on UI Using (2.42) and (2.41), we have for any (empty or nonempty) UI constructed by Algorithm 2.5, IISdl (1 + (1 + C1)2) (1 + + C1)) where C and c are constants from (2.40) and (2.41), respectively. By construction, NI = I except on Ker (SI)' Thus, using Lemma 2.16, -12 -12 IISI(QI Pl+1QI+1)ulllRnl IISINI(QI Pl+1QI+1)UlllRnl 2 -I 2 < IISdl IIN I(QI P I+1QI+1)ulhRn l < with C1 = JCA (1 + (1 + C1)2) (1 + + C1)). Convergence of the method for the model problem now follows from Lemmas 2.11, 2.17, and orthogonality of the tentative prolongator, by the abstract convergence Theorem 2.4 108

PAGE 126

Theorem 2.18 When used to solve the model problem on a domain with a single high-coefficient region, with components constructed by Algorithm 2.5, and using appropriate Rl Algorithm 2.1 converges, and satisfies where Al = D-1 / 2 AD-I / 2 and x is the solution of (2.3), with CI (l) = 1 + C1(l-1), c2(l) = CI + C2CI(l), and co(L) = (1 + Cl + C2CR)2(L -1), for MG as in Algorithm 2.l. Moreover, the preconditioner P defined by the action of MG(O, .) is symmetric with respect to (., )rrlnl and cond(AI' P) co(L). 2.8 MUltiple sub domains with high-coefficients As suggested in Remark 2.9, in the case of multiple subdomains with high-coefficients .ok each sub domain is treated independently by Algorithm 2.5. That is, when an individual subdomain is represented by a single node, it is processed without smoothing until its respective diagonal entry in Al is large enough in comparision to the (decreasing with coarsening) spectral radius of AVi. The node is then eliminated. The subspaces and V/.l and Ul are now potentially multi-dimensional. Define the index set of nodes eliminated at level l by 1ifof = {j : Wl,j at Stage 3.} Associated with a 1ifof, define the index set of high-coefficient areas eliminated at level l, by 1iP = {k : f2k is represented by the single node Wl,j for some j E 1ifof}. Algorithm 2.5 should yield aggregates with the property: For every j E 1ifof, there is a k tj. 1iP, If < l such that all nodes of A{ belong to f2k The domain B; n .ok is shape-regular. Lemmas 2.6 and 2.11 are formulated to include the possibility of mul tiple high-coefficient subdomains. It remains to be shown that the tentative prolongators provide the needed approximation properties. As with a single 109

PAGE 127

high-coefficient subdomain, the weak approximation property can be satis fied without reference to eliminated high-coefficient regions, as Lemmas 2 .19 and 2.20 will demonstrate. To eliminate a high-coefficient node, define a filter on level l by a diagonal matrix Nl where (2.57) Algorithm 2.5 applied to the model problem with multiple regions of high-coefficient yields the following environment: Dl BI+1 7\T Bl rl+1 HI = NIP/+1 pl = N1Pi NI_IP/-I (2.58) (2.59) (2. 60) These attributes of the tentative prolongator are used below to show an approximation property of the coarse spaces with respect to the energy norm of the scaled matrix AI. The bound is without reference to high-coefficient regions corresponding to eliminated nodes. Lemma 2.19 Assume there is a positive constant CA such that for every v E IRnl and all levels l = 1, ... ,L -1. (2.61) where wf {UWk,j, k = 1, ... l, j E k,j Then there is a sequence of linear mappings 110

PAGE 128

such that I 2 C 2 IINI( Ql Pl+1 Ql+1)UlllRnl ::; 5)luIiD-l/2AD-l/2 (2.62) holds for all u E IRnl, l = 1, ... L -1. Proof: For every levell = 1, ... L, define an nl xnl diagonal matrix (Nl)ii = {o for i E .Wl,k, 1 otherwIse. Clearly, NI = NI Further, we show that By the definition of Nl one gets for every x E IRnl, where ei is the lh canonical basis vector of IRnl. The nonzero structure of P/ then gives proving (2.63). To determine the range of P/, we first use (2.58) and (2.63): DIEI pI pi-lEI nl N El -I N(Dl E I I ) .Ll 1-1 I .Ll-l 1-1 1-1 .Ll-I Hence, by induction with Nl = NI it follows that 1 I 1 PI E = NI Nt-IE The product Nl ... Nt -l is an nl x nl diagonal matrix, {o for i E wl-I (Nl ... Nl-l)ii = 1 otherwise. 111 (2.63) (2.64) (2.65)

PAGE 129

Similarly, the zero structure of the rows of P/ is apparent from This is shown using (2.59) and (2.63) to give, PI pI pl-l pI liT pl-l lI-T pI pl-l I = 1-1 I = l-llYl-l I = lYI-l 1-1 I ; where by the same argument, nl 7\T nl pl-2 t d pI I .Ll-l lVl-2.L1_2 1-1' e c an 1 = For each composite aggregate define a vector W1,i E IRnl by { 1 -1-1 w1,i = (N1 ... N1 -1B)k for k E Ai o otherwise. (2.66) (2.67) Set WI to be an nl x nl matrix consisting of the columns w1,i i = 1 ... nl Then P/ei is the ith column of P/, and from (2.67), (2.64) and the nonzero structure of P"1, one gets Hence, Range WI C Range p? (2.68) For each level l < L, introduce a seminorm Since Dl/2 B = Bl by construction, and setting u = Dl/2v, we have from (2.65) and (2.67), for every r = (Tl' ... ,TnJT 112

PAGE 130

So, from (2.68) and (2.66) nl+l n;!n IIDI / 2(v > Let us set Ql = (Pl)T, l = 1, ... L. Since (P/)TP/ = I, the mappings PlQl are projections onto Range f>t1, orthogonal with respect to the Euclidean inner product. We now estimate using (2.69), well-known properties of orthogonal projections, Range C Range Pl, and (2.66), (2.63), and (2.59): nl+l IIDI / 2(v > 11(1 -Q/+r)lVl ... lV/ullinl II (1 -P/Ql)lVl ... lV/ulli n 1-1 2 +IIP l (Ql Pl+IQlH)N I ... NlUlllRnl > IIP/[(NI ... lVlf>tIf 1 -I T 2 -Pl+I (NI ... NIP IH ) ]ulllRnl IIP/[(p/ Nlf NIP/ H I -1 2 IIPI Nl(Ql PIHQl+I)ulllRnl 1 2 IINI (Ql -Pl+I QlH)UlllRnl Since = the conclusion now follows. For the model problem, the relevant kernel consists of the vector of ones. Algorithm 2.5 constructs tentative prolongators which satisfy the assumption of the preceeding Lemma 2.19. Lemma 2.20 Under Assumptions 2.6.1 and 2.6.2, there is a constant C > 0 such that for every u E lRnl and every levell = 1, ... ,L 1 it holds that 113

PAGE 131

Proof: Let be an aggregate of the first kind. Then, there is a domain Ok with coefficient Ck such that all nodes of belong to Ok' Hence Djj ::::; Cc-;:;2 hd -2 for all j E \ wf = and, IID1 / 2(V IID1 / 2(V < Cc-;:;2hd -21Iv < Cc-;:;2h2I1IIv r i lli 2 (B;nnv Choosing ri to be an integral average of IIv over n Ok, 1 1 ri = I (IIv)dx, meas(Bi n Of) Blnn' the scaled Poincare inequality (2.46) together with Assumption 2 6.13 gives Therefore, rr;!n IID1 / 2(V < < (2. 70) If is an aggregate of the second kind, all nodes j E \ wf are located inside the part of 0 where the coefficient is equal to one Hence IID1 / 2(V < Chd -21lv < Ch-211IIv rilli2 (Blnn')' Again setting ri to be an integral average of IIv over Bi n Ok, using the scaled Poincare inequality (2.46), and Assumption 2.6.13, we get rr;!n IID1 / 2(v < < (2. 71) From (2.70) and (2.71), the conclusion follows from the bounded intersections of the balls 2.6.12, (2.52), and (2.43) 114

PAGE 132

Lemma 2.17 is seen to hold for the multiple region of high coefficient case, although Nt can now have a multidimensional kernel. As for the case of the single high-coefficient subdomain, convergence of Algorithm 2 1 for the model problem now follows from Lemmas 2.11, 2 17, and orthogonality of the tentative prolongator, using abstract convergence Theorem 2.4 Thus we have: Theorem 2.21 When used to solve the model problem on a domain with mul tiple high-coefficient regions, using components constructed by Algorithm 2 5 and appropriate Rt Algorithm 2 1 converges, satisfying where Al = D-I / 2 AD-1 / 2 and x is the solution of (2.3), with C1 (l) = 1 + DI(l-1), C2(l) = C1 + C2C1(l), and co(L) = (1 + C1 + C2CR)2(L -1), for MG as in Algorithm 2.1. Moreover, the pre conditioner P defined by the action of MG(O,') is symmetric with respect to (', ')lRn1 and cond(A1' P) :S co(L). 2.9 Computational Experiments In this section, we attempt to demonstrate on several model examples the efficacy of the multilevel method based on the smoothed aggregation ap proach. The problem solved is (2.16) with domain n a unit cube. A Dirichlet boundary condition will be considered over the x = 0 face of the cube. All experiments were run using V-cycles. The iteration of the method was terminated when the Euclidean norm of the initial residual had been re duced by c = 10-12. In order to observe the behavior of the true error, zero right-hand side and a random initial approximation were chosen All problems were discretized on 125 subdomains. We show results for AMG used as a solver, and as a pre conditioner. The meshsize is varied so that in each table, each experiment is shown using 20 x 20 x 20, 40 x 40 x 40, and 80 x 80 x 80 elements. Experiments are shown with coefficients varying from 10-0" to 100" for each problem size, with 115

PAGE 133

Table 2.1. Checkerboard pattern, coefficients IOu, lO-u I dof. I a I iter I conv. ratio I iter PCG I conv. ratio I 9,261 0 12 8 847e-02 10 5.991e-02 9,261 3 11 7.026e-02 11 6.16ge-02 9,261 6 14 1.335e-01 12 9.377e-02 68,921 0 15 1.42ge-01 10 5.301e-02 68,921 3 12 8 551e-02 10 5.783e-02 68,921 6 13 1.130e-01 13 1. 110e-01 531,441 0 15 1.462e-01 10 5.980e-02 531,441 3 14 1.328e-01 12 8.281e-02 531,441 6 18" 2.104e-01 18 1.913e-01 the Laplace problem, a = 0, listed first in each case. The pattern of the coefficients is a checkerboard in table 2 .1, cubes touching at a node in table 2 .2, and randomly distributed in table 2.3. Convergence rates are the average of the rates measured on each level. We note that the method practically used did not strictly adhere to all the assumptions of the theory. The method tends to coarsen separately over regions with very different coefficients, although this is not explicitly enforced. 116

PAGE 134

Table 2.2. Two cubes touching at a node. Coefficient=10a in dark subregions (see Figure 2.1). dof. I (J I iter I conv. ratio I iter PCG I cony. ratio I 9,261 0 12 8.847e-02 10 5.991e-02 9,261 3 13 1.066e-Ol 11 6.61ge-02 9,261 6 14 1.344e-Ol 11 7.493e-02 68,921 0 15 1.42ge-Ol 10 5.301e-02 68,921 3 16 1.670e -Ol 11 6.354e-02 68,921 6 16 1. 720e-Ol 11 6.255e-02 531,441 0 15 1.462e-Ol 10 5.980e-02 531,441 3 15 1.440e-Ol 11 7.63ge-02 531,441 6 15 1.444e-Ol 11 7.874e-02 2.10 Conclusion The variant of smoothed aggregation AMG described here is shown to converge for the scalar model problem with a single high-coefficient subdomain. The method is also demonstrated to converge when applied to the problem with multiple high-coefficient 8ubdomains, under certain assumptions Com putational results are given, indicating that even without strict enforcement of the strategies presented here, smoothed aggregations AM G can perform well on problems with discontinuous coefficients. 117

PAGE 135

Figure 2.1. Configuration with two high-coefficient subregions touching at a single node. Table 2.3. Element coefficients random in (10-0",100"). dof. I a I iter I cony. ratio I iter PCG I cony. ratio I 9,261 0 12 8.847e-02 10 5.991e-02 9,261 3 12 9.871e-02 11 6.415e-02 9,261 6 12 9.871e-02 11 6.430e-02 68,921 0 15 1.42ge-0l 10 5.301e-02 68,921 3 15 1.511e-01 10 5.506e-02 68,921 6 15 1.511e-0l 10 5.530e-02 531,441 0 15 1.462e-Ol 10 5.980e-02 531,441 3 15 1.440e-01 11 7.365e-02 531,441 6 16 1.624e-Ol 11 6 921e-02 118

PAGE 136

References [1] R.A. Adams. Sobolev Spaces. Academic Press, New York, 1975. [2] T. Arbogast, A. Chilakapati, and M. F. Wheeler. A characteristic mixed method for contaminant transport and miscable displacement. In T. Russell et al., editor, Computational Methods in Water Resour ces IX, volume 1, chapter Numerical Methods in Water Resources, pages 77-84. Computational Mechanics Publications Publications and Else vier Applied Science, London, New York, 1992. [3] T. Arbogast and M F Wheeler. A characteristics-mixed finite element method for advection-dominated transport problems. SIAM J. Numer. Anal., pages 404-424, 1995. [4] N S. Bakhvalov and A. V. Knyazev. Fictitious domain methods and computation of homogenized properties of composites with a periodic structure of essentially different components. In Gury I. Marchuk, ed itor, Numerical Methods and Applications, pages 221-276 CRC Press, Boca Raton, 1994. [5] N. S. Bakhvalov and A. V. Knyazev. Preconditioned iterative meth ods in a subspace for linear algebraic equations with large jumps in the coefficients. In D. Keyes and J. Xu, editors, Domain Decomposi tion Methods in Science and Engineering, volume 180 of Contemporary Mathematics, pages 157-162. American Mathematical Society, Provi dence, 1994 Proceedings of the Seventh International Conference on Domain Decomposition, October 27-30 1993 held at the Pennsylvania State University. [6] N. S. Bakhvalov, A. V. Knyazev, and G. M. Kobel'kov. Iterative meth ods for solving equations with highly varying coefficients In Roland Glowinski, Yuri A. Kuznetsov, Gerard A. Meurant, Jacques Periaux and Olof Widlund, editors, Fourth International Symposium on Domain Decomposition Methods for Partial Differential Equations, pages 197205, Philadelphia, PA, 1991. SIAM. 119

PAGE 137

[7] N. S. Bakhvalov, A. V. Knyazev, and R. R. Parashkevov An efficient iterative method for solving Lame equations for nearly incompressible media and Stokes equations with highly discontinuous coefficients. Numerische Mathematik. Submitted. Published as a technical report UCDCCM 120, 1997, at the Center for Computational Mathematics, Univer sity of Colorado at Denver. (8] R. Bank, J. Burger, W. Fichtner, and R. Smith. Some upwinding tech niques for finite element approximations of convection diffusion equa tions. Numer. Math., pages 185-202, 1990. [9] R. E Bank, J. Mandel, S. F. McCormick, and R. E: Bank. Variational multi grid theory. In Multigrid Methods, pages 131-178. SIAM, Philadel phia, 1987. [10] J. W. Barrett and K. W. Morton. Approximate symmetrization and Petrov-Galerkin methods for diffusion-convection problems Comput. Methods Appl. Mech. Engrg., pages 97-122, 1984. [11] J. P Benque and J. Ronat. Quelques difficulties des modeles numeriques en hydraulique. Comput. Methods Appl Mech. Engrg., pages 471-494, 1996. [12] P. J. Binning. Modeling unsaturated zone flow and contaminant in the air and water phases. PhD thesis, Dept. of Civil Engineering and Op erational Research, Princeton University, Princeton,New Jersey, 1994. Ph.D. thesis. [13] P. J. Binning and M. A Celia A finite volume Eulerian-Lagrangian lo calized adjoint method for solution of the contaminant transport equa tions in two-dimensional multi-phase flow systems. Water Resources Research, pages 103-114, 1996. [14] Albrecht Bottcher and Bernard Silberman. Introduction to Large Truncated Toeplitz Matrices. Springer-Verlag, New York, 1999. [15] E T. Bouloutas and M. A. Celia. An analysis of some classes of Petrov Galerkin and optimal test function methods. In M. A. Celia et al., edi tor, Proceedings of Seventh International Conference on Computational Methods in Water Resources, volume 2, pages 15-20 Computational Mechanics Publications, Southampton, England, 1988. 120

PAGE 138

[16] E. T. Botiloutas and M. A. Celia. An improved cubic Petrov-Galerkin method for simulation of transient advection-diffusion processes in rect angularly decomposable domains. Comput. Methods Appl. Mech. En grg., pages 289-308, 1991. [17] J.H. Bramble, J.E. Pasciak, J. Wang, and J. Xu. Convergence estimates for multigrid algorithm without regularity assumptions. Math. Comp., 57:23-45, 1991. [18] A. Brandt, S. F. McCormick, and J. W. Ruge. Algebraic multigrid (AMG) for sparse matrix equations. In D. J. Evans, editor, Sparsity and Its Applications. Cambridge University Press, Cambridge, 1984. [19] Marian Brezina and Petr Vanek. A black-box iterative solver based on a two-level Schwarz method. Computing, 63:233-263, 1999. Dec 07, 1999. [20] A. N. Brooks and T. J. R. Hughes. Streamline upwind Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., pages 199-259, 1982. [21] R. D. Burnett and E. O. Frind. Simulation of contaminant transport in three dimensions, 2. dimensionality effects. Water Resources Research, 23:695-705, 1987. [22] C. S. Carrano and G. T. Yeh. A fourier analysis and optimaztion of the Petrov-Galerkin finite element method. In Alexander Peters et al., editor, Computational Methods in Water Resources X, pages 191-198. Kluwer Academic Press, London, 1994. [23] M. A. Celia. Eulerian-Lagrangian localized adjoint methods for contami nant transport simulations. In Alexander Peters et al., editor, Computa tional Methods in Water Resources X, pages 207-216. Kluwer Academic Press, London, 1994. [24] M. A. Celia and L. A. Ferrand. A comparison of ELLAM formulations for simulation of reactive transport in groundwater. In Wang, editor, Advances in Hydro-Science and Engineering, i(B), pages 1829-1836. University of Mississippi Press, Jackson, MS, 1993. 121

PAGE 139

[25] M A Celia, I. Herrera, E. T. Bouloutas, and J. S Kindred A new nu merical approach for the advective-diffusive transport equation. Numer Methods Partial Differential Equations, pages 203-226, 1989. [26] M. A. Celia, T. F. Russell, I. Herrera, and R. E. Ewing. An Eulerian Lagrangian localized adjoint method for the advection-dispersion equa tion. Advances in Water Resources, 13:187-206 1990. [27] M. A. Celia and S. Zisman. An Eulerian-Lagrangian localized adjoint method for reactive transport in groundwater. In G. Gambolati et al., editor, Computational Methods in Water Resources VII, pages 383-392. Springer, New York, 1990. [28] I. Christie, D. F. Griffiths, A. R. Mitchell, and O. C. Zienkiewicz. Finite element methods for second order differential equations with significant first derivatives. Internat. J. Numer. Methods Engrg., pages 1389-1396, 1976 [29] R. A. Cox and T. Nishikawa. A new total variation diminishing scheme for the solution of advective-dominant solute transport. Water Re sources Research, pages 2645-2654, 1991. (30J H. K. Dahle, M. S. Espedal, R. E. Ewing, and O. SlEvareid Character istic adaptive subdomain methods for reservoir flow problems. Numer. Methods Partial Differential Equations, pages 279-309, 1990. [31] H. K. Dahle, R. E. Ewing, and T. F. Russell Eulerian-Lagrangian localized adjoint methods for a nonlinear convection-diffusion equation. Comput. Methods Appl. Mech. Engrg., pages 223-250, 1995. [32] L. Dempkowitz and J. T. Oden. An adaptive characteristic Petrov Galerkin finite element for convection-dominated linear and nonlinear parabolic problems in two space variables Comput. Methods Appl. Mech. Engrg., pages 63-87, 1986. [33] R. E.Ewing. Operator splitting and Eulerian-Lagrangian localized ad joint methods for multiphase flow. In The Mathematics of Finite Ele ments and Applications, pages 215-232. Academic Press, London, 1991. 122

PAGE 140

[34] R. E.Ewing and H. Wang. Eulerian-Lagrangian localized adjoint meth ods for linear advection equations. In Comput. Mech. '91, pages 245250. Springer, Berlin, 1991. [35] R. E.Ewing and H. Wang. An Eulerian-Lagrangian localized adjoint method for variable-coefficient advection-reaction problems. In Wang, editor, Advances in Hydro-Science and Engineering, 1 (B), pages 20102015. University of Mississippi Press, Jackson, MS, 1993. [36] R. E.Ewing and H. Wang. An Eulerian-Lagrangian localized adjoint method with exponential-along-characteristic test functions for variable coefficient advective-diffusive-reactive equations. In U. Choi, D. Kwak, and J. Yim, editors, Proceedings of KAIST Mathematical Workshop. Analysis and Geometry, 8, pages 77-91. Teajon,Korea, 1993. [37] R. E.Ewing and H. Wang. Eulerian-Lagrangian localized adjoint methods for variable-coefficient advective-diffusive-reactive equations in groundwater contaminant transport. In Gomez and Hennart, editors, Advances in Optimization and Numerical Analysis, Mathematics and Its Applications, Vol. 275, pages 185-205. Kluwer Academic Publishers, Dordrecht, Netherlands, 1994. [38] R. E.Ewing and H. Wang. Ail optimal order error estimate for Eulerian Lagrangian localized adjoint methods for variable-coefficient advection reaction problems. SIAM J. Numer. Anal., pages 318-348, 1996. [39] K. Eriksson and C. Johnson. Adaptive streamline diffusion finite element methods for stationary convection-diffusion problems. J. Math. Comp., pages 167-188, 1993. [40] M. S. Espedal and R. E. Ewing. Characteristic Petrov-Galerkin subdo main methods for two-phase immiscible flow. Comput. Methods Appl. Mech. Engrg., pages 113-135, 1987. [41] R. E. Ewing and M. A. Celia. Numerical methods for reactive trans port and biodegradation. In T. F. Russell et al., editor, Computational Methods in Water Resources IX, pages 51-58. Computational Mechanics Publications, Southampton, England, 1992. [42} R. E. Ewing, T. F. Russell, and M. F. Wheeler. Simulation of miscible 123

PAGE 141

displacement using mixed methods and a modified method of character istics. Society of Petroleum Engineers of AIME, pages 71-81, 1983. [43] R. E. Ewing and H. Wang. Eulerian-Lagrangian localized adjoint meth ods for linear advection or advection-reaction equations and their con vergence analysis. Comput. Mech., pages 97-121, 1993. [44] J. Fish and V. Belsky. Generalized aggregation multilevel solver. Inter nat. J. Numer. Methods Engrg., 40(23):4341-4361, 1997. [45] Jacob Fish and Vladimir Belsky. Multigrid method for periodic hetero geneous media. II. Multiscale modeling and quality control in multidi mensional case. Comput. Methods Appl. Mech. Engrg., 126(1-2):17-38, 1995. [46] L. P. Franca, S. L. Frey, and T. J. R. Hughes. Stabilized finite ele ment methods: I. application to the advective-diffusive model. Comput. Methods Appl. Mech. Engrg., pages 253-276, 1992. [47] A. O. Garder, D. W. Peaceman, and A. L. Pozzi. Numerical calculation of multidimensional miscible displacement by the method of character istics. Soc. Petroleum Eng Jour., 4:26-36, 1964. [48] A. O. Garder, D. W. Peaceman, and A. L. Pozzi. Numerical calculations of multidimensional miscible displacement by the method of character istics. Soc. Pet .. Engrg. J., pages 26-36, 1964. [49] Karl E. Gustafson and Duggirala K. M. Rao. Numerical Range. Springer-Verlag, New York, 1996. [50] P. Hansbo. The characteristic streamline diffusion method for the time independent incompressible Navier-Stokes equations Comput. Methods Appl. Mech. Engrg., pages 171-186, 1992. [51] P. Hansbo and A. Szepessy. A velocity-pressure streamline diffusion finite element method for the incompressible N avier-Stokes equations. Comput. Methods Appl. M ech. Engrg., pages 107-129, 1990. [52] A. W. Harbaugh and M G McDonald. A modular three-dimensional finite-difference ground-water flow model. Techniques of Water Resources Investigations Book 6, Chapter A1, U.S. Geological Survey, 1988. 124

PAGE 142

[53] A. W. Harbaugh and M. G. McDonald. Programmer's documentation for MODFLOW-96, an update to the U.S. Geological Survey modu lar finite-difference ground-water flow model. Open-File Report 96-486, U.S. Geological Survey, 1996. [54] A. W. Harbaugh and M. G. McDonald. User's documentation for MODFLOW-96, an update to the U.S. Geological Survey modular finite difference ground-water flow model. Open-File Report 96-485, U.S. Geological Survey, 1996. [55] R. W. Healy and T. F. Russell. A finite-volume Eulerian-Lagrangian localized adjoint method for solution of the advection-dispersion equation Water Resources Research, 29:2399-2413, 1993. [56] R. W. Healy and T. F. Russell. Solution of the advection-dispersion equation in two dimensions by a finite-volume Eulerian-Lagrangian localize d adjoint method. Advances in Water Resources, pages 11-26, 1998. (57] C. I. Heberton, T. F. Russell, L. F. Konikow, and G. Z. Hornberger. A three-dimensional finite-volume Eulerian-Lagrangian localized adjoint method (ELLAM) for solute-transport modeling. Water Resources In vestigations Report OO-xxxx, U.S. Geological Survey, 2000. awaiting publication. (58] I. Herrera, R. E. Ewing, M. A. Celia, and T. F. Russell. Eulerian Lagrangian localized adjoint methods: The theoretical framework. Nu mer. Methods Partial Differential Equations, pages 431-458, 1993. [59] J. M. Hervouet. Application of the method of characteristics in their weak approximation to solving two-dimensional advection equations on mesh grids. In C. Taylor, J. John, and W. Smith, editors, Computational Techniques for Fluid Flow, Recent Advances in Numerical Methods in Fluids 5. SIAM Pineridge Press, Swansea, 1986. (60] Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, Cambridge, New York, and Oakleigh, Victoria, Aus tralia, 1985. [61] P. A. Hsieh. A new formula for the analytical solution of the radial dispersion problem. Water Resources Research, 22:1597-1605, 1986. 125

PAGE 143

(62] T. J. R. Hughes. Multiscale phenomena: Green functions, the dirichlet to-neumann formulation, subgrid scale models, bubbles, and the origins of stabilized methods Comput Methods Appl. Mech. Engrg., pages 387401, 1995. (63] T. J. R. Hughes and A. N Brooks. A multidimensional upwinding scheme with no crosswind diffusion. In Finite Element Methods for Convection Dominated Flows 34. ASME, New York, 1979. [64] T. J. R. Hughes and M. Mallet. A new finite element formulation for computational fluid dynamics: Iii. the general streamline operator for multidimensional advective-diffusive systems. Comput. Methods Appl. Mech. Engrg., pages 305-328, 1986. (65] Jr. J. Douglas and T. F. Russell Numerical methods for convection dominated diffusion problems based combining the method of charac teristics with finite element or finite difference procedures. SIAM J Numer. Anal, pages 871-885, 1982. (66] C. Johnson. Numerical Solutions of Partial Differential Equations by the Finte Element Method. Cambridge University Press, Cambridge, UK, 1987. [67] C. Johnson and U. Navert. An analysis of some finite methods for advection-diffusion equations. In O. Axelsson, L. Frank, and A van Sluis, editors, Analytical and Numerical Approaches to Asymptotic Prob lems in Analysis. North-Holland, Amsterdam, 1981. [68] C. Johnson and A. Szepessy Adaptive finite element methods for con servation laws based on a posteriori error estimates. Comm. Pure Appl Math pages 199-234, 1995. [69] C. Johnson, A. Szepessy, and P. Hansbo. On the convergence of shock capturing streamline diffusion finte element methods for hyperbolic con servation laws. J. Math. Comp., pages 107-129, 1990. (70] K. L. Kipp, L. F. Konikow, and G Z. Hornberger. An implicit disper sive transport algorithm for the U.S. Geological Survey MOC3D solute transport model. Water Resources Investigations Report 98-4234, U.S. Geological Survey, 1998. 126

PAGE 144

[71] L. F. Konikow, D. J. Goode, and G. Z. Hornberger. A three-dimensional method-of-characteristics solute-transport model (MOC3D). Water Resources Investigations Report 96-4267, U.S. Geological Survey, 1996. [72] Erwin Kreyszig. Introductory Functional Analysis with Applications. John Wiley & Sons, Inc. New York, 1978. [73] D .A. Lavis and B.W. Southern. The inverse of a symmetric banded toeplitz matrix. Reports on Mathematical Physics, 39:137-146, 1997 [74] G. Lube and D. Weiss. Stabilized finite element methods for singularly perturbed parabolic problems. Appl. Numer. Math., pages 431-459, 1995. [75] Jan Mandel, Marian Brezina, and Petr Vanek. Energy optimization of algebraic multigrid bases. Computing, 62:205-228, 1999. June 22, 1999. [76] K. W. Morton, A. Priestly, and E. Siili. Stability of the lagrangian galerkin method with nonexact integration. Model. Math. Anal. Numer., pages 625-653, 1988. (77] S. P. Neuman. Adjoint Petrov-Galerkin method with optimal weight and interpolation functions defined on multi-dimensional nested grids In G. Gambolati et al., editor, Computational Methods in Water Resourc e s VII, pages 347-356 Springer, New York 1990. [78] G. F. Pinder and H. H. Cooper A numerical technique for calculating the transient position of the saltwater front. Water Resources Research, pages 875-882, 1970. [79] O. Pironneau. On the transport diffusion algorithm and its application to the Navier-Stokes equations. Numer. Math, pages 309-332, 1982. [80] J. W. Ruge and K. Stiiben. Algebraic multigrid (AMG). In S. F. McCormick, editor, Multigrid Methods, volume 3 of Frontiers in Applied Mathematics, pages 73-130. SIAM, Philadelphia, PA, 1987. [81] T. F. Russell. Eulerian-Lagrangian localized adjoint methods for advection-diffusion problems In D Griffiths and G. Watson, editors, Proceedings of the 13th Dundee Conference on Numerical Analysis, 1989 Pitman Res. Notes Math Ser. 228, pages 206-228. Longman Scientific and Technical, Harlow, UK, 1990. 127

PAGE 145

[82] T. F. Russell and R. V. Trujillo. Eulerian-Lagrangian localized adjoint methods with variable coefficients in multiple dimensions. In G. Gam bolati, editor, Computational Methods in Surface Hydrology, pages 357-363. Springer, Berlin, 1990. [83] T. F. Russell, M. F. Wheeler, and C. Chiang. Large-scale simulation of miscible displacement by mixed and characteristic finite elements. In W. E. Fitzgibbon, editor, Mathematical and Computation Methods in Seismic Exploration and Reservoir Modeling. SIAM, Philadelphia, 1986. [84] M. Stynes and T. F. Russell. An optimal-order estimate for the finite volome Eulerian-Lagrangian localized adjoint method for advection problems. [85] J. E. Vag, H. Wang, and H. K. Dahle. Eulerian-Lagrangian localized ad joint methods for systems of nonlinear advective-diffusive-reactive equa tions. Advances in Water Resources, pages 297-315, 1996. [86] P. Vanek. Acceleration of convergence of a two-level algorithm by smoothing transfer operator. Applications of Mathematics, 37:265-274, 1992. [87] P. Vanek, J. Mandel, and M. Brezina. Algebraic multi grid by smoothed aggregation for second and fourth order elliptic problems. Computing, 56:179-196, 1996. [88] P. Vanek, J. Mandel, and M. Brezina. Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems. Computing, 56:179-196,1996. [89] Petr Vanek, Marian Brezina, and Jan Mandel. Convergence of algebraic multigrid based on smoothed aggregation. 2000. To appear in Numer. Math. [90] Petr Vanek, Marian Brezina, and Radek Tezaur. Two-grid method for linear elasticity on unstructured meshes. SIAM J. Sci Comp., 21(3):900-923, 1999. [91] Petr Vanek, Jan Mandel, and Marian Brezina. Algebraic multigrid on unstructured meshes. Technical Report 34, UCD/CCM, 1994. [92] Petr Vanek, Jan Mandel, and Marian Brezina. Two-level algebraic 128

PAGE 146

multigrid for the Helmholtz problem. In Jan Mandel, Charbel Farhat, and Xiao-Chuan Cai, editors, Domain Decomposition Methods in Sci ence and Engineering, pages 349-356. American Mathematical Society, Providence, RI, 1998. Proceedings of the 10th International Symposium on Domain Decomposition Methods, Boulder, Colorado, 1997. [93] E. Varoglu and W. D. L. Finn. Finite elements incorporating character istics for one-dimensional diffusion-convection equation. J. Compo Phys, pages 371-389, 1980. [94] H Wang. Eulerian-Lagrangian localized adjoint methods: analyses, nu merical implementations and their applications. PhD thesis, University of Wyoming, Laramie, Wyoming, 1992. Ph.D. thesis. [95] H. Wang and R. E.Ewing. Optimal order convergence rates for ELLAM for reactive transport and contamination in groundwater. Numer. Methods Partial Differential Equations, pages 1-31, 1995. [96] H. Wang, R. E. Ewing, and T. F. Russell. Eulerian-Lagrangian localized adjoint methods for convection-diffusion equations and their convergence analysis. IMA J.Numer. Anal., pages 405-459, 1995. [97] H. Wang, R. C. Sharpley, and S. Man. An ELLAM scheme for advection diffusion equations in multi-dimensions. In Aldama et al., editor, Com putational Methods in Water Resources XI, pages 99-106. Computa tional Mechanics Publications, Southampton, England, 1996. [98] J. J. Westerink and D. Shea. Consider higher degree Petrov-Galerkin methods for the solution of the transient convection-diffusion equation. Internat. J. Numer. Methods Engrg., pages 1077-1101, 1989. [99] E. J. Wexler. Analytical solutions for one-, two-, and three-dimensional solute transport in ground-water systems with uniform flow. Techniques of Water-Resources Investigations Book 3, Chapter B7, U.S. Geological Survey, 1992. [100] M. F. Wheeler and C. N. Dawson. An operator splitting method for advection-diffusion-reaction problems. In MAFELAP Proceedings 6, pages 463-482. Academic Press, New York, 1988. [101] D. Yang. A characteristic-mixed method with dynamic finite element 129

PAGE 147

space for convection-dominated diffusion problems. J. Comput. Appl. Math., pages 343-353, 1992. [102] G. Zhou. An adaptive streamline diffusion finite element method for hy perbolic systems in gas dynamics. PhD thesis, University of Heidelberg, Heidelberg, Germany, 1992 Ph.D. thesis. [103] G. Zhou. A local error analysis of the streamline diffusion method for nonstationary convection-diffusion systems. In Mathematical Modeling and NumericalAnalysis, pages 577-603. RAIRO, 1995. 130