E PA/600/2-91/014
April 1991
DEVITRIFICATION IN NONHOMOGENEOUS LABORATORY SCALE AQUIFERS:
4. HYDRAULICS, NITROGEN CHEMISTRY, AND MICROBIOLOGY IN A SINGLE LAYER
by
F. T. Lindstrom, L. Boersma, D. Myrold, and M. Barlaz
Department of Soil Science
Oregon State University
Corvallis, OR 97331
Cooperative Agreement CR 814502-01-1
Project Officer
Thomas E. Short
Processes and System Research Division
Robert S. Kerr Environmental Research Laboratory
Ada, OK 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OK 74820

-------
DISCLAIMER
The information iri this document has been funded wholly or in part
by the United States Environmental Protection Agency under Assistance
Agreement No. CR-814502-01-1 to Oregon State University. It has been
subjected to the Agency's peer and administrative review, and it has
been approved for publication as an EPA document. Mention of trade
names or commercial products does not constitute endorsement or
recommendation for use.
ii

-------
FOREWORD
EPA is charged by congress to protect the Nation's land, air and
water systems. Under a mandate of national environmental laws focused
on air and water quality, solid waste management and control of toxic
substances, pesticides, noise, and radiation, the agency strives to
formulate and implement actions which lead to a compatible balance
between human activities and the ability of natural systems to support
and nurture life.
The Robert S. Kerr Environmental Research Laboratory is the
agency's center of expertise for investigation of the soil and
subsurface environment. Personnel at the laboratory are responsible for
management of research programs to: (a) determine the fate, transport
and transformation rates of pollutants in the soil, the unsaturated zone
and the saturated zones of the subsurface environment; (b) define the
processes to be used in characterizing the soil and subsurface
environment as a receptor of pollutants; (c) develop techniques for
predicting the effects of pollutants on groundwater, soil and indigenous
organisms; and (d) define and demonstrate the applicability and
limitations of using natural processes, indigenous to the soil and
subsurface environment, for the protection of this resource.
Laboratory scale, physical, test aquifers are increasingly being
used for the study of transport and fate processes. It is generally
less expensive to evaluate hypotheses for pollution fate and transport
using laboratory scale aquifers than to work under field conditions.
Furthermore, numerical models of the fate and transport of chemicals in
aquifers are now rapidly coming within the reach of environmental
scientists. These models are even cheaper and faster to operate than
are laboratory systems. The current investigation will analyze the
ability of laboratory aquifers to validate mathematical models. This
report presents a mathematical analysis of two dimensional horizontal
transport and fate processes in aquifers. This analysis led to
development of a computerized model. In further development (later
report), results obtained from the mathematical model will be compared
to the results obtained in a laboratory test aquifer used to study in
situ denitrification as the first test case.
Clinton U. Hall
Director
Robert S. Kerr Environmental
Research Laboratory
iii

-------
ABSTRACT
A two dimensional mathematical model for simulating the transport
and fate of organic chemicals in a laboratory scale, single layer
aquifer is presented. The aquifer can be nonhoraogeneous and anisotropic
with respect to fluid flow properties. The physical model for which
this mathematical model has been developed is assumed to have open inlet
and outlet ends and to be bounded by impermeable walls on all sides.
The mathematical model allows placement of fully penetrating injection
and/or extraction wells anywhere in the flow field. The inlet and
outlet boundaries have user prescribed hydraulic pressure fields. The
steady state hydraulic pressure field is obtained first, by using the
two-dimensional Darcy flow law and the continuity equation.
Separate dynamic transport and fate equations are then set up for
each of four dissolved chemicals, which include a substrate, nutrients,
oxygen, and nitrate. Two equations, modeling the local growth and decay
of two microbial populations, one operating with either oxygen or
nitrogen, the other only with oxygen, are coupled to the transport and
fate equations. The four chemical transport and fate equations are then
solved in terms of user prescribed initial conditions. Boundary condi-
tions are zero flow at the top, bottom, and sidewalls and accounting of
mass at the inlet and exit ports. The model accounts for the major
physical processes of dispersion and advection, and also can account for
linear equilibrium sorption, four first order loss processes, including
irreversible chemical reaction and/or dissolution into the organic
phase, and irreversible binding in the sorbed state. The loss of
substrate, nitrate, nutrient, and oxygen is accounted for by modified
Monod kinetic type rate rules. The chemical may be released internally
by distributed sources that do not perturb the flow field, or from fully
penetrating injection wells. Chemical compound may also enter at the
inlet boundary. Chemical mass balance type inlet and outlet boundary
conditions are used. The solution to the field equation for hydraulic
pressure is approximated by the space centered finite difference method
using the strongly implicit procedure (SIP) with a user specified
heuristic for choosing the iteration parameter. A solution to the
transport and fate equations is approximated with a forward in time
Euler-Lagrange time integrator applied to the chemical transport and
fate semi-discretization. A user manual for the model is available in a
separate document.
iv

-------
CONTENTS
Disclaimer		ii
Foreword		iii
Abstract		iv
List of Figures		vii
1.	Introduction		1
2.	Statement of the Problem		3
1.	Introduction		3
2.	Physical Aquifers at RSKERL		3
3.	Long-term Goals of the Mathematical Modeling Effort	4
4.	Assumptions Underlying Model LT3VSI		5
3.	Nomenclature and Notation used for the Model		8
1.	Fluid Flow Field		8
2.	Chemical Transport and Fate Field		9
3.	Chemical Transformations		11
4.	Fluid Flow Field		15
1.	Statement of the Problem		15
2.	Approximation of the Fluid Flow Equations		17
3.	Boundary Conditions		21
4.	Solution of the Linear System of Equations		22
5.	Velocity Components		23
v

-------
CONTENTS (continued)
5.	Chemical Transport and Fate Model		26
1.	Transport and Fate Equations		26
2.	Initial Data, Boundary Data, and Total Chemical Mass	35
3.	Approximate Solution for the Transport and Fate
Equations		40
h. Boundary Conditions for Approximate Concentration
Field		45
5.	Approximate Microbial Populations		51
6.	Stability Considerations		53
6.	Conclusions for the Model		55
7.	Literature Cited		57
Appendix A		61
vi

-------
FIGURES
Figure 1. Schematic diagram of physical models used for
the studies of aquifer restoration at Robert S.
Kerr Environmental Research Laboratory	4
Figure 2. Schematic diagram of the physical model used for
studies of transport and fate problems	7
Figure 3. Local interior region of lattice points showing
subregion	18
Figure 4. Schematic diagram of the physical aquifer used for
testing the mathematical model, showing mixing
chambers at inlet and outlet ends	37
Figure 5. Details of the use of the "Gaussian pill box"
concept in the mixing chambers, showing flow
orientation, and positions where mathematical
terms were defined	37
Figure 6. Relations of the various subregions of the
overall space-time cylinder and the superimposed
set of lattice points	41
Vll

-------
DENITRIFICATION IN NONHOMOGENEOUS LABORATORY SCALE AQUIFERS:
4. HYDRAULICS, NITROGEN CHEMISTRY, AND MICROBIOLOGY IN A SINGLE LAYER
1. INTRODUCTION
This is the fourth report in a series describing mathematical
models for denitrification processes in nonhomogeneous, laboratory scale
aquifers. This report expands the mathematical model developed in
Volujne 1 (Lindstrom and Boersma, 1990) of the series. The nitrogen
chemistry and microbiological processes that occur in a single layer,
saturated, aquifer have been added to this model. This development is
an intermediate step towards completing a fully three dimensional model
that applies to an aquifer made up of three layers. Such an aquifer is
used at the USEPA Laboratory at Ada, Oklahoma for testing bioremediation
processes for denitrification. Details of the bioremediation processes
being investigated and the justification for this research were de-
scribed in detail in Volume 1 of this series (Lindstrom and Boersma,
1990).
The success of biological denitrification methods depends on avail-
ability of a fundamental understanding of the transport and fate
processes. In addition, knowledge about the important limiting factors,
or limiting system properties must be acquired. Therefore, the immedi-
ate objective was to develop a preliminary mathematical model (Lindstrom
and Boersma, 1990) and associated computer code (Bachelor et al., 1990)
to describe substrate injection into a single layer, laboratory scale
aquifer and to use the model in a sensitivity manner to assess the
magnitude of the physical and biological factors controlling aquifer
denitrification processes and identify those which can be manipulated to
1

-------
enhance the process (Boersma and Lindstrom, 1990). The sensitivity
analysis is reported in the associated documents Volume 2 (Bachelor et
al., 1990) and Volume 3, which describes applications to 2-D non-
homogeneous flow fields (Boersma and Lindstrom, 1990).
The preliminary mathematical model addresses aquifer fluid trans-
port phenomena including perturbations due to the presence of injection
and withdrawal wells. Solute movement includes advection, dispersion,
molecular diffusion, and simple first order chemical and/or biological
reactions in the aquifer. The mathematical modeling effort presented
here is part of an ongoing study of aquifer restoration being conducted
by the U.S. Environmental Protection Agency at the R. S. Kerr Environ-
mental Research Laboratory in Ada, Oklahoma. One of the long-term goals
of this study is development of a mathematical model of aquifer denitri-
fication processes enhanced by stimulation of microbial populations.
The experience gained from developing the preliminary model was used to
develop the two-dimensional model for the simultaneous transport and
fate of nitrate, and oxygen, substrate, e.g. methanol and inorganic
nutrients, in the single layer aquifer, described in this report. Two
independently operating microbial populations are included in the model,
both using modified Monod kinetics. The model, called LT3VSI, is
described in this document.
2

-------
2. STATEMENT OF THE PROBLEM
2.1 Introduction
Laboratory scale, physical models of aquifers are increasingly
being used for the study of aquifer processes. Often it is less
expensive to evaluate hypotheses for restoration procedures using
laboratory scale models than to work under field conditions. Aquifer
restoration methods are not without hazards. Biorestoration processes
may alter the hydraulic properties of the aquifer.
Furthermore, numerical models of the transport and fate of
chemicals in aquifers are now rapidly coming within the reach of
environmental scientists. These models, once validated for the systems
they are designed to simulate, are often cheaper and much faster to
operate in real time response than the physical models. Thus the
combination of mathematical models and laboratory scale models presents
a cost effective and time efficient method for the study of bioremedia-
tion of contaminated aquifers.
2 . 2 Physical Aquifers at RSKF.RL
Two large (4 ft wide, 4 ft high, 16 ft long) physical aquifers
(Figure 1) were constructed at the USEPA Robert S. Kerr Environmental
Research Laboratory in Ada, Oklahoma. Each aquifer contains three
horizontal layers of material, with each layer assumed to be homogeneous
and isotropic with respect to water flow. These systems can be used for
validation of mathematical models that simulate the hydrodynamic pres-
sure distribution (Yates, 1988a,b), for the study of transport and fate
of chemicals, and for evaluation of the growth characteristics of
3

-------
indigenous microbial populations. The physical aquifers are also used
for the study of proposed physical and biological remediation schemes.
No flow walls
Open bp
No Row
prescribed
total hydraulic
hwd
partial ly
.penetrating
well (extraction)
No flow

prescribed
total hydraulic
head .
J	
Aquidude
(bottom slab)
— x - L
No now
No now bottom
— x - 0
Figure 1. Schematic diagram of physical models used for the studies of
aquifer restoration at Robert S. Kerr Environmental Research
Laboratory
2.3 Long-term Goals of the Mathematical Modeling Effort
The goal of the present mathematical modeling effort is to describe
the fate and transport of contaminants in the physical models. This
work includes two-dimensional mathematical modeling of the steady state
hydraulics and simultaneous transport and fate of the dissolved oxygen,
nutrients (such as phosphorus), a carbon based substrate (for example
methanol), and dissolved nitrate. Also included are two microbial
populations which change in space and time. The mathematical model will
be used to study scenarios for biorestoration of aquifers contaminated
with nitrates.
k

-------
The procedures being followed to achieve these goals include
several steps. The first step was development of a preliminary model of
the transport and fate of chemical compounds, with constant first-order
loss processes and linear equilibrium sorption assumptions. This model
is for two space dimensions and simulates only the "aquifer" slab of the
physical model. The second step was to use the preliminary, two-
dimensional model, called LT2VSI, for preliminary numerical studies of
several scenarios for injection and/or extraction well placement. The
third step was to expand LT2VSI by including four chemical compounds and
two microbial populations. This model, referred to as LT3VSI, is the
subject of this report.
2.h Assumptions Underlying Model LT3VSI
The RSKERL, physical models have been constructed in such a way that
homogeneous and isotropic soil slabs were obtained. They have imperme-
able (no flow) side walls, an open top, partially open ends, and an
impermeable lower or bottom boundary. The assumptions regarding the
walls and the bottom made in the preliminary mathematical model LT2VSI
(Figure 2) reflect these conditions with the exception that the top
boundary was assumed to be a "no flow" boundary. The hydraulic head
distributions at the completely open inlet and exit boundaries are
prescribed in model LT2VSI and fully penetrating injection and extrac-
tion wells may be present. Model LT2VSI was developed for a single-
layer of soil representing the aquifer part of the three soil layers
making up the RSKERL aquifers.
5

-------
This report describes an extension of the preliminary model which
is referred to as LT3VSI. The hydraulics in LT3VSI are the same as in
LT2VSI. The important features of LT3VSI are:
1)	Two-dimensional, horizontal, steady state, fluid flow field
defined by a hydraulic head field which depends on appropriate
Dirichlet and Neumann boundary conditions and characterizing
the spatial dependency of the longitudinal and transverse
components of the hydraulic conductivity tensor at saturation,
2)	Simultaneous two-dimensional transport and fate of four
dissolved chemicals, i.e.: oxygen, substrate, nutrient, and
nitrate in the nonhomogeneous aquifer. The distribution of
chemicals is affected by
i)	advection and dispersion in both the longitudinal and
transverse directions,
ii)	linear equilibrium adsorption/desorption processes on each
of the porous medium fractions,
iii)	three different first-order loss processes
a)	chemical reaction with other soil components in the
free phase,
b)	other irreversible processes in the free phase,
c)	chemical reaction in the sorbed phase,
iv) modified Monod microbial degradation kinetics of the
substrate system including oxygen, nutrients, and nitrate
in addition to two microbial populations,
v) the presence of zero order sources of the chemicals,
6

-------
vi) appropriate Dirichlet and no flux Neumann boundary
conditions with a provision for nonzero initial
distribution of the chemicals,
vii) presence of fully penetrating injection and/or extraction
wells.
It is assumed that most, if not all, of the chemical and biological
process coefficients are at least once continuously differentiable
functions of the transverse coordinate x and the longitudinal coordinate
y (Figure 2) .
z = L,
prescribe
hydraulic
head
H(x,0) (meters) Jr- —
/ '
/ x « Lx
.1^	
z = 0
/
x-o
o

Fully penetrating
injection or extraction
well
i i
o
OUTLET
prescribed
total
hydraulic
— head - — — -
H(xtLy) (meters)
/
y - Ly
Figure 2. Schematic diagram of the physical model used for studies of
transport and fate problems
7

-------
D
H
H
3. NOMENCLATURE AND NOTATION USED FOR THE MODEL
3.1. Fluid Flow Field
A real nonsymmetric positive definite
(Nx-1) • (Ny-1) by (Nx-1) • (Ny-1)
array used for computing H
vector of boundary data in the hydraulic
head field
(x,y) in Rz| (0 < x < L,) X (0 < y < Ly) }
hydraulic head field
vector of approximate hydraulic heads
m"1 day 1
day"1
dimensionless
m water
m water
Hir.
^out
K^yy
K
h
K
q
q*
qy
approximate hydraulic head field at
location (xj, yj)
hydraulic head field at y - 0
hydraulic head field at y = 1^
transverse component of saturated
hydraulic conductivity tensor
longitudinal component of saturated
hydraulic conductivity tensor
transverse dimension (width) of aquifer
longitudinal dimension (length) of aquifer
vertical dimension of aquifer
number of subintervals in interval [0, L^] .
There are Nx-1 internal nodes on the
transverse coordinate and Nx+1 total
nodes along this same coordinate.
number of subintervals in interval [0, Ly] .
There are Ny-1 internal nodes on the
longitudinal coordinate and Ny+1 total
nodes along this same coordinate.
Darcy velocity vector
x-component of Darcy velocity vector
y-component of Darcy velocity vector
m water
m water
m water
m day"1
m day'1
m
m day"1
m day"1
m day"1
8

-------
Q	mass density rate function for sources	kg m"3 day"1
and/or sinks; these rules must be small
relative to the inlet flow rates
Qinj	mass density rate of injection well	kg m"3 day"1
Qout	mass density rate of extraction well	kg nT3 day~2
U	average seepage velocity vector	m day"1
Ux	x-coraponent of average seepage velocity	m day-1
vector
-l
Uy	y-component of average seepage velocity	m day
vector
V^v	representative elementary volume	m3
x	transverse spatial component	m
y	longitudinal spatial component	m
QAj	ijth correction factor used for correcting
the area estimate for adjacent injection
and/or pumping wells
e	volume porosity of porous medium	m3 m"3
density of water	kg m"3
V	"del" or nabla vector operator	nf1
V2	Laplace operator	m"2
3.2. Chemical Transport and Fate Field
The superscript symbol * in the equations denotes one of the four
chemicals of interest in the transport and fate fields, as follows: Nu
for nutrient; S for substrate; 0 for oxygen; and Ni for nitrate.
C*	distribution of chemical concentration	kg m~3
in the free phase
*	i
Cin	chemical concentration in the inlet end	kg m
mixing tank
C.,ut	chemical concentration in the outlet end	kg m~3
mixing tank
*	-1
C30ur	chemical concentration of the injection	kg 1
wells
9

-------
*
o
*
cin
cout
XX
*
yy
*
xy
in
^*yin

out

M* (t)
out
so
*
V,
Qdispx
'dispy
chemical concentration in feed stream
average dispersion coefficient in the
inlet tank
average dispersion coefficient in the
outlet tank
molecular diffusion coefficient of
compound in free solution
x-component of dispersion
y-component of dispersion
cross component of dispersion
length of the inlet end mixing tank
length of the outlet end mixing tank
distance from the plane y-0 to the
center of the inlet end mixing tank
distance from the plane y=Ly to the
center of the outlet end mixing tank
mass of compound * in the inlet tank
at time t days
mass of compound * in the outlet tank
at time t days
source rate density from buried sources
retention parameter
elapsed time in days from commencing
injection and/or pumping
magnitude of x-component of the seepage
velocity
magnitude of y-component of the seepage
velocity
average flow velocity in inlet tank
average flow velocity in outlet tank
x-component of dispersivity
y-component of dispersivity
10
kg m"3
m2 day'1
mz day"1
m2 day"1
m2 day"1
m2 day"1
m2 day"1
m
m
m
m
kg
kg
kg m"3 day"1
dimensionless
days
m day"1
ra day"1
m day"*
m day"1
m
m

-------
atort	tortuosity factor	diraensionless
At	time increment t^+j = t^ + At	days
3.3 Chemical Transformations
*>t	-
Korg	linear equilibrium distribution constant	m kg~'org
for organics expressed in volume per
unit mass of organic matter
*	i	n
K„e	linear equilibrium distribution constant	m" kg
for mildly sorbing surfaces, expressed
as volume per unit mass of mildly sorbing
particles
¦Jr
Ks	linear equilibrium distribution constant	mJ kg 3
for strongly sorbing surfaces, expressed
as volume per unit mass of strongly sorbing
particles
Ks0	one half maximum saturation constant for	kg m~3
population 1 utilizing the substrate
Kq	one half maximum saturation constant for kg ra~3
population 1 utilizing oxygen
Koru	one half maximum saturation constant for	kg m~3
population 1 utilizing nutrients
2	,
Kso	one half maximum saturation constant for	kg m
population 2 utilizing the substrate
2
K0	one half maximum saturation constant for	kg m~3
population 2 utilizing oxygen
2
Konu	one half maximum saturation constant for	kg m~3
population 2 utilizing nutrients
1	,
Ksni	one half maximum nitrate based saturation kg m J
constant for population 1 utilizing
substrate
Kjji	one half maximum nitrate based saturation kg m~3
constant for population 1 utilizing
nitrate
K^lnu	one half maximum nitrate based saturation kg m~3
constant for population 1 utilizing
nutrients
K^ni	one half maximum nitrate based inhibition kg m~3
constant for population 1 utilizing nitrate
instead of oxygen as its electron acceptor
11

-------
k}oii1	one half saturation coefficient for	kg m"3
maintenance in population 1
2
Kjot,	one half saturation coefficient for	kg ro~3
maintenance in population 2
N:(t)	population 1 cell mass in the unit	kg kg"1
pore volume at elapse time t days in
microbial mass per unit soil mass
N2(t)	population 2 cell mass in the unit	kg kg"1
pore volume at elapse time t days in
microbial mass per unit soil mass
rsor.u	oxygen utilization rate by microbial	kg kg"1 day"1
population 1 in substrate mass per
unit cell mass per day
rEr.ir.u	nitrate utilization rate by microbial	kg kg"* day"1
population 1 in substrate mass per
unit cell mass per day
2
rsonu	oxygen utilization rate by microbial	kg kg"1 day"1
population 2 in substrate mass per
unit cell mass per day
S	chemical concentration distribution in	kg m"J
the sorbed phase for linear equilibrium
Y^0	yield coefficient for microbial	kg kg"1
population 1 utilizing oxygen in cell
mass per unit mass of substrate
2
Yso	yield coefficient for microbial	kg kg"1
population 2 utilizing oxygen in cell
mass per unit mass of substrate
Yjni	yield coefficient for microbial	kg kg-1
population 1 utilizing nitrate in cell
mass per unit mass of substrate
M0/M.	mass fraction of organic matter in	kg kg"1
porous medium
M^/Mt mass fraction of weakly sorbing particles kg kg"1
in porous medium
Ms/Mt	mass fraction of strongly sorbing	kg kg-1
particles in porous medium
12

-------
a	oxygen maintenance rate coefficient based kg kg"1
on population 1 using oxygen under
aerobic conditions in mass of oxygen per
unit mass of cells
2
a	oxygen maintenance rate coefficient based kg kg *
on population 2 using oxygen under
aerobic conditions in mass of oxygen per
unit mass of cells
nitrate maintenance rate coefficient based kg kg"1
on population 1 using nitrate under
aerobic conditions in mass of nitrate per
unit mass of cells
51	constant for death rate of population 1	day"1
52	constant for death rate of population 2	day"1
coefficient for nitrate use based on	kg kg"1
1
a .
ni
population 1 and anaerobic conditions
7^	coefficient for oxygen use based on	kg kg"1
population 1 using oxygen under aerobic
conditions in mass of oxygen per unit
mass of substrate
2
7	coefficient for oxygen use based on	kg kg"1
population 2 using oxygen under aerobic
conditions in mass of oxygen per unit
mass of substrate
A	first-order free phase loss rate	day"'
constant for irreversible loss
processes, e.g. chemical reactions
S*	*
A.	first-order sorbed phase loss rate	day
irr	r	J
constant for irreversible processes
*
A	overall first-order loss processes rate	day""
coefficient
A (x,y,t) loss rate function of compound * due	kg* m"3 day"1
to microbial processes
nutrient use coefficient based on	kg* m~3 day"1
population 1 and oxygen conditions in
nutrient mass per unit mass of substrate
2
nutrient use coefficient based on	kg* nf3 day 1
population 2 and oxygen conditions in
mass of nutrient per unit mass of substrate
13

-------
pk	average bulk density of porous medium	kg m 3
pws	average particle density of weakly sorbing kg m*3
fraction
ps	average particle density of strongly	kg m"3
sorbing fraction
¦'org
ni
average particle density of organic	kg nf3
matter fraction
nutrient use coefficient based on	kg m"3
population 1 using nitrate under anaerobic
conditions in nutrient mass per unit mass
of substrate
^	maximum specific growth rate for	day"1
heterotrophic population 1
2
H	maximum specific growth rate for	day"1
heterotrophic population 2
.	maximum specific nitrate based on growth day"1
rate for heterotrophic population 1
14

-------
4. FLUID FLOW FIELD
4.1. Statement of the Problem
Figure 2 shows a schematic diagram of the single layer aquifer of
the physical model, with water flowing into and out from the soil via
open ends. Water can also be extracted or injected through fully
penetrating injection and/or extraction wells. The transverse coordi-
nate is x (m) and the longitudinal coordinate is y (m). The third
coordinate L^ (m) represents aquifer thickness. The porous medium may
be nonhomogeneous and anisotropic with respect to fluid flow properties
and has impervious walls on all sides. These conditions allow
evaluation of two-dimensional transport and fate.
The assumptions concerning the fluid flow field are:
1)	The fluid flow field operates at steady state conditions at
all times.
2)	Any fluid flow perturbations introduced at the flow boundaries
propagate extremely rapidly throughout the flow field, so that
a new steady state is achieved instantly. Under these condi-
tions the fluid storativity term in the fluid flow model may
be neglected as shown for aquifers of the size considered here
(Bodvarsson, 1984).
3)	The aquifer material can be nonhomogeneous as well as aniso-
tropic, with the principal components Ksxx and Ksyy of the
saturated hydraulic conductivity tensor assumed to be once
continuously differentiable over the interior of the flow
domain.
4)	Dirichlet boundary conditions hold at both the inlet and
outlet ends. H1ti (m) at the inlet end and Hout (m) at the
15

-------
outlet end are specified and assumed to extend across, 0 < x <
L*-
5)	No-flow Neumann or flux type boundary conditions are specified
along the walls, i.e. the x = 0 and x => L* planes.
6)	The Darcy velocity field components of the flow vector q(m
day"1) are defined by
q = (f U , e U ) = (- K	, -K ~~) ,	(A.1.1)
n	x y	sxx Sx syy ay
where Ksxx and Ksyy are the transverse and longitudinal compo-
nents (m day"1) of the saturated hydraulic conductivity in the
aquifer, e is porosity, and H is the hydraulic head (m). The
coordinate system has been chosen to coincide with the two
principal directions of flow (Bear and Verruijt, 1987).
Combining equation (4.1.1) with the continuity equation for steady
state conditions in a representative elementary volume (REV) in the
aquifer, obtains
v .  ,	(4.1.2)
where pw is the density of water (kg m~3) is assumed constant, Qir_j is
the fluid mass injection rate (kg m"3 day"1) and Qout is the fluid mass
extraction rate (kg m"3 day"1).
Substitution for (e U) in equation (4.1.2) leads to the well known
Poisson problem in potential field theory:
'w
d(-K 6s) 3 (-K f*)
	sxx ox +	svv dv
3x	dy
ra(£ U ) 3(e U )]
Pw ^ 3x + 3y '
+ (^inj " W'	(413)
16

-------
where e, KSJQC, and Ksyy are assumed to be once continuously differentia-
ble functions of x and y inside the flow domain.
Boundary conditions to be satisfied at all times, are written
according to assumptions 4 and 5 as follows:
i) along y - 0 (inlet end): H(x,0) - Hln(x), 0 < x < L,,	(4.1.4)
ii) along y = Ly (outlet end): H(x,Ly) = Hout(x), 0 < x < 1^, (4.1.5)
iii) along x - 0:	fx ~ ^
iv) along x = L :	= 0 .	(4.1.7)
x	dx
Since Ksxx and Ksyy are usually not constant and may vary by several
orders of magnitude over the interior of the flow domain, analytic
solutions of equation (4.1.3), subject to boundary conditions (4.1.4)
through (4,1,7), are not obtainable. For isotropic and homogeneous
porous media, Lindstrom et al. (1989) give an analytical representation
in terms of a double sura infinite series of trigonometric functions,
i.e. an eigenfunction solution procedure, However, even in the very
special isotropic case of KSXI = Ksyy - Ks, where Ks is a constant, the
solutions are usually very slow to converge. Thousands if not millions
of terms are necessary to achieve the required nuraber of significant
digits in each velocity component. Thus, the hydraulic head field on D,
the interior of the flow domain, is usually approximated by means of
finite difference or finite element methods.
4.2. Approximation of the Fluid Flow Equations
Figure 3 shows a local region of D on which a lattice of nodal
points has been superimposed. D is the set of all (x,y)'s in R2 such
17

-------
that { [0 < x < L*] X 0 < y < Ly] ) , with nonhoraogeneous nodal spacing
allowed in both coordinates.
Ax
Ax
*1-1
*i
Figure 3. Local interior region of lattice points showing subregion
Bij-
We choose the integrated steady state hydraulic head field equation
method developed by Varga (1962, pp. 181-191). The details of the
development were given by Lindstrom and Boersraa (1990). The following
system of linear finite difference equations is obtained.
For i = 1,2,... Nx-1, j = 1,2,... Ny-1
18

-------
a.S.t.2. . H, , . +	adtl. . H. . .. + (adt2x. . + adt2y. .) H
i,j i-l.J i.j i,J-l	i.J	i.J i.j
+ adt3 . H	+ aut2. . H . = fln] ——1 . ..
i.J i.J+1	i.J i+l.J I P„ J1-J
(4.2.1)
where H^j is the approximate hydraulic head evaluated at nodal point
(xif yj) . The coefficients given in equation (4.2.1) are defined as
follows:
fK
sxx.
i-1/2,j
ait2.
^ Axi-i J
i.j ¦	t + flx.j •
(^•2.2)
syy;
i.1-1/2
adtl

i.j

(4.2.3)
Ksyyi,j+l/2 + Syyi.j-1/2
adt2y. . =
i.J
iyi

ilu	±
i + ayi
(4.2.4)
K
K
sxx.
t+1/2,1 ^ 1-1/2,1
Ax.
adt2x. .=
!.J
Ax. -
l-l
[ftVi + ax<]
(4.2.5)
adt3.
K
Syyi,j+l/2
Ay,
i,j "fAyi-i + Ayi
l	9
(4.2.6)
19

-------
KsXXi+l/2,j
Ax.
aut2. . - -7-	1—;—7 .	(4. 2. 7)
i.J fA" ' A" 1
fAx. , + Ax.1
—1
Equations (4.2.2) through (4.2.7) show that the components of the
hydraulic conductivity tensor are evaluated at mid-nodal points. The
hydraulic conductivity field is not currently available at these points.
However, the values of KEXX and Ksyy are known at all the nodal points of
the domain D. Because of the C1,1(E) smoothness of Ksxx and Ksyy and the
low order of the finite difference approximations to the first
derivatives themselves, the arithmetic mean between adjacent nodal
points is used as the approximation to the mid-nodal point values of
^SXX	Kgyy
K	+ K
sxx. , . sxx. .
K - 	^^ ,	(4.2.8)
i+1/2,j
K	+ K
sxx.	. sxx.
K - 	r	,	(4.2.9)
i-1/2,j
K	+ K
Syyi,j+l/2 =
K	+ K
Syyi,j-l/2 =
syyi i+i ^ ^
K	- 	 ,J ,	^ ,	(4.2.10)
syyt 1 syy
K	= 	^_	(4.2.11)
20

-------
4.3. Boundary Conditions
Equations (A.1.4) and (4.1.5) are Dirichlet conditions, because the
total hydraulic field along the planes y = 0 and y - Ly is specified.
Hence,
H(x,0) = H. (x.) = H(x.,0), 0 < x < L ,
in i	i	x
(4.3.1)
Similarly at y = Ly define
H. kT = H (x. ) , 0 < x < L .
l,Ny out l	x
(4.3.2)
Equations (4.1.6) and (4.1.7) are Neumann conditions, since the flux in
the direction normal to the surface is specified. In this case, imper-
meable walls along the planes x = 0 and x = L* exist, so that a zero or
no flux boundary condition obtains, i.e. along x — 0, where 3H/3x = 0,
approximate for j - 1,2, Ny.j,
ox
- 0 -
x=0
Ax_ -	Ax	Ax,
(2 + —') Hn . + (—*• + 2 + —A) H. .
Ax1 0,j	Axt	Ax2 1,J
^1 ~
Ax H2,j
- / (Axx + Ax2) ,
(4.3.3)
a formula obtained by differentiating the Lagrange degree 2 interpolat-
ing polynomial,
V>" Vj
H. . - H. .
-Lj	U
Ax,
(X - x0)
2.i ' "1.1 1.1 ' 0.1
Ax,
Ax,
Ax^ + Ax^
(x -
Xq) (x - xQ - Axx)
(4.3.4)
21

-------
setting x = xc (=0) with P2'(x0) = 0, and simplifying the result. Since
Ax! and Ax2 are positive, possibly nonuniform nodal spacings, the defi-
nition of H0 j is obtained as a linear combination of Hjj and H2 j :
0, j
Ax^ Ax"
2 + —- + —-
Ax^ Ax^
hi.j •
r x
AX1 ~
7— H-
Ax2 2¦J

Ax "
2 + A
Ax^
k. >

It can be shovm in like manner that at x - L where ~
x	ox
- 0,
x=L
(4.3.5)
H.
Nx, j
Nx
+
\. ixSx-ll
hk. n .
Nx-1,j
Ax . nNx-2,j
Nx-1 J
2 + +
AxM . Ax.,
Nx-1 Nx

AX M 1
2 + n*'1
AXNx
V J

(4.3.6)
for j - 1,2,... Ny - 1.
A standard Taylor series based error analysis leads to the errors:
at x = 0, 01(Ax1 Ax2) ,	(4.3.7)
which if Axj — Ax2 - Ax is 0-(Ax2) and
at x = Lx, 02(AxNx Ax^-i) ,	(4.3.8)
which if AxNx = Ax^ = Ax is 02(Ax2) .
It is good practice to choose Ax1( Ax2, AxI)x_1, and AxNx equal to about
1/10 the intradomain node spacing of the lattice.
4.4. Solution of the Linear System of Equations
A vertical or horizontal ordering of the nodes may be chosen in
equation (4.2,17). In either case, when using the boundary data, the
system may be written in matrix form
A H = b ,	(4.4.1)
22

-------
where A is a real, nonsymmetric, Nx-1 • Ny-1 by Nx-1 • Ny-1, 5-band array
which is irreducible, weakly diagonally dominant, and positive definite
array, an M-array as defined by Varga (1962, pp. 186-188), and b is the
vector of boundary data and source/sink data. H then exists uniquely as
H = A'1 b .	(4.4.2)
Because storage requirements for a standard direct LU decomposition
method are high, a rapidly converging robust iterative method was chosen
for finding H. The method chosen is the Strongly Implicit Procedure
(SIP) described by Stone (1968) . A heuristic for choosing the
"cancellation of terms" parameter is included in the SIP subroutine
(Bachelor et al., 1990).
4.5. Velocity Components
Once all components of H are known, the x and y components of the
Darcy fluid velocity field, and also the effective pore velocities, can
be estimated using equation (4.1.1). However, since the hydraulic head
field H(x,y) is known only approximately and only on a finite set of
grid points, the two velocity components must be numerically estimated
at each (Xi ,y_j) .
Two 3-point Lagrange polynomials, one in x and the other in y,
similar to those stated by equation (4.3.4), are differentiated at the
intranode points. Differentiation is done first along the x-coordinate
to provide an estimate of oH/3x|j j and then along the y-coordinate to
provide an estimate of di\/5y\i 3. For the nodes inside D and adjacent
to the impervious walls as well as those nodes inside D but adjacent to
23

-------
the inlet and outlet boundaries, velocity estimates are obtained using
the formulas:
at y = 0, for i = 1,2,... ,
fly
y-°
Hi.i • Hi.o
(4.5.1)
5H
dx
Ax.
H.
y=0
Ax. , (Ax. . + Ax.) i-1,0
l-l l-l	l
+ (Ax. " Ax.'' "i,0
l-l	l
AXi-l
H.
Ax.(Ax. , + Ax.) i+1,0 '
l l-l	l
(4.5.2)
at y - Ly, for i - 1,2,... ,
£H
dy
H- V' " H- » 1
l.Nv i.Nv-1
Ay.
y=L,
Ny
(4.5.3)
1H
dx
Ax.
H.
y=K
Ax. ,(Ax. , + Ax.) i-l,Ny
l-l l-l	l	J
+ (Ax.	Ax.^ "i,Ny
l-l	l	J
Axi-1
H.
Ax^(Ax^ ^ + Ax^) i+l,Ny '
(4.5.4)
at x - 0, for j - 1,2,... Ny - 1
5H
dx
— 0, according to equation (4.1.6)
x=0
(4.5.5)
24

-------
3H
9y
x=0
Ayi	~
Ayj _x(Ay _x + Ay^) H0,j-1
+ (Ayj-l ~ Ayj) "°,J
^1.
A>'j (Ayj -i + A>'j ^ °'j+1 '
(4.5.6)
at x = L,, for j - 1,2,... Ny.!
dx
0, according to equation (4.1.7)
x=L
(4.5.7)
£H
ay

x=L
Ayj"l(Ayj-l + A/j) ^'j'1
+ ( 1	H
Ayj-1 " Ayj NX,j
"ii
Ay^(Ay ^ + Ay^.) Nx,j + 1
(4.5.8)
For the remainder of the internal nodes write:
£H
<9x
Ax,
H.

Ax. ,(Ax. . + Ax.) i-l,i
l-l l-l	l	J
-	Ri,j
i-i	i	J
Ax. ,
l-l
Ax^(Ax^ ^ + Ax^) i+l,j'
(4.5.9)
and
25

-------
an
ay
+ Ay.) Hi.j-1
+
J •*•
^yj(Ayj.L + ay ) i,j+i
(4.5.10)
Formulas (4.5.1) and (4.5.3) which approximate oH/3y along the open
portions at y - 0 and y = Ly are 0(Ay) and 0(AyNy) accurate. The
remaining approximations are 0(AXi Ax^) and 0(Ay.j Ayj_x) accurate.
These results are easily confirmed by means of standard Taylor series
analysis. Close to the walls, at the inlet and exit ports, and near the
injection and extraction wells, the nodal spacing must be chosen to be
very small for accurate approximations. In preliminary work it has been
found experimentally that to obtain Ux and Uy to three or four signifi-
cant digits, H(x,y) must be computed to about eight significant digits.
This is so because of the loss of significant digits which occurs with
numerical differentiation of tabular data.
This section describes the basic equations for chemical transport
and fate in the single layer aquifer. The assumptions, which form the
basis for the transport and fate model and which hold for each one of
the four chemical compounds, i.e. substrate, S; nutrient, Nu; oxygen, 0;
and Nitrate, Ni, are now listed.
1) Mass transport is via advection (convection) and dispersion.
5. CHEMICAL TRANSPORT AND FATE MODEL
5.1. Transport and Fate Equations
26

-------
The x and y dispersion components are linearly dependent upon the
moduli of the velocity components of the flow field for two-
dimensional flow in an isotropic and nonhomogeneous aquifer (Bear
and Verruijt, 1987, p. 161-164)
2 *  °*disPx) as	as the
tortuosity (atort,) are assumed to be once continuously differentiable
functions of x and y inside the flow domain D.
27

-------
The porous medium can be partitioned into three distinct fractions
as follows:
i) sorbing particles (clay minerals and small silt particles),
ii) weakly sorbing particles (large silt and sand particles) , and
iii) strongly sorbing organic matter,
with the following linear equilibrium sorption rule assumed to hold
for the porous medium
S* ' = (M /M p K* + M /M p K*
ws t ws ws s t s s
+ M /M P K* )C*	(5.1.5)
o t org org
The assumed linear equilibrium sorption rule of equation (5.1.5)
gives rise to the set of retention coefficients (Lindstrom et al.,
1989)
R* = 1 " 6 ((M /M ) p	+ (M /M ) p K*
e	ws t ws ws	s t s s
+ (M /M ) p K* ) .	(5.1.6)
o t org org
The three porous medium fractions (M.rfS/Mt, MsMt, Me/Mt) are continu-
ous functions of x and y in D.
Chemicals can be introduced into the aquifer with the feed stream at
the inlet end (y = 0) or from constantly emitting sources in the
aquifer. Fluids added by these methods must have a low volumetric
concentration and the flow rates must be low enough so that the
previously established fluid flow field is not disturbed. It is
assumed that density gradients, density stratification, or local
changes in the transport and/or fate properties of the porous medium
do not occur in time.
28

-------
5)	Water containing chemicals can be introduced via fully penetrating
injection wells or extracted from similar wells by pumping.
6)	Loss of chemical can occur via first order irreversible loss proc-
esses such as chemical transformations and precipitation in both the
free and sorbed phases in addition to loss via microbial degrada-
tion. The overall first order irreversible process law for each
chemical is characterized by the constant
* * * *
A - A. + R. X
irr irr
(5.1.7)
All the first order loss process coefficients are assumed to be
continuous functions of x and y inside B.
7) Microbiological processes are modeled using process laws described
by Molz et al. (1986) and Widdowson et al. (1988), who constructed
biodegradation models using the carbon assimilation and oxidation
assumptions of Herbert (1958). The model developed here includes
two microbial populations, namely N^t) (number of cells per kg of
dry soil), which utilizes substrate under both aerobic and anaerobic
conditions, and N2(t) (number of cells per kg of dry soil), which
utilizes substrate only under aerobic conditions. The utilization
rate laws adapted from Widdowson et al. (1988) are
sonu
sonu
1
^o
'
cs
1

r c° i

'
cnu
>
J
J

.¦S



,0

.,1

nu
V
K
+
c

K
+
i;

K
+
i;
so
^ so


k 0


v onu

?







f




cs



c°



cnu

Y2
so
K2
^ SO
+
cs.

K2
^ 0
+
c°

K2 +
onu
cnuJ
(5.1.8)
(5.1.9)
29

-------
where the units of r are (kg kg"1 day"1) . The growth of the two
microbial populations adapted from Widdowson et al. (1988) are
°N1	ll	111
TT1 - [ (i r - SA) + Y1 . r . ] N" , K. (0) - Nin(>0) , (5.1.11)
3t	so sonu	sni sninu 1	1	10
3^2 r ,„2 2	r2
d
^ [(Yso rsonu ' 5 IV N2(0) -N20(£0) "	(5"1'12)
The two maintenance rate coefficients 61 and S2 at time t - 0 are
defined by
61 = (Y1 r1 + Y1 . r1 . ) ,	(5.1.13)
so sonu sni sninu
t-0
2 2 2
6 - (Y r ) .	(5.1.14)
so sonu „
t-0
Consideration of the balance of chemical mass in the REV leads to
the coupled system of four nonlinear, two-dimensional, transport and
fate equations.
Nutrient
a r-nu	f-	QpHU
_nu. o C .. lo _nu dC
Vo™ «(1 + R )	D,	
REV	9t	REV lox	xx dx.
, &— / ^nu dCRU 3 . _nu dCnU, d . _nu dCnu.
ax xy dy * Sy <( yx 3x ' + dy (< yy 3y 1
- f- (« U cnu) - f- (« U Cnu)
dx x	dy y
" Vd™ MCtf1 rl + O1. r1 . ) N (x,y, t)
REV b o sonu m sninu 1 J
30

-------
12 2 __ *	i. __ . ~ _nu	j>.riu
+ i> r N„ (x,y, t)) + V IQ. . C + Q
o sonu 2 J	REV inj sour so
Q
_2U£ c„u ( Anu cnu
Substrate
VREV e(1 + R ) ft" " VREV fe (t °xx fx_)
d . _s 3CS. 9 . s <3CS. 3 . s 3CS.
+ — (« D a—) + ~ (e D t—) + T" (c D -—)
5x xy 3y 5y yx 3x dy yy 3y
d	 . Tt „s. d	 . T.
- T~ (e U C ) - ~ (e L C)
ox x	5y y
}
Vnpn»	+ ^ ) N.(x,y,t)
REV b sonu sninu 1
+ r2 Nn(x,y,t)) + V^JQ. . CS + QS
sonu 2	REV inj sour so
Q
cs - , AS CS)
'w
Oxygen
VREV 'd + R°> ff " VREV {fc 
-------
, 2 2
+ (7 r + a
o sonu
o 52 ^77' »2u.y.«}
som
+ VREVlQinj Csour + Qso
out -O	.0 _o.
	k C - e A C ]
(5.1.17)
Nitrate
_ni
rT	„ni. dC __
REV + * dt = REV I5x ^ "xx 9x
{:
L (f „nl
- f « D"1	+ f- (« D"1	+ f « Dni
dx	xy ay	dy	yx ox	a y yy 0 y
f" 
-------
^ + d* + u* + <^-> c* + A*
ot	X ax	y ay	1 + R	1 + R*
= (e D* jp~) + J" (e D* *p~)
^3x xx 5x dx xy 5y
+ J" (e D* ^-) + f— (e D* ^")| / [ (c (1 + R*))]
dy yx ox dy yy 3y J
* * *
Q + Q. .(C - C ,
so inj sour ~p~)
+ 			,	(5.1.19)
c(l + R )
where
U* = X and U* - y . ,	(5.1.20)
x i + r y 1 + R
and A*(x y t) is the loss rate function due to microbial processes
operating at point (x,y,t) in the space-time cylinder D* = D U [0,T].
The A*(x>y>t) for each chemical is:
nutrient
»nu,	.	,,.11	.11..,,	.
A (x.y.t) - Kh l(0orsonu + «ni rsnlnu> N^x.y.t)
+ rsonu V'"'01 ' < '	(5.1.21)
substrate
AS(x,y,t) ¦= p. ((r1 +r* . ) N (x,y,t) + r2 N (x ,y, t) )/e , (5.1.22)
^	b sonu sninu 1 J	sonu 2 J
33

-------
oxygen
A (x,y,t)
.{
11 ^ 1 r1
7 r	+ a S
o sonu o
K1 + C01
sora
N1(x,y,t)
[2 2	2 2
-y r	+ a 6
o sonu o
,,2 -o
K + C
som
N2(x,y,t)j-
/ e
(5.1.23)
nitrate
xni/	L 1 1
A (x,y, t) = p, A (t) . r
b [ ni sninu
i A
+ a . 6
ni
,ni
k\ + cni
ni
K
oni
,,1
K . + C
oni
} N^x.y,
t) / e .
(5.1.24)
Since the "effective" material or total derivative
it	~	•&	+
d_C_ „ flC_ * ££_ * £C_
dt dt x 3x y dy
d C
dt
is
substitution in equation (5.1.19) obtains
dC
dt

* " * A*fx.v.tl
*
1 + R
C +
f— (t D* ^-) +2~ (c D* Y~)
ay	xx 3x 3x	xy dy
*
* ec
+ (£ D.... 7T") + t: (« d;;„ / [«a + r*)

-------
functions do not exist for the system of equations implied by
equation (5.1.25). Therefore an approximate solution must be obtained
using numerical procedures. The procedure used here is a type of finite
difference Euler-Lagrange procedure (Cheng et al., 1984; Baptista et
al., 1984; Baptista, 1986; Baptista, 1987), which is a modification of
the method of characteristics. Equation (5.1.25) was written in the
mixed Euler-Lagrange form in anticipation of choosing this method of
approximation.
5.2. Initial Data. Boundary Data, and Total Chemical Mass
Completion of the conditions on C* over the space-time cylinder D*
requires C" to be prescribed at time t - 0, i.e.,
Equation (5.2.2) holds for each of the four chemicals.
Recall from Figure 2 that at each of the two planes y - 0 and y -
Ly the total hydraulic head is prescribed rather than the fluid fluxes.
The fluid velocity field can vary across these planes, depending upon
the presence or absence of injection and/or extraction wells in the
interior. The chemical mass in the system is also desired. A complete
accounting of all the chemical mass which enters and exits the aquifer
via the inlet and outlet ports is required for each chemical. Fixing
C*(x,o,t) and C*(x,Ly,t), 0 < x < L*, for all time t > 0 will not satisfy
the accounting condition (Parker and van Genuchten, 1984). Furthermore,
attempting to fix the concentration at y — Ly is not possible in most
"At	^	^
C (x,y,o ) = g (x,y) ,
(5.2.1)
where (x,y) is in D so that
L^ J J C (x,y,o+)dxdy <
(5.2.2)
35

-------
situations. Iri order to satisfy the mass accounting requirement, a flux
type concept along the inlet and exit port boundaries is introduced.
This concept allows conservation of mass better than would be possible
by fixing the chemical concentrations at the inlet and exit port. Fig-
ures h and 5 show aspects of this concept of the mass conserving
boundary conditions at the inlet and exit ports.
A detailed discussion of the development of the mass accounting
boundary conditions at the inlet and outlet ports, induced by the finite
size of the ports, was given by Lindstrom and Boersma (1990) for a
single dissolved chemical. Thin boundary layers are introduced at each
of the porous medium inlet and outlet port interfaces. Assuming that
the four dissolved chemicals are noninteracting in either the inlet or
exit ports, the inlet and outlet port concepts presented by Lindstrom
and Boersma (1990) for one compound can be used for each of the four
compounds. The boundary conditions at the inlet and outlet ports can
then be written as
C (x,0,t)
~
(£ D )
yy
* *
C (x.iy t) ,D
" ~ + w~ * « b-->
v—0 &yl
in
y y=°J
*
c.
m
*
(£ D )
yy
v-0
Ayi
+ « Vy-0
cin

in
(5.2.3)
and
C (x, L , t)
*
D
c u + -p-
y Ny
*
C (x,L
cout _*
Ay,t) + ~	C
J	Ay	out
J f
y-K
out
t D
e U + —
y Ay
* >
)
jj.
cout
Ny . t
y=L„
Ay
out
(5.2.4)
36

-------
open lace ts	outlet end
mixing cnamber
.njediorv'eirtraaion well
l-ilel end
mixing cnamber
CC
x = Lx
Figure. 4. Schematic diagram of the physical aquifer used for testing
the mathematical model, showing mixing chambers at inlet and
outlet ends.
wwll atirr»d outtrt
mixing t»oh	
CO
CO
Covti DquI> ^owt

—H
Figure 5. Details of the use of the "Gaussian pill box" concept in the
mixing chambers, showing flow orientation, and positions
where mathematical terms were defined.
37

-------
It is easily seen from these two equations that in the limit as D*cir and
D*cout large, then C*(x,0,t) -+ C*lr) and C*r(x,Lyit) -» C*0.Jt,( i.e. with the
assumption of well stirred conditions in both the inlet and outlet
chambers, the diffusion gradients within the ports themselves vanish.
The differential equations defining C\n and C*out were derived in
detail by Lindstrom and Boersma (1990). They find that at the inlet
port
*
dC.

dt
.* 1	r x		1
¦o rr J <€yu )dx + rr
in x 0	y y=0 in x 0
L £ D
Ayi
C(x,Ay^,t)dx
y=0
L. L
in x
J
€ D
-2X.
Ay^
+ e U
dx
y=0
*
C.
m
(5.2.5)
and at the outlet port
dC
out
dt
*
1 L r E D	-V
— r x I—XX- +	e u I
LoutLx 0 } AyNy	y^y=L
C U'Ly * ^Ny'^
dx
1 I rLx fe Dvy + £ y ]
^out LX [0	^ ^yNy	y^y=L
dx
out '
(5.2.6)
where C*in and C*out at time t = 0 must be specified.
Along the planes x = 0 and x = Ljj the zero chemical flux (Neumann)
conditions
dC
dx
*
= 0 =
dx
(5.2.7)
x-0
x-L
38

-------
are required. These two conditions obtain directly from applying a
Gaussian pill box concept to the impervious, bounding, walls and
recalling that Ux — 0, along both walls as well.
The total chemical mass in the system at time t (days) following
initiation of the transport and fate processes is defined by
L L
*	y x * *
M (t) - L J J «(1 + R ) C (x,y,t)dxdy ,	(5.2.8)
Ati	w 0 0
where L„ is the vertical dimension (m) of the aquifer.
The total cumulative chemical mass entering (+) or leaving (-) the
inlet port is defined by equation (5.2.9),
L
Minlet(t) " Lv / / 
y=0
The total cumulative chemical mass leaving (+) or entering (-). the
exit port is defined by equation (5.2.10),
MoutUt(t) " K J J <"j * 'I* | )dxdt '	(5.2.10)
y-Ly
The total cumulative chemical mass "lost" by all first-order processes
operating in the system is defined by equation (5.2.11),
*	**	* *	Q r*
M-. C^(c) ~ L J J J {£ A C(x,y,r) + ~2iLE Jdxdydr. (5.2.11)
0	0 0 0	"w
The total cumulative chemical mass gained from distributed zero
order sources, e.g. leaking sources, is defined by
t L Lx
Curce(t) = Lw S I y S {^so(x'y'r) + Qinj C*'dxdyd' • (5-2.12)
0 0 0	J
39

-------
The five scalar mass quantities, which are computed in the model LT3VSI
give a complete mass accounting for the aquifer.
5.3. Approximate Solution for the Transport and Fate Equations
The approximation is initiated by recalling the subregion in space
Dij U B, which contains the i,jth nodal point of both the flow field and
the chemical transport field (Figure 3). The time coordinate [0,T] is
partitioned into a set of adjoining time subintervals so that time tn =
ta-i + Atn, n = 1,2,3,...	Define
fi*. - fi.. U'O.T)	(5.3.1)
ij	ij '
to be a space-time subcylinder so that
D* = u {D* }	(5.3.2)
1 » J 1J
is the global space-time cylinder. Then,
D*-
40

-------
y
= T
x = 0
x = L
I *_ *n+l
^n,n+1 = MCDij U[t n,t n^_] ])
x
Figure 6. Relations of the subregions of the overall space-time
cylinder to each other and the superimposed set of lattice
points.
Lindstrom and Boersma (1990) give the necessary details of the
approximation procedure where the nodal points are chosen in the mesh
rather than in the center of each block. N'ow assume that the four
compounds do not alter the flow field and that they may be treated
independently of one another. Thus, using the previously developed
equations, the generic equation approximating the transport and fate of
each compound can be written as
41

-------
*
C- 4
i,j,n+l
C (P.W J + 	 +	 C, • +
lj(n) 1 + R ,J'n
ij
At
*
1 + R. .
A (xi>yj,tn)
- aitl*. C*	+ ait2*. C*	+ ait3*. C* , . ,
lj l -1, j -1, n	ij i-l,j,n	lj l-1,j+1,n
+ adtl,. C, . . + adt2,. C, . + adt3.. C. . ..
ij i, J -1 ,n	ij i.j.n	ij i,j+l,n
* tV	•* ¦*
+ autl. . C, , . .. + aut2.. C, .. . + aut3.. C, , . ,
ij i+l,j-l,n	ij i+l,J,n	ij i+l,j+l,n
£ij(1
:i + r..)1
ij
*	r *
Q	+ Q. . C
so. .	mi . . I sour
ij	Jijv	ij
[c* -	,
4 I sour	p JJ
(5.3.4)
where the nine coefficients a^tlx j ait2j j, etc. are defined as follows:
ai tl. . =
ij
(5.3.5)
(e D . /0 . + e D . . -)/4
	XV1-1/2.1	XVI.1-1/2	
At
{(AxiL + Ax^/2} ( (Ay j __ 1 + Ay )/2)
/ [<..(1 + R*.)]"1
ij	ij
ait2. . =
J-J
(5.3.6)
At-
"A"	^ ^
rr&y. i + Ay.^ eD . - ,0 .	(eD . . - /0-eD . . . ,n)-i
I ] -1	_j.| xxi-1/2. i	xvi. 1-1/2 xvi. 1+1/2
2 Axi_i	+ 4
j(AXi_i + Ax1)|j(Ay._1 + Ay.)
/[4ij(1+Rij)1'
ait3,
ij
(5.3.7)
•At
(c D . . + t D . , /0 .)
xvi .1+1/2 	xvi-1/2 . i
f(Ax. + Ax.
i—^

1 + AV
/ i-ijd ~ Ry)] ¦
42

-------
*
adtl.. =
ij
(5.3.8)
At
"k	~ *¦
r/-Ax. , + Ax.-v cD . , ,	(eD . , ,-eD . , .
i-l	i1 vvi. 1 -1/2	xvi-l/2.i xvi+l/2.i
2 J Ay. ,	4

/[.^(URy)]
*
adt2..
l-j
(5.3.9)
'[Ayi-1 * Ayjj [£ Dxxi+l/2,j
-At
e D
xxi-1/2
Ax.
	I
AXi-l

{/[e. . (1+R. .)] ,
' lj	lj
autl.. -
ij
(5.3.11)
-At'
(e Pxvi+l/2.i + £ Dxvi.i-l/2)
{Axi-1 + ^V^i-l +
/ [ e i j (1 + Rij)] '
43

-------
¦*
aut2..
(5.3.12)
~	XX
r ,-Ay. , + Ay.A eD . , .	(«D . . , ,„-(D . . , ,„K
f i -1 Ji] xxi+1/2.i	xvi.i+1/2 xvi.i-1/2'
^ 0 J Av	/.
At'
Ax.
l
|(Ax. ! 4- A»t>^(ayi-1 + Ay.)j
/Kyd^j)].
*
aut3..
ij
At'
U D?cyi+l/2,j + f Dxyi,j+l/2)

(Ay^
i-i + ayi'}
/Kijd-Ry)
(5.3.13)
P*ij(I1j is one of four points within , which is viewed as the position,
at time tn of that particle of compound * which lies on the trajectory
or fluid flow path, passing through (x^yj) at time tn+1. A detailed
explanation of how the node P"iJ(n) is found is given in Appendix A.
With a good estimate of P*ij(r) now available C*(P*1j n) is estimated,
using the interpolation formulas given in Appendix A.
The mid-nodal point estimates of the various transverse, cross-axis,
and longitudinal dispersion coefficients are defined to be
* *
.	e D .... + € D . .
f D	i —xxi±Lj	ma	(5 3 14)
xxi+1/2, j	2	'
~
e D . 0 .
xxl-1/2 , j
*
€ Dxyi+1/2,j
*
e D .	.
xyi-1/2,j
+ *
e D + £ D . .
xxi. 1	xxi-1.1
* *
e D . . . + e D . .
XV 1 + 1 . i	xvi. ]
* *
£ D . . + e D - n .
	xyi, j	xy i-1, j
(5.3.15)
(5.3.16)
(5.3.17)
44

-------
* *
.	.	c D . . + c D . .
< D* . . , „ - < D* , . , „ - 	SWiL	,	(5.3.18)
yxi,j+l/2	xyi,j+l/2	2
*Je
.	.	£ D . + £ D . . .
( D	= e D	= 	X^1' 3	xyL '1	(5 3 19)
yxi.j-1/2 xyi.j-1/2	2	'	O.J.i*;
£ D*	+ £ D* . ,
e D	^ —xyiu±l	izLd	(5 3 20)
yyi,j+l/2	2	'
•k	*
«	t D . + e D .
t D . . 1 - 	 —7	(5.3.21)
yyi,j-l/2	2
5.4. Boundary Conditions for Approximate Concentration Field
To complete the description of the approximate chemical concentra-
tion fields, C'ij n, the initial and boundary data must be specified.
5.4.1.	Initial data
At time t - 0 days
C* . - C* (x., y., 0+) ,	(5.4.1)
i,j,o	i Jj'	'
completely specifies the initial data. Recalling the geometric inter-
pretation of C'ij j, this means that over each region D-j at time t = 0,
the discrete approximation begins with the concentration of each
chemical being constant in the space-time cylinder, with the value of
the true concentration distribution * evaluated at nodal point (xiPyj).
5.4.2.	Inlet and outlet end data
Recalling the definition of the inlet boundary concentration,
C*in(t), given by equation (5.2.5), define the constants
45

-------
Uinl
r s (t vdx
x 0

(5.4.2)
U
in2
1 x
tj <
x 0
t ^ + e U )dx ,
A-yi yy-0
(5.4.3)
and function
U!„3
, x D	^
r J <£ a^~) C (x,Ay t)dx ,
x 0	yl y-0
(5.4.4)
so that the differential equation defining C"in(t) becomes
U. , U* ,(t) u*
¦*
in _* "inl "in3
		 = c -	 + 		
dt	o L.	L.
in	in
in2 C*
L. in
in
(5.4.5)
This linear,, first order, Bernoulli class ordinary differential equation
can be shown to have the unique time convolution solution between the
two time levels tn and tn+1,
expfU. „t .,/L. ] C. (t ,)- exp[U. t / L. ] C. (t )
Kl in2 n+1 in in v n+1	in2 n' inJ in v n
. U.	U* (t)
r	._* inl in3
- J	^	 + 		) exP
"n+1
o L.
in
r *
U. _ t
in2
L.
in
L.
in
dt
(5.4.6)
Assume that over [tn,tn+1], U*ln3(t) is continuous and one-signed, and
exp [U*ir2t/Liri] is continuous and dif f erentiable. This yields, upon
multiplication of both sides by the expression exp[- L'*in2t/Lin] . the
formula
C. (t ) - exp[-U. At/L. ' C. (t )
m n+1	rL m2	in" in n
¦k ic	•-k
[C U. + U. (f )1
1 o inl	i n3 n	-	*
in2
1 " exp[',;iI,2 "'''"in1] '
(5.4.7)
46

-------
where f*n is some value of time t in [t^tj,^]. Since fn is not known,
its value is approximated by setting f*n = t^,, thereby making equation
(5.4.7) an explicit formula, which is in line with the explicit
computing formula for the approximate free phase chemical concentration
distribution, given by equation (5.3.4). The formula actually encoded
is then based upon the recursive scheme
*in(n+l) = eXp["Uin2 At/Lin] Xin(n)
[Co Uinl + Uin3 (tn)]	*	.
+ 	rZ	 (1 ¦ exp:"Uin2 in ' (5'4-8)
in2
where Xin(n)	approximation to C*ln at time tj,. The initial
condition required to get this recursive scheme going is x*in(o) =
C*in(0) where the zero argument of Cin stands for the zero point in
time.
Since Uy(x,0) is not known, except at each nodal point, each of the
three integrals stated by equations (5.4.2) through (5.4.4) is approxi-
mated by the trapezoidal rule as follows:
L	N -1
x	x
; (e U ) dx - £ I1" ,	(5.4.9)
0	yy=0 k=l
where
Ax,
Jil ^ ~T [f uy(V0) + £ VVr0)1-	(5.4.10)
and
Lx D*	"x'1
f (£ Zf" + f Uv)dX = S J2k '	(5.4.11)
0	Dyl	yy-0 k-1
47

-------
where
in
2k
Ax. r £ D
[t u
(	*X"
Ayl
£ D
+ e u )	+ (
y (x ,0) Ayl

+ £ U )
(xk+l'0)
(5.4.12)
and
X £ D
S (—»
0
N -1
X
-) C(x,Ay t)dx = Z I
yl y=0	k=l
(5.4.13)
where
Ax.
. m 	k
"3k = 2
£ D
JO_

£ D
C (VAylft) + Ay
(xk>°)
JOL
C (xk+1)Ay1,t)
(Vi-0)
(5.4.147)
U*ini and U*ll)2 are calculated only once in time, but U*iIl3 must be recal-
culated at every time step.
For the outlet tank and/or boundary condition, use equation (5.2.6)
to define the constant
L
x ft D
u* - J- r f—
°utl Lx 0 ^ %y
y-K
+ (e U )1 dx
y )y=K
(5.4.15)
and the function
*	1	r
Uout  - - /
L
x ff D
yy
x o
Ay,
Ny
+ (e U )	C (x,L -Ay t)dx,
y=L	y }y=L	y y
y	y
(5.4.16)
so that the differential	equation defining CQut ( t) becomes
•k -k	*
dC U	U n(t)
out _ outl	out2
dt L out	L
out
out
(5.4.17)
48

-------
This first order Bernoulli ordinary differential equation (ODE) can be
treated in the same way as the ODE for the inlet tank. Over the two
time level interval [tn,tT1+1]
exp[U . t ./L ] C (t .. ) - exp[U .. t /L ]C (t )
outl n+1 out out n+1	outl n out out n
,Vl Uout2(t)
= r _2M£^	 ry t/L j d	(5.4.18)
J	L	outl ' out
t	out
n
Applying the mean value theorem again obtains, upon subsequently
multiplying through both sides by exp [ -U*ou.^L^t) tr+1] ,
C* (t .) - exp[-U* . At/L } C* (t )
out n+1	outl	out out n
( ^ )
_ou^2—n (1 . expr.u At/L 1) ,	(5.4.19)
U ,	Kl outl ' out '
outl
where is a time point in the interval [tn.t,,.!]. As before, rjn* = tn,
is chosen to generate an explicit formula, approximating
X*out(n+i) so that
*out(n+l) exPt Uoutl At^Lout^ *out(n)
Uout2(tn)
+ —	 (1 - exp[-Uoutl At/L^]) .	(5.4.20)
Uoutl
Initially, set x'outco) ~ C*cut(0) , and apply equation (5.4.20) recursive-
ly, in conjunction with the evaluation of U*0._,
As was done in the case of the inlet boundary, approximate
L	*	N -1
x rt D	x
dx - S I™* ,	(5.4.21)
y=L k=l
¦ y
X f(
/ t
0 ^ %y	y
49

-------
where
rout
lk
2
re D
> Ny	y
)l + ('

-------
AX»7	"j.	t	AX„,	fiX,
Nx
i c*	+ (2 * N* + ^SsJ,] C*
c*	axN«.l Nx-2.J,n [ Ax^, A«„y J Nx-l.J.n ^
Nx.J.n _	^ ^ ^,.1	' '
AxKx
for j = 1,2,... N,
y-l ¦
5.5. Approximate Microbial Populations
The microbial populations at each point in space and time satisfy
the first order process laws stated by equations (5.1.11) and (5.1.12)
3N
^T~ (x,y , t) = ui^(x,y, t)N^	(kg kg"1 day"1)	(5.5.1)
and
aN
— (x,y,t) = w2(x,y,t)NT2	(kg kg"1 day*1)	(5.5.2)
where the two first order rate coefficient functions are defined by
w, - Y1 r1 - 51 + Y1 . r1 .	(5.5.3)
1 so sonu	sm sninu
and
2 2	2
u„ = Y r - 6 .	(5.5.4)
L so sonu
Assuming that at time t = 0 days, N1(x,y,0) and N2(x,y,0) are both known
initial distributions, the formulas
t
Nl(*.y.Vl) =	exPl f n+ (X, y, r) dr ]	(5.5.5)
t
n
and
51

-------
N2(x,y,tn+i) - N'j (x ,y, t ) exp[ J n+1 (^(x.y, r)dr ]	(5.5.6)
t
n
are easily obtained. Application of the mean value theorem to the
integrals in equations (5.5.5) and (5.5.6) yields
N1(x,y,tn+1) = N1(x,y,tn) exp[At ^(x,y,rln)]	(5.5.7)
and
N2(x.y.tn+i) = ^2 (x' y * tn) exp[At u)2^X,y,T2n^ *	(5.5.8)
where both Tln and r2n are values of time t in	Up to this
point no approximations have been made in the development of the
representations of the two microbial populations over time. By choosing
both Tla and r2n equal to tn, the approximating populations are defined
recursively to be
N1(x,y,tn+1) = N1(x,y,tn) exp(it u^x.y,^)]	(5.5.9)
and
N2(x,y,tn+i) = N2(x,y,tn) exp [ At u>2 (x, y, t^) ] .	(5.5.10)
It is clear from the definitions of wj and u>2 and the mathematical
structure of formulas (5.5.9) and (5.5.10) that whenever and u2 ar©
greater than zero the two populations are increasing and unstable.
Whenever and u>2 are equal to zero, the populations are constant, i.e.
stable at their respective values whenever ci^ and to2 were first zero.
If and oj2 are less than zero the populations are decreasing (stable)
and will eventually cease to exist. Equations (5.5.9) and (5.5.10) are
used in the model LT3VSI.
52

-------
5.6. Stability Considerations
Since the computing scheme defined by equation (5.3.4) is explicit,
a stability criterion arises. Richtmyer and Morton (1967), Varga
(1962), and Reddell and Sunada (1970) have established that by choosing
0.5	
*
At0 *
(5.6.1)
rD
max
i.j
D
xxij + yyij
Ax.
Ax2
J
9
the computing scheme will be stable when the local Courant numbers are
less than or equal to 1.0 and there are no losses or chemical injections
present. However, Cheng et al. (1984) claim that a sufficient condition
for stability is
	1	
*
A^ <,
max
i.j
*
|U . .
' xii
*
|U . .
(5.6.2)
Ax.
Ay
^ + 2(^U + ^U]
l	v Ax.	Av,
i	j
which is a tighter restriction than that given by equation (5.6.1).
Because the flow field is at steady state, it is a simple matter to
calculate this value of At before the start of the chemical transport
and fate computation. The At value, obtained initially so that the
local Courant numbers are less than or equal to 1, is defined by
At < min
i.j
Ax.
-1

2 IU
2 IU
(5.6.3)
This condition guarantees that the stream line and the particle of
interest "riding" on that stream line, going through (xi,y1) at time
tn+1 must also pass through P*ij(n) at time t^,. The addition of distrib-
uted first order losses and the presence of injection wells further
53

-------
reduce the magnitude of At. These important restrictive conditions will
be satisfied by choosing At so that
* <
At.
0	+	*	*	*
vini..	21U ..I 2IU ..I rD ..	D
. 	ikJ	 ,	. VII , 7\ XX11 yyij
Ax- Ay- L 2 2
i	1	vAx.	Ay.
i	'j
max<, *	,, * .
.. ]1+R.. e(1+R..)p
1J I	w
(5.6.4)
1
The presence of the complex rates for utilization of chemical compounds
by microbes rules also affects the choice of the maximum allowable time
step size. Careful inspection of equations (5.1.21) through (5.1.24),
which define the loss rate functions, and equations (5.1.8) through
(5.1.10), which define the utilization rates, clearly shows that the
chemical concentration C* appears linearly in the numerator of each loss
function, A*. Since each C* is nonnegative it is easy to define a set
of positive upper bounds upon A*, one for each compound *, which
depends linearly upon C*. For example, the bound on A™
0 < AnU < BnU(K1, N2)Cnu
(5.6.5)
where
.nu,,.	. b
( 1' 2 = T
,1 to- -L-
0 Y1 K1
so onu
,1 "ni
ni ,^rl
sni ninu'
+ jfj
2
° Y2 K2
so onu
(5.6.6)
The other three bounds have similar forms. For stability purposes then
define that
54

-------
At,
raaXil + R*
B (N N )
+ 	:	 +
Qinj..
	L±J-
2|U*..I 2|U*..
+ xij + yi]
ij
ij
1 + R
ij
rD . . D .
*¦ A v	Air *
Ax.
l
Ayj
f 'w(1 + V
-1
Ax.
1

(5.6.7)
Since B* is a function of the mass of the microbial population per unit
mass of aquifer material which changes with time, we must be very
careful to recompute it/, periodically. A rule for periodically
checking on the magnitudes of the B* values is included in the
Fortran 77 code, LT3VSI.
We then choose
At < m^n (At^)
(5.6.8)
which will be used as the sufficient stability criterion.
6. CONCLUSIONS FOR THE MODEL
Methods were developed for solving equations that describe trans-
port and fate of chemicals in laboratory scale models of aquifers. The
mathematical model is for aquifers consisting of a single layer of
material, which can be either heterogeneous or homogeneous and
anisotropic or isotropic with respect to the water flow field and
heterogeneous or homogeneous but isotropic with respect to the chemical
transport field properties.
The two-dimensional transport and fate model can be used for study
of the important aspects of bioremediation of aquifers contaminated with
55

-------
nitrogen. A broad range of aquifer remediation scenarios may be con-
sidered. These scenarios could include studies of placement of
injection/extraction wells to induce plume spreading or plume shaping
and the effects of regions of varying hydraulic conductivity on the
shape of the plumes. The comprehensive treatment of the inlet and exit
port induced boundary conditions, included with the analysis represents
a significant step forward in modeling the transport and fate of
chemicals in laboratory scale physical aquifers.
56

-------
7. LITERATURE CITED
Bachelor, G. A., D. E. Cawlfield, F. T Lindstrom arid L. Boersma. 1990.
Denitrification in nonhomogeneous laboratory scale aquifers: 2.
User's manual for the preliminary mathematical model LT2VSI,
Special Report No. 865, Agric. Exp. Sta., Oregon State University,
Corvallis, OR. 158 pp.
Baptista, A. M. de M. 1986. Accuracy analysis of the backward method
of characteristics. In: VI. Intnl. Conf. on Finite Elements in
Water Resources, Lisbon, Portugal, June 1986.
Baptista, A. M. de M. 1987. Solution of advection dominated transport
by Eulerian-Lagrangian methods using the backward method of
characteristics. Ph.D. Thesis. Dept. of Civil Eng., MIT,
Cambridge, MA. 260 pp.
Baptista, A. M. de M., E. E. Adams and K. D. Stolzenbach. 1984.
Eularian-Lagrangian analysis of pollutant transport in shallow
water. Ralph M. Parsons Lab Rept. #296, Dept. of Civil Eng., MIT,
Cambridge, MA 02139.
Bear, J. and A. Verruijt. 1987. Modeling Groundwater Flow and Pollu-
tion. Reidel Pub. Co., Dordrecht. 414 pp.
Bloomfield, P. 1976. Fourier Analysis of Time Series: An
Introduction. Wiley, New York. 258 pp.
Bodvarsson, G. 1984. Linearization techniques and surface operators in
the theory of unconfined aquifers. Water Resour. Res. 20:1271-
1276.
57

-------
Boersma, L. and F. T. Lindstrom. 1990. Denitrification in nonhomogene-
ous laboratory scale aquifers: 3. Simulated application of the
prelimninary model (in preparation).
Cheng, R. T., V. Casulli and S. N. Milford. 1984. Eulerian-Lagrangian
solution of the convection-dispersion equation in natural coordi-
nates. Water Resour. Res. 20:944-952.
Herbert, D. 1958. Some Principles of Continuous Culture. In: Recent
Progress in Microbiology. Edited by G. Tunevall, Blackwell
Scientific, Boston, MA.
Konikow, L. F. and J. D. Bredehoeft. 1978. Computer model of two-
dimensional solute transport and dispersion in ground water. In:
Techniques of Water-Resources Investigations of the United States
Geological Survey, Book 7, Chapter C2.
Kreyszig, E. 1979. Advanced Engineering Mathematics, 4th Ed. Wiley,
New York. 939 pp.
Lindstrom, F. T. and L. Boersma. 1990. Denitrification in
nonhomogeneous laboratory scale aquifers: 1. Preliminary model for
transport and fate of a single compound. USEPA RSKERL (in press) .
Lindstrom, F. T., L. Boersma, M. A. Barlaz and F. Beck. 1989.
Transport and fate of water and chemicals in laboratory scale,
single layer, aquifers. Vol. 1: Mathematical model. Spec. Rept.
No. 845, Agric. Exp. Sta., Oregon State University, Corvallis, OR.
52 pp.
Molz, F. J., M. A. Widdowson and L. D. Benefield. 1986. Simulation of
microbial growth dynamics coupled to nutrient and oxygen transport
in porous media. Water Resour. Res. 22:1207-1216.
58

-------
Parker, J. C. and M. Th. Van Genuchten. 1984. Flux-averaged and
volume-averaged concentrations in continuum approaches to solute
transport. Water Resour. Res. 20:866-872.
Prenter, P. M. 1975. Splines and Variational Methods. Wiley, New
York. 323 pp.
Reddell, D. L. and D. K. Sunada. 1970. Numerical simulation of disper
sion in groundwater aquifers. Colorado State Univ. Hydrology Pape
No. 41. 79 pp.
Richtmyer, R. D. and K. W. Morton. 1967. Difference Methods for
Initial-Value Problems. Interscience Publishers, Inc., New York.
405 pp.
Stone, H. L. 1968. Iterative solution of implicit approximations of
multidimensional partial differential equations. Soc. Ind. Appl.
Math. J. on Num. Anal. 5:530-558.
Varga, R. S. 1962. Matrix Iterative Analysis. Prentice-Hall,
Englewood Cliffs, NJ. 322 pp.
Vemuri, V. and W. J. Karplus, 1981. Digital Computer Treatment of
Partial Differential Equations. Prentice-Hall, Englewood Cliffs,
NJ. 449 pp.
Widdowson, M. A., F. J. Molz and L. D. Benefield. 1988. A numerical
transport model for oxygen- and nitrate-based respiration linked t
substrate and nutrient availability in porous media. Water Resour
Res. 24:1553-1565.
Yates, S. R. 1988a. An analytical solution to saturated flow in a
finite stratified aquifer. Ground Water 26:199-206.
Yates, S. R. 1988b. Seepage in a saturated-stratified aquifer with
recharge. Soil Sci. Soc. Am. J. 52:356-363.
59

-------
APPENDIX A
The methods used to find the particle tracking point P*i ^ and to
subsequently estimate the chemical concentration C(P*ij n) are described
below. This presentation requires discussion of the double Fourier
integral representation of a specialization of the generic transport and
fate equation. The two-dimensional procedures used in this analysis are
extensions of the one-dimensional procedures developed and discussed by
Baptista et al. (1984), Baptista (1986), and Baptists (1987).
Figure A.l shows the subregion enclosing the node (x^yj) and the
eight surrounding nodes with their positions defined. Ux* and Uy* are
signed quantities and are once continuously differentiable functions of
the space coordinates everywhere, except perhaps on top of the wells.
The procedure for finding P*\jn is described first.
A.l Locating, point P*;j„
The selection rules for finding point P*ijn are defined as follows:
i) if Ux* > 0, Uy* > 0,	lies in subregion I
ii) if Ux* < 0, U/ > 0,	P*ij(n) lies i-n subregion II
iii) if Ux* > 0, Uy* < 0,	P*i.(n) lies in subregion III
iv) if U*x < 0, Uy* < 0,	lies in subregion IV
The details of the estimate for the position of P*i;(ri; are described for
subregion 1, but similar arguments apply to subregions II, III, and IV.
Since Ux* and Uy* have been computed at every nodal point in £5, Ux*
and Uy" are known at (x^yj) and the eight surrounding nodal points (Fig-
ure A.l). Figure A.2 shows an enlargement of region I with flow lines
60

-------









A _
j Flo*
< u, 3
u,!
orientation





y
~



AJl-l
1
r\
11
71



¥
ID
Y.
r>
N




^i-i
AZ|
I













*I.J	*1.1	*1	*1.1	*U)
Figure A.l. Subregion Did divided into four separate subregions, one of
which contains point P*ij(n). Flow orientation is also
shown.
.\
\

jj.j
t
Ayi-i
I
pT
l\ Ui^

y
•K
h
|<	Aljj	*•

*n	*i
Figure A. 2. Details of procedure used to find position of point P ij(n)
in region I.
61

-------
sketched in. The effective seepage velocity components U*x and U*y are
approximated by the following bilinear interpolation formulas:
U* - aJ + bJ(x - Xi) +' cj(y - y ) + dJ(x - x^y - y ) ,	(A.1.1)
(A.1.2)
= A2 + B2(X " V + C2(y " yj} + °2(X " Xi)(y " yj} '
where the coefficients are defined by:
T *	I *
A. - U . . , A, - U . . ,	(A.1.3)
1	xi,j	2 yi,j
* * * *
U..-U...	_ U..-U...
xi.i xi-l.i	_I vi.1 vi-l.i	. ..
=	Ax ,	'	2 °	Ax. ,	 •	(A-1-4)
i-1	1.-1
u . . - u . . .	_ u . . - u . .
_		^ul	(A-1.5)
Ayj-1	2 Ayj-1
*	"k	*	*
U..+U....-U...-U...
_ xi. i	xi -1. i -1	xl -1. i	xi. i -1
Ax. . Ay.
1-1 J-1
&	-k	"k	-k
u . . + u . . . . - u . . . - u . . .
d1 - 	yi-i.j-i	yi-i,3	yi.j-i	(A1 6)
2	Ax. , Ay. ,	'
l-l J -1
so that both Ux" and Uy* interpolate Ux* and Uy* exactly at the four
corner points.
Since
f - U* and & - t* ,	(A. 1.7)
dt x	dt y
integrating both sides over [tn,cn+1j yields
tn+i *
X(tn+1} = X(Cn} + S Ux dt '	(A.1.8)
62

-------
tn+i +
y(tn+i) - y(tn) + J Uy dt ,	(A.1.9)
where now [x(tn+1) , y(tn+1)] = (x^yp and [x(tn) ,y(tj ] = P\j(n) - [x*iCl0,
The relationships (A.1.6) and (A.1.7) are approximated in a
manner analogous to the low order accuracy, time quadrature, approxima-
tion used in the main transport and fate equation (equation 5.3.4):
x. = x.	+ At U	(A.1.10)
11	x
and
y. - y.	+ At U .	(A.1.11)
J l J l	y
Substituting for Ux* and Uy* leads to the coupled system of nonlinear
algebraic equations in (x^y^*):
fl(xi'yj) ~ Xi " Xi + At[Al + 3l(xi " xi) + Cl(yj " yj)
+ dJ(x* - x.)(y* - yj)] - 0	(A.1.12)
T * &	*	TT^"	T
f2(xi.yj) = yj ' yj + At[A2 + 2(xi ¦ xi) + c2(yj - yj>
+ ^(XjL " xi^yj " yj^ = 0 '	(A.1.13)
The method of Newton-Raphson-Kantorovic (Vemuri and Karplus, 1981) is
used to find (x^, yj*), with (x^y^ as the initial guess and a relative
error tolerance of 10~*. The theorem, which states that the nonlinear
system of equations (A.1.12) and (A.1.13) has a unique solution inside
subregion I, for the condition that A1!, AI2, B1!, BI2, C'j, Ciz, D1!, and
Dxz are all finite and that At is chosen so that
63

-------
At < min-
Ax. ,
. 2 IU I
Ay
a
2 IU
(A.1.14)
y
can readily be proven. This proof is not shown here.
A.2 Fourier Component Analysis Under Special Hypothesis
Suppose that instead of working with a finite size domain D, we let
the side walls and the inlet and outlet ports recede to infinity. This
is an acceptable assumption for conditions where regions influenced by
the injection/extraction wells and/or regions with initial distributions
of compound are small relative to the longitudinal and transverse
dimensions of the aquifer. For a homogeneous and isotropic aquifer with
a velocity field that is constant in time and space, the generic
transport equation (5.1.25) reduces to the form
dC
dt
2	2	2 1
*	a c * s c * a c	*
D 7T~ + 2D 6-7- + D *-?-}> / (1 + R )
XX dx Xy y yy dy J
~	*
+ ^so + ^ini Csour
«(1 + R*>
(A.2.1)
under the assumption that loss processes are not operating in the
aquifer. Q*ao. Qin< > an<^ Qout may be considered to have a zero value
everywhere in the space time cylinder D+, where
D+ = D u [0,T).
without loss of concept.
The next step is to define the double Fourier transform of some
function f(x,y) which satisfies Dirichlet conditions (Kreyszig, 1979) to
be
64

-------
+as +o°
F(f ,»j) - / / f(x,y) e"i(fX+,?y:) dxdy	(A.2.2)
where i - 7-1, and to specifically assume that the initial distribution
of chemical * is double Fourier transformable, i.e.
C*(r,»?,0) - J J C*(x,y,0+) e*l(rX+,?y) dxdy .	(A.2.3)
-CO -w
The double Fourier transformation of equation (A.2.1) obtains the form
9 4-	0 *3e
(f D +2fr?D +T) D )t
XX	XV	vv	^	^
*	-ifV t -ir?V t
*	*	+	1 + R	1 x	' v
C (f,u,t) - C (f.^,0 ) e	K	e X e y .(A.2.4)
Since the inversion double Fourier integral can be written as
f(x.y) - J J F(f.rj) el(fx+*y) dfdr?	(A.2.5)
4tt -«o -oo
the solution to the sourceless form of equation (A.2.1) can be written
C*(x,y,t) = J J C*(fl,,t) el(rX+f?y) drdr, .	(A.2.6)
4 n -00 - ®
Substitution of definition (A.2.4) into (A.2.6) yields the double
Fourier integral (series) representation of the solution to the infinite
2-space Cauchy problem (A.2.1) as
0 -A-	0 &
(f D +2fr?D +r? D )t
XX	XV	vv
C*(x,yft) - J J C*(r ,»|f ,0t e	1+R
4tT - a> - oo
if (x-V*t) + i»j(y-V*t)
e	y dfdr? .	(A.2.7)
65

-------
Equation (A.2.7) is the linear combination of a continuum of two
dimensional right moving damped waves. The Fourier variables f and rj
have the dimensions of inverse wavelength and are referred to as the
wave numbers in the x and y directions, respectively. The Fourier
transform F(f,fj) relates directly to the energy spectrum of the wave
train f(x,y) and is therefore referred to as the "power or energy
spectrum" by the electrical engineers (Bloomfield, 1976). If	D*xy,
and	all equal zero, then equation (A.2.7) is the solution to the
embedded two dimensional hyperbolic (wave) equation obtained by setting
these same three dissipation parameters equal to zero in equation
(A.2.1).
In the remainder we are interested in the double Fourier
representation evaluated at an arbitrary grid point (x, ,yn) in D. With
time partitioned into discrete finite length subintervals of length At
and with a uniform partitioning of the x-axis into subintervals of
length Ax and the y-axis into subintervals of length Ay, we have
2 *	* 2 *
(C D +2fnD +77 D )nAt
xx	xv yv	
C*(xrym'cn) = ^2 f f C*(
-------
where the amplitude per time step factor RE(J",r/) is defined by
'(f2D* +2fr?D* +r?2D* )nAt|
XX	XV	vv
RE(r.»?.) = e
*
1 + R
and the total phase angle per time step	is defined by
$E(r.»7) = crx$"Ax + crynAy =	+ 4>Ey(r,'?)	(A.2.10)
with Crx being the x-coordinate Courant number
+
V	At
C =	.	(A.2.11)
rx Ax
and Cry being the y-coordinate Courant number
*
V	At
C - -J£— .	(A.2.12)
ry Ay
Quantities $Ex/J"At and $Ey/rjAt have dimensions of effective seepage
velocity in the x and y direction, respectively. The reason for
rearranging the exact solution into the form given by equation (A.2.8)
is that much the same form is obtained for each numerical approximation
method chosen (both finite difference and finite element), with the
differences coming from the definitions of the amplitude per time step
and the phase angle per time step. Distribution (A.2.8) is the standard
or reference distribution when the amplitude and phase angles per time
step are defined by the equations (A.2.9) and (A.2.10). The definition
of the two Courant numbers are universal and hold for both the exact and
approximate solutions.
Two alternative methods for estimating C*(P1jn) have been encoded in
the numerical model LT3VSI. An in-depth Fourier analysis of the first
67

-------
of these methods is now considered, using the above outlined specialized
hypothesis.
A. 3 Alternatives for Estimating C*(P.j..)
A.3.1 Alternative 1: Bilinear Lagrange Polynomials
Since point P*ljn is rarely, if ever, a nodal point, C*(P11n) must "be
estimated using an interpolation formula. The logical first choice for
an interpolation formula is the bilinear Lagrange polynomial. This is
shown by continuing with the assumption that P*^ lies in subregion 1 in
Figure A.l.
Approximate C*(Pftijn) using the formula
_* *	 _*
(*i " *i-l)
Formula (A.3.1.1) has been obtained by applying the product basis
assumption and the Lagrange canonical basis functions for the set of all
piecewise defined polynomials of degrees one or less (Prenter, 1975).
Formula (A.3.1.1) produces estimates which are easily shown by standard
multivariate Taylor series analysis to be OCAx2^ + Ay2j-i) accurate,
which is in line with the low order accuracy of the approximate discrete
68

-------
transport and fate equation. The other three subregions yield similar
formulas with the particle tracking concept being carried over from one
region to the next. The Fortran 77 encoded program LT3VSI contains the
formulas for all four subregions induced about every interior node.
Unfortunately the bilinear Lagrange interpolation procedure,
defined as alternative 1, carries with it a great deal of numerical
diffusion (Baptista et al., 1984), which can be reduced only be using
very small Ax, Ay combinations. But, close spacing of the nodes results
in a very large number of nodal points which means that large systems of
equations must be solved. Thus the procedure has the disadvantage of
being slow and limited by computer storage and speed. The problem of
numerical diffusion is magnified in regions of high seepage velocities.
Reasons for the large numerical diffusion associated with this
interpolation procedure are discussed in detail by Baptista (1987) for
the one-dimensional case.
Since the flow field is assumed to be constant, the coordinates of
the point P*ijn are defined by the equations
*
x
i
(A.3.1.2)
and
*
V At
y
(A.3.1.3)
Substitution of these quantities into the bilinear form (A.3.1.1)
obtains
C C. . . . + (1 - C )C C. , .
rx ry i-l,j-l,n	rx ry i,j-l,n
(A.3.1.4)
69

-------
At the arbitrary point (x, ,ym) the approximate transport and fate
equation (5.3.4) operating under the special hypothesis becomes
ct	= c c C* 1 + (1 - C )C C* .
i,n,n+l rx ry i-l,m-l,n	rx ry i,n-l,n
+ C (1 - C ) C*	+ (1 - C )(1 - C ) C*
rx	ry iy,m,n	rx	ry i,m,n
*
D At
XV
C,
~
D At
xx
C
. ,,	i-l,m-l,n A 2/t _* i-l,ra,n
2AxAy(1 + R )	Ax (1 + R )
*
D At
_x*	
c
*
D At
w	
—x
c
o. . /t	i-l.m+l.n A 2... *.
2AxAy(1 + R )	Ay (1 + R )
- 2
D At D At
_xx	 i 	
Ax
Ay
/(I + R )
—•A"
Q
I ,m, n
*
D At
yy	
_*
c,
*
D At
—
_*
C,
, 2... _* -2,ir.+l,n _A A ...	i+l.m-l.n
Ay (1 + R ) '	2AxAy(1 + R )
*
D At
xx
Ax (1 + R )
* ^i+1,m,n +
*
D At
*y	
c,
n. » '/i n*-. i+l,m+l.n
2AxAy(l + R )
(A.3.1.5)
A standard multivariate Taylor series analysis applied to equation
(A.3.1.5) shows the total diffusion coefficient in each coordinate to be
of the form
*	*9
* * V*	Vx "
}XX DXX + 2 "	2
tot
(A.3.1.6)
and
D
yy
tot
* V?
D +
yy 2
*2
V At
_J£	
(A.3.1.7)
70

-------
The cross term D*,^ is not altered.
The	Fourier component can be written as
r* . . 	 .	n
Z, m,n(f,,?) Aje,m e
(A.3.1.8)
where A, n is a nonzero amplitude parameter.
Substitution of definition (A.3.1.8) into equation (A.3.1.5) and
subsequent simplification yields
0 ~ C C e-i(fAx+r7Ay) + (1 . c )c e"i^y
rx ry	rx ry
+ C (1 - C ) e"irAx + (1 - C )(1 - C )
rx	ry	rx	ry
*
D At
xv
2AxAy(l + R )
*
-i(^Ax+rjAy) + xx
2	*
Ax (1 + R )
¦if Ax
*
D At
	M	
2AxAy(l + R )
e_ i (f Ax - ryAy) +
D At
—yj—
2	*
Ax (1 + R )
-ir/Ay
¦Jr	*
D At D At
XX
JQt-
A x
*
D At
—¥2	
2	*
Ay (1 + R )
Ay
irjAy
/(I + R )
*
D At
	*2	
2AxAy(l + R )
i(fAx-rjAy)
*
D At
XX
2	*
Ax (1 + R )
iCAx
e * +
¦k
D At
XV	
2AxAy(l + R )
i(fAx+rjAy)
(A.3.1.9)
Simplifying this representation by using the classical Euler formulas
relating complex exponential functions and the circular functions or
trigonometric functions obtains the complex number
- f]_(f .»?) - i f2(r,*7)
(A.3.1.10)
where
71

-------
= (1 - crx)(i - cry) + crx)(i - Cry) cos(fAx)
+ (1 - Crx)(Cry) cos(rjAy) + crx) (cry) cos($"Ax + r?Ay)
2D* At	2D* At
-—^—— [1 - cos(fAx)] - 	22	— sin(fAx) sin(r?Ay)
Ax (1+R )	AxAy(l+R )
*
2D At
- /y ^ [1 - cos(f Ay) ]	(A.3.1.11)
Ay (1+R )
and
f0(r.'?) = c (1 - C ) sin (J" Ax) + (1 - C )C sin(rjAy)
2	rx	ry	rx ry	J
+ C C sin($"Ax + rjAy) .	(A.3.1.12)
rx ry	J
In polar form we may write
-i$ (f ,r?)
6 = RA(5",r?) e *	(A.3.1.13)
where
Ra » + f22)1/2	(A. 3.1.14)
and
*A - T'n"1(i^] " *Ax(C''> +	•	(A.3.1.15)
Clearly, from complex number theory
-in$ (f.rj)
9 = RA(r.»?) e	.	(A.3.1.16)
Back substitution of 8r' into (A. 3.1.8) obtains the frjth Fourier
component
72

-------
-*	n
S m	= Af ,n RA(f-'?) e	•	(A.3.1.17)
x , m, n	X , in A
The approximate solution then assumes the form
^ +® +« ^
c  - — ; f C(r,v,0+) R^r.ij) •
Z17T - oo - oo
	$
Ax	A.v
e	e	dCd rf	(A.3.1.18)
Comparing the structure of RA, defined by equation (A.3.1.14), with the
more simple structure of REi defined by equation (A.2.9) shows that even
with the single 4-point bilinear interpolation scheme, the amplitude for
each time step function is extremely complicated. The phase angle for
each time step function is also a complicated function of and rf, not
simply a linear function of f and rf such as E is. Careful inspection
of equations (A. 3.1.11) and (A. 3.1.13) reveals that whenever $"Ax and r)Ax
are multiples of 2n radians, fj and f2 repeat themselves, regardless of
the Courant numbers or the other parameters which appear in these
equations. This periodicity stands in contrast to the monotonic
function of RE as defined by equation (A.2.9). At the points where $"Ax
and r;Ax are multiples of 2-n radians, RA = 1 and $A = 0 so that the
Fourier components suffer no attenuation and do not propagate in space,
as time evolves. Because both RA and A are periodic functions of and
Tj only the basic rectangular region in IR2 which consist of the set of
all (f ,f|) such that 0 < $* < 27t/Ax and 0 < r/ < 2?r/Ay needs to be
considered. Figure A.3 shows the placement of the basic regions in the
Eucledian plane, IR2. For evaluation of RA and A all other values of
(£",»?). are aliases of this basic region (Bloomfield, 1976).

-------
Since both RA and $A are periodic functions of f and tj , the long
time character of the approximate distribution C*(x,y,t) depends upon
the character or shape of the double Fourier transform of the initial
distribution C* (£ , tj , 0+) , for initial value problems or for the double
Fourier transform of the source function, if any sources exist. Table
A.l shows the values of the functions RE, 3>E, RA, and $A per time step
function RA(f,rj) over the basic region of definition for the following
geometric and transport parameter values: Ax - Ay = 1.0 meter, Vx - 0.0
m/day, Vy = 0.55 m/day, D„ — 0.002 m2/day, - 0.0 m/day, Dyy = 0.002
m2/day, t - 0.001 day. Table A.l also shows the results of evaluating
the phase angle per time step function	over the basic region of
definition. For comparison purposes the exact amplitude per time step
Re and the exact phase angle $E are included in Table A.l. Using the
value of At - 0.001 day requires 20,000 time steps for a simulation 20
days long. Thus, n = 20,000 at 20 days. Table A.2 shows the values of
the amplitudes and phase angles at time t — 20 days for the basic
region. Recall that (Figure A.3) the approximate amplitudes and phase
angles outside the basic region are found from aliasing. The exact
amplitudes and phase angles do not follow this simple aliasing
procedure.
Careful study of Table A.2 shows that by the time 20 days have
elapsed the approximate amplitude factors have decayed to negligible
value except in tight clusters around integer multiples of 2ir in each
coordinate. That is, only in the neighborhood of those (fAx.rjAy) =
(pjr,q7r), p or q or both even integers, is RAn still significant when n =
20,000 time steps. However, these are precisely the Fourier components
which do not propagate at all. Negative signs indicate waves

-------
Basic
Region
2rL'Ax
2TL/Ax
¦2UAy
211
Ay
2Ay
n
Ay
n
2Ay
_CL
n n 3n 2n
2 Ax Ax 2 Ax Ax
->» c
Figure A.3. Placement of the basic region of definition in IR
2
75

-------
Table A.l. Evaluation of Rt, $E, Ra, and $A in the illustrative data corresponding to
a simple constant longitudinal Darcian flow field and a uniform mesh
spacing. The units for RE and RA are (1/time step) and for $E and $>A they
are (rad/time step).
J" Ax

Re


Ra



(time step)"1
(rad/time step) L
(time step)"1
(rad/time step)
0
0
1.000000
0.000000
1.000000
0.000000
0
tt/2
0.999995
8.63938
E-04
0.999448
5.50304 E-04
0
7T
0.999980
1.72788
E-03
0.998896
0.000000
0
37T/2
0.999956
2.59181
E-03
0.999448
-5.50304 E-04
0
27T
0.999921
3.45575
E-03
1.000000
0.000000
t/2
0
0.999995
0.00000

0.999996
0.000000
tt/2
7r/ 2
0.999990
8.63938
E-04
0.998888
5.50307 E-04
t/2
7T
0.999975
1.72788
E-03
0.998888
0.000000
¦x/2
3?r/2
0.999951
2.59181
E-03
0.999442
-5.50307 E-04
tt/2
27T
0.999916
3.45575
E-03
0.999996
0.000000
7T
0
0.999980
0.00000

0.999992
0.000000
77
tt/2
0.999975
8.63938
E-04
0.999438
5.50309 E-04
7T
7T
0.999961
1.72788
E-03
0.998884
0.000000
7T
3tt/2
0.99936
2.59181
E-03
0.999438
-5.50309 E-04
7T
27T
0.999901
3.45575
E-03
0.999992
0.000000
3tt/2
0
0.999956
0.00000

0.999996
0.000000
3tt/2
tt/2
0.999951
8.63938
E-04
0.999442
5.50307 E-04
3tt/2
7T
0.999936
1.72788
E-03
0.998888
0.000000
3tt/2
37t/2
0.999911
2.59181
E-03
0.999442
-5.50307 E-04
3tt/2
2-n
0.999877
3.45575
E-03
0.999996
0.000000
27T
0
0.999921
0.00000

1.000000
0.000000
2 n-
7r/2
0.999916
9.63938
E-04
0.999446
5.50305 E-04
27T
7T
0.999901
1.72788
E-03
0.998892
0.000000
In
3tt/2
0.999877
2.59181
E-03
0.999446
-5.05305 E-04
2n
2ir
0. 999842
3.45575
E-03
1.000000
0.000000
76

-------
Table A. 2. Evaluation of REn, n$E, RAn, and n$A for the data given in
Table A.1.
$" Ax	r?Ay	REn	n$E	RAn	n$A
0
0
1.000000
0.00000
1.000000
0.000000
0
n/2
0.904837
17.27876
0.000015
11.00610
0
7T
0.670317
34.5576
0.000000
0.000000
0
3tt/2
0.414775
51.8362
0.000015
-11.00610
0
2?r
0.205962
69.1150
1.000000
0.000000
tt/2
0
0.904837
0.00000
0.923116
0.000000
ir/2
tt/2
0.818729
17.27876
0.000014
11.00614
tt/2
7T
0.606526
34.5576
0.000000
0.000000
tt/2
3tt/2
0.375302
51.8362
0.000014
-11.00614
w/2
2tt
0.186360
69.1150
0.923116
0.000000
7T
0
0.670317
0.00000
0.852143
0.000000
7T
tt/2
0.606526
17.27876
0.000018
11.00618
7T
7T
0.458399
34.5576
0.000000
0.000000
7T
3tt/2
0.278025
51.8362
0.000013
-11.00618
?T
2tt
0.138055
69.1150
0.852143
0.000000
3tt/2
0
0.414774
0.00000
0.923116
0.000000
3tt/2
it/2
0.375302
17.27876
0.000014
11.00614
3tt/2
7T
0.278025
34.5576
0.000000
0.000000
3tt/2
3tt/2
0.168624
51.8362
0.000014
-11.00614
3tt/2
27T
0.085422
69.1150
0.923116
0.000000
2tt
0
0.205962
0.00000
1.000000
0.000000
27T
tt/2
0.186360
17.27876
0.000015
11.00610
2tt

0.138055
34.5576
0.000000
0.000000
2 7T
3tt/2
0.085422
51.8362
0.000015
-11.00610
2 7T
2tt
0.042415
69.1150
1.000000
0.000000
77

-------
propagating in the leftward direction. It appears that only those
Fourier components whose frequencies fAx and r?Ay fall in [0,7r/2]
propagate reasonably well. The disparity between the phase angles, n4>E
and n$A for in between values of f and ?? is obvious. Clearly as time
evolves, shrinking sized island of area centered about (pjr.qir) being to
appear as the dominant areas which contribute to the actual approximate
distribution C*(xp,yq,tn). The extent of the wiggles depends largely
then upon the magnitude of the tail regions (mid to high frequency
components) of the double Fourier transform of C*(x,y,0+); namely,
C*(f,»?>0+). If distribution C*(,,0+) has considerable high frequency
information in it, e.g. the double Fourier transform of a vertically
sided rectangular distribution of constant amplitude, Amp, and of
dimensions 6xs by 6ys referenced to corner point (xs,yE), i.e.
5x	6 y	5x	6 y
-if(x +-r-) -irKy +-r~) sin(f_—) sin(rj-r~)
F(T , 17) = ^Amp e	e	. 	. 	 , (A. 3 .1.19)
f	V
then a large number of terms including those well outside the basic
region must be retained to faithfully represent C*(xt ,ym,tn) . If
distribution C* (f , rj, 0*) consists of mainly low frequency components,
e.g. the double Fourier transform of a two-dimensional Gausshill
centered at (x5,yE) and of amplitude, Amp, with shaping done by the two
"variance" parameters ox2 and ay2, i.e.
2 2 2 2
o f o T]
-i(fx + rjy ) -(-f— + "f—)
F(f,i7) = Amp e	e	,	(A.3.1.20)
then only a few low frequency terms may be necessary to retain a
faithful representation C*(x, ,ym, tn) , especially if ax2 and ay2 are large
variances.
78

-------
A.3.2 Alternative 2: Biquartic Lagranpe Polynomials
Baptista (1987) gives an extensive discussion of the linear,
quadratic, cubic, quartic, and cubic-Herraite polynomial based methods of
interpolation or estimation of C*(Pi n) for the one-dimensional transport
case, using tables and graphs showing the amplitudes and phases of the
Fourier components. The underlying solution theory for linear parabolic
transport and fate problems readily allows the extension of these one-
dimensional concepts to two and three dimensional linear transport and
fate problems. The one-dimensional quartic (4th degree) Lagrange
interpolation structure presented by Baptista (1987) is extended to a
biquartic (4th degree in each variable) Lagrange interpolation
structure. Note that a 4th degree polynomial is needed in order to be
able to uniquely Lagrange interpolate five distinct points of a
univariate function defined on the interval [x1_2,xi+2] , for i =
2,3,...Nx-2. Therefore, for a two dimensional domain of definition, 25
distinct points are required for a biquartic interpolation formula. The
Lagrange bicanonical basis form is the easiest structure to work with in
this case. Define the biquartic Lagrange approximation to be
C* (P*. )	= E. Z. 2 (x*	y* ,)C* . . .	(A.3.2.1)
ljn	q=l p=l p,q i(n)	¦'j (n) l - 3+p, j - 3+q, n
where
.	, P (x. . . ) P (y. . . )
1 <*T, yt, >> -	Lf"' ¦, ~ „qY, ] ^	(A.3.2.2)
pq t(n)' >j(„)	qy - 3+q
for each i,J combination in the superinterior of the lattice of nodal
points, where superinterior refers to those interior points which lie
beyond at least one line of points inside the boundary line of points.
79

-------
Figure A.4 shows the acceptable region when points are close to the
domain boundaries. Concentrations C*(P*i
-------
Inlet end
ra
$
*55
&


&

re
5
ir
Outlet end
Figure A.4. Two-dimensional lattice of nodal points superimposed on IR2
showing the region of acceptability for the 25 point
interpolation alternative.
81

-------
expected since the one space dimensional version was shown to have no
numerical diffusion (Baptista, 1987). The price paid for the
elimination of numerical diffusion is the occurrence of small
"oscillations" or "wiggles" near the regions of rapid transition or
change in concentration. This numerical problem was expected to occur,
based on the knowledge about the one space dimensional case. The
oscillations or wiggles near a leading front and a trailing edge which
occur in the multidimensional transport case, where multidimensional
Lagrange interpolating polynomials of degree 2 or higher in each space
variable are used, are hypothesized to be caused by the different
amplitude factors and the different propagation velocities of the
multidimensional Fourier harmonic components in the exact solution of
the approximating finite difference (discrete) equations when compared
to the exact solution of the continuous parabolic transport and fate
equation. As was done in alternative 1, the bilinear Lagrange
interpolating polynomials, using the above outlined specialized
hypothesis on the flow domain, it is possible to double Fourier analyze
the approximate distribution C*(x,y,t) at domain point (xf,ym) at time
tn. The approximation for C*(Plj„), analogous to equation (A.3.14), is
now a linear combination of 25 coefficient functions of Crx and Cry, each
multiplying one of the 25 concentrations in the advection molecule
centered at nodal point (xt,ym). Then the form (A.3.1.6) can be
substituted into the analogous expression to (A.3.1.5), and 6 obtained
for the biquartic procedure. Finally, 8 can be expressed in polar form
with the final result being the biquartic approximate distribution
analogue of equation (A.3.1.16). The only changes will be in functions
fi(f»^). definition (A.3,1.9), and f 2 (S"»*7) > definition (A.3.1.10).

-------
These functions will be much more complicated functions of Cr3t, Cj-y, fAx,
and vAy than the bilinear derived forms. We choose not to list these
extremely awkward functions here, and simply rely on the results of the
one-dimensional analysis done by Baptista (1987). However, further work
is planned since in the multidimensional case this is still a largely
unexplored area of research.
Experience gained with the biquartic method shows that to use the
method, the investigator must be very careful in representing his
initial data and his source data. The use of steep sided distributions
or source regions produced wiggles of about 5% the maximum concentration
early on in the above mentioned special hypothesis two-dimensional flow
problem. Using a smooth sided two-dimensional Gausshill initial
distribution, wiggles, whose amplitude grew asymptotically in time to
about 0.2% of the peak or maximum value, were produced. This is a much
more satisfactory situation.
83

-------
TECHNICAL REPORT DATA
{Please read Instructions on the reverse before comple*
1. REPORT NO. 2.
EPA/600/2-91/014
•
4. TITLE AND SUBTITLE
DENITRIFICATION IN N0NH0M0GENE0US LABORATORY SCALE AQUIFERS
4. HYDRAULICS, NITROGEN CHEMISTRY, AMD MICROBIOLOGY IN A
SINGLE LAYER
5. REPORT DATE
7\pril 1991
6. PERFORMING ORGANIZATION CODE
7. AoTHORISI
F.T. Lindstrom, L. Boersma, D. Myrold, and M. Barlaz
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Department of Soil Science
Oregon State University
Corvallis, OR 97331
10. PROGRAM ELEMENT NO.
CBWD1A
11 contract/grant no
CR-814502
12. SPONSORING AGENCY NAME AND AODRESS
U.S. Environmental Protection Agency-
Robert S. Kerr Environmental Research Laboratory
P.O. Box 1198
Ada, OK 74820
13. TYPE OF REPORT AND PERIOD COVERED
Final Report
14. SPONSORING AGENCY CODE
15. SUPPLEMENTARY NOTES
Project Officer: Thomas E. Short FTS: 743-2292
16. ABSTRACT
A two-dimensional mathematical model for simulating the transport and
fate of organic chemicals in a laboratory scale, single layer aquifer is
presented^ The aquifer can be nonhomogeneous and anisotropic with respect
to its fluid flow properties. The physical model has open inlet and outlet
ends and is bounded by impermeable walls on all sides. Fully penetrating
injection and/or extraction wells can be placed anywhere in the flow field.
The inlet and outlet boundaries have user prescribed hydraulic pressure
fields. The steady state hydraulic pressure field is obtained first by
using the two-dimensional Darcy flow law and the continuity equation. The
chemical transport and fate equation is then solved in terms of user
stipulated initial and boundary conditions. The model accounts for the
major physical processes of storage, dispersion, and advection, and also
can account for linear equilibrium sorption, first—order loss processes,
microbial denitrification, irreversible sorption and/or dissolution into
the.organic phase, metabolism in the sorbed state, and first order loss in
the sorbed state. The chemical may be released internally via distributed
leaks, sources that do not perturb the flow field, or fully penetrating
injection wells. Chemical compounds may also enter at the inlet boundary.
Chemical mass balance type inlet and outlet boundary conditions are used.
J7. KEY WORDS ANO DOCUMENT ANALYSIS
a. DESCRIPTORS
b.IDENTIFIERS/OPEN ENOED TERMS
c. COSaTi Field/ Group



IB. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS i This Rtpnrrt
UNCLASSIFIED
21. NO. OF PAGES
91
20. SECURITY CLASS (This pi
UNCLASSIFIED
22. PRICE
EPA Form 2230.1 (R«». 4-77) previous edition is obioliti
i

-------