United States
Environmental Protection
Agency
National Risk Management
Research Laboratory
Ada, OK 74820
Research and Development
EPA/600/SR-97/052
August 1997
&EPA Project Summary
User's Manuals for Two-(2DFATMIC) and
Three-(3DFATM!C) Dimensional
Subsurface Flow, Fate and Transport of
Microbes and Chemicals Models
Gour-Tsyh (George) Yeh, Jing-Ru (Ruth) Cheng, and Thomas E. Short
Two new computer models were
developed for the simulation of
Subsurface Flow, FAte and Transport
of Microbes andChemicals Model using
a Lagrangian-Eulerian adapted zooming
and peak capturing algorithm.
2DFATMIC is a two-dimensional model,
and 3DFATMIC is a three-dimensional
model. The two models are designed to
obtain the density-dependent fluid
velocity field, and to solve the ad vective-
dispersive transport equation coupled
with biodegradatipn and microbial
biomass production. Water flow
through saturated-unsaturated media
and the fate and transport of seven
components (one substrate, two
electron acceptors, one trace element,
and three microbial populations) are
modeled in each program. The input
data include the control indices,
properties of the media either in tabular
or analytical form, the geometry in the
form of elements and nodes, initial
conditions and boundary conditions for
flow and transport, and microbe-
chemical interaction constants.
Principal outputs include the spatial
distribution of pressure head, total head,
moisture content, Darcy velocity
component, concentrations, and
material fluxes at any desired time step.
For each model, fluxes through various
types of boundaries are shown in the
mass balance table. In addition,
diagnostic variables, such as the
number of non-convergent nodes and
residuals, may be printed if desired for
debugging purposes.
Each model can completely eliminate
peak clipping, spurious oscillation, and
numerical diffusion; i.e., solve the
advective transport problems exactly,
within any prescribed error tolerance,
using very large mesh Courant numbers.
The size of the mesh Courant number is
limited only bytheaccuracy requirement
of theEulerian step. Since these models
also include diffusion zooming in
solving diffusion elemental matrix, the
accuracy is improved by specifying the
number of local subelements in every
global element. In other words, the
more subelements zoomed in the
diffusion step, the more accuracy at the
Eulerian step. To sum up, the larger the
time-step size used, the better the
solution that can be obtained with
respect to advective transport. The time-
step sizes are only limited by the
accuracy requirement with respect to
diffusion/dispersion transport and
chemical reaction terms. However, the
limitation of time-step size imposed by
diffusion/dispersion transport is
normally not a very severe restriction.
This Project Summary was developed
by EPA's National Risk Management
Research Laboratory's Subsurface
Protection and Remediation Division,
Ada, OK, to announce key findings of
Printed on Recycled Paper
-------
the research project that is fully
documented in a separate report of the
same title (see Project Report ordering
Information at back).
Introduction
The two new models (2DFATMIC and
3DFATMIC: 2-Dimensional and 3-
Dimensional Subsurface flow, FAte and
Transport of Microbes and Chemicals) can
be used to investigate saturated-
unsaturated flow alone, contaminant
transport alone, combined flow and
transport, or the fate and transport of
microbes and chemicals in a ground-water
environment. The two models are identical
in terms of their background theories and
capabilities except for the mode!
dimensionality. For both models, the
Galerkin finite element method is used to
discretize the Richards' equation for the
flow module, and for the transport module,
the hybrid Lagrangian-Eulerian approach
with an adapted zooming and capturing
algorithm is used to discretize the transport
equation. This approach can completely
eliminate spurious oscillation, numerical
dispersion, and peak clipping due to
advective transport. Large time-step sizes
as well as large spatial-grid sizes can be
used and still yield accurate simulations.
The only limitation on the size of the time
steps is the requirement of accuracy with
respectto dispersive transport, which does
not pose severe restrictions.
The special features of the models are
their flexibility and versatility in modeling
as wide a range of problems as possible.
The models can handle: (1) heterogeneous
and anisotropic media consisting of as
many geologic formations as desired; (2)
both spatially distributed and point sources/
sinks that are spatially and temporally
dependent; (3) the prescribed initial
conditions by input or by simulating a steady
state version of the system under
consideration; (4) the prescribed transient
concentration over Dirichlet nodes; (5) time
dependent fluxes over Neumann nodes;
(6) time dependenttotal fluxes over Cauchy
nodes; (7) variable boundary conditions
for evaporation, infiltration, or seepage on
the soil-air interface for the flow module
and variable boundary conditions of inflow
and outflow for the transport module
automatically; (8) two options for treating
the mass matrix - consistent and lumping;
(9) three options (exact relaxation, under-
and over- relaxation) for estimating the
nonlinear matrix; (10) automatic time-step
size reset when boundary conditions or
sources/sinks change abruptly; (11) two
options, Galerkin weighting or upstream
weighting, for the advection term in the
transport module; (12) two options for the
Lagrangian numerical scheme in the
transport module, enabling and disabling
adapted zooming scheme; (13) two options
for solving the Eulerian step including the
enabling and disabling of diffusion
zooming; (14) mass balance checking over
the entire region for every time step; and,
(15) modification of the program if different
conditions are used.
Description of Models
Mathematical Statement of
2DFATMIC and 3DFATMIC
2DFATMIC and 3DFATMIC are
designed to solve the following system of
governing equations along with initial and
boundary conditions, which describe flow
and transport through saturated-
unsaturated media. The governing
equations for flow, which describe the flow
of a variable-density fluid, are basically
the Richards' equation.
Governing Flow Equation
pw dh ch
Pw
Pw
or —— q
(1)
The saturated hydraulic conductivity Ks
is given by
(P/Pw)
(MO
(2)
where h is the referenced pressure head
defined asp/pwg in whichp is pressure (M/
LT2),KS is the saturated hydraulic
conductivity tensor (L/T), Kr is the relative
hydraulic conductivity or relative
permeability, z is the potential head (L), q
is the source and/or sink (L3/T), and 9 is
the moisture content, p and p. are the
density (M/L3) and dynamic viscosity
(M/LT), and Ksw, pw and nw are the
referenced saturated hydraulic
conductivity tensor, density, and dynamic
viscosity, respectively. The strength of
the source/sink is the discharge or
withdraw flow rate q, and p' is the density
of the injected fluid. These referenced
values are usually taken as the saturated
hydraulic conductivity at zero microbial
and chemical concentrations.
The Darcy velocity is calculated as
follows.
(3)
To solve the governing flow equation,
2
the prescribed initial conditions and
boundary conditions are required. In terms
of initial conditions, the model can take
input initial values or obtain initial values
by simulating a steady state version of the
system under consideration. In terms of
boundary conditions, the model can handle:
(1) the prescribed transient concentration
over Dirichlet nodes; (2) time dependent
fluxes over Neumann nodes; (3) time
dependent total fluxes over Cauchy nodes;
and (4) variable boundary conditions of
evaporation, infiltration, or seepage on the
soil-air interface for the flow module and
variable boundary conditions of inflow and
outflow for the transport module
automatically (see the original reports for
full description).
Governing Equations for Transport
The governing equations for transport
are derived based on the continuity of
mass and flux laws. The major processes
are advection, dispersion/diffusion,
adsorption, decay, source/sink, and
microbial-chemical interactions. Transport
of the substrate, oxygen, nitrate and nutrient
in the bulk pore fluid is expressed by
advection-dispersion equations that are
coupled sink terms that account for
biodegradation. The biodegradation is
modeled with the modified Monod kinetics.
The transport equations for the three
microbes can be derived based on mass-
balance of microbes. The equation for
substrate transport is presented here (for
complete list of equations, see the original
reports).
(4)
pfr cs T cn T
if k> +d*f +Q
-------
where 9 is the moisture content, pb is the
bulk density of the medium (M/L3), t is time
(T), V is the discharge (L/T), V is the del
operator, D is the dispersion coefficient
tensor (L2/T), where As, is the trans-
formation rate constants, Kds, Kd1, Kd2, Kd3
are the distribution coefficient of the
dissolved substrate, microbe No. 1,
microbe No. 2 and microbe No. 3, qin is the
source rate of water, Csh is the source
concentration of substrate, and Cs,C0,Cn,Cp,
are the dissolved concentrations of
substrate, oxygen, nitrogen and nutrient.
l(Cg) = [ 1 + C0/KJ-1 is an inhibition function
which is under the assumption that
denitrifying enzyme inhibition is reversible
and noncompetitive, where K is inhibition
coefficient (M/L3). n0<1>, nn<2>, n0<3> and nn<3>
are the maximum specific oxygen-based
growth rates for microbe No. 1, the
maximum specific nitrate-based growth
rate for microbe No. 2, the maximum
specific oxygen-based growth rate for
microbe No. 3, and the maximum specific
nitrate-based growth rate for microbe No.
3 (1/T), respectively; Y0<1', Yn<2>, Y0<3>, and
Yn(3) are the yield coefficient for microbe
No. 1 utilizing oxygen, the yielding
coefficient for microbe No.2utilizing nitrate,
the yielding coefficient for microbe No. 3
utilizing oxygen and nitrate, in mass of
microbe per unit mass of substrate (MM);
K <1> K <3> K <2> K <3> K <1> K <3> K <2)
'So ' so ' 'Sn ' sn ' *>> ' ^po ' pp '
K <3> are the retarded substrate saturation
constants under aerobic conditions with
respect to microbes No. 1 and No. 3,
respectively, the retarded substrate
saturation constants under anaerobic
conditions with respect to microbes No. 2
and No. 3, respectively, the retarded
nutrient saturation constants under aerobic
conditions with respect to microbes No. 1
and No. 3, respectively, and the retarded
nutrient saturation constants under
anaerobic conditions with respect to
microbes No. 2 and No. 3, respectively;
K0(1), K0<3>, Kn<2>, Kn<3> are the retarded oxygen
saturation constants under aerobic
conditions with respect to microbes No. 1
and No. 3, respectively, and the retarded
nitrate saturation constant under anaerobic
conditions with respect to microbes No. 2
and No. 3 (M/L3), respectively.
The dispersion coefficient tensor D in
Eq. (4) is given by
where |V| is the magnitude of V, 5 is the
Kronecker delta tensor, aT is lateral
dispersivity, aL is the longitudinal
dispersivity, am is the molecular diffusion
coefficient, andt is the tortuosity. To solve
the transport equations, appropriate type
initial and boundary conditions (Dirichlet,
Neumann, Cauchy, and Variable Boundary
Conditions: see the mairi report for detail)
are required.
Numerical Approximation
For the flow module, the Galerkin finite
element method is used to discretize the
Richards' equation which governs the flow
in variably saturated media. For the
transport module, the hybrid Lagrangian-
Eulerian approach with an adapted
zooming and peak capturing algorithm
(LEZOOMPC) is used to discretize the
transport equation. Basically, LEZOOMPC
is a modified method of the zoomable
hiddenfine-mesh approach (LEZOOM) and
exact peak capturing and oscillation-free
scheme (EPCOF) to solve advection-
dispersion transport equations. The
detailed algorithm of the LEZOOMPC can
be found elsewhere. This approach can
completely eliminate spurious oscillation,
numerical dispersion, and peak clipping
due to advective transport. Large time-
step sizes as well as large spatial-grid
sizes can be used and still yield accurate
simulations. The only limitation on the size
of the time steps is the requirement of
accuracy with respect to dispersive
transport which does not pose severe
restrictions.
Figure 1 is the basic conceptual structure
for solving transport by 2DFATMIC and
3DFATMIC. It contains two main steps,
Lagrangian and a Eulerian step. Readers
are referred to the main reports for detailed
discussion. 2DFATMIC consists of a MAIN
program and 77 subroutines whereas
3DFATMIC consists of a MAIN program
and 120 subroutines. The MAIN is utilized
to specify the sizes of all arrays. The
functions of these subroutines are
described in the main reports.
Sample Problems
To verify the two models, several
hypothetical sample problems are
presented in each report. In the 2DFATMIC
user's manual, three illustrative examples
were presented. The first one is a one-
dimensional flow problem through a soil
column. The second one is a one-
dimensional transport problem through a
soil column. The third one is two-
dimensional coupled flow and transport
with biodegradation. In the 3DFATMIC
manual, eight illustrative examples are
presented. The first three examples,
originally designed for 3DFEMWATER, are
flow only problems. The fourth and fifth
examples, originally designed for
3DLEWASTE, are transport only problems.
Example six is a two-dimensional
biodegradation problem which is used to
verify the flow and transport coupling loop
and show the effects of biodegradation.
Examples seven and eight illustrate the
behavior of dissolved organic and oxygen
plumes undergoing natural biodegradation
in a uniform flow field. In this project
summary, two example problems are
introduced; one from 2DFATMIC and
another one from 3DFATMIC program.
Two-Dimensional Coupled Flow
and Multicomponent Transport
Problem for 2DFATMIC program
This example is presented in the
2DFATMIC user's manual as problem
number 3. The example aquifer used for
this example is 1.4 meters long in the X
direction, 1.6 meters thickintheZdirection,
and is shown in Figure 2. This aquifer has
8 different materials, one injection well,
and one extraction well. The initial condition
of the flow field is obtained by simulating a
steady state flow field without sources and
sinks. Then the flowfield and concentration
distribution are updated at each time step.
There are two types of microbes included
in the system, say microbe No. 1 and
microbe No.Swith 1.77xio-4 Kg/m3 initially.
The initial concentrations of the chemicals
are 5* 1O'3 Kg/m3 of substrate, 5* 10'3 Kg/m3
of oxygen, SxlO"3 Kg/m3 of nitrate, and
3* 10-3 Kg/m3 of nutrient. The daily injection
and withdrawal rates of water are 3.75* 10-3
and 3.75x10-3 m3, respectively. The total
hydraulic head is 1.0 m at the upstream
boundary AB and 0.0 m at the downstream
Lagrangian Step
Eulerian Step
Figure 1. The basic structure for coding
transport part of 2DFATMIC and 3DFATM1C.
-------
boundary CD (Figure 2). For transport
simulation, a variable boundary condition
is implemented at the downstream
boundary CD (Figure 2.) and 1.77X10-4
Kg/m3 of microbe No. 1,1.77* 10-" Kg/m3 of
microbe No. 3,1.5*10'2 Kg/m3 of substrate,
5.0*10-3 Kg/m3 of oxygen, 5.0*10'3 Kg/m3
of nitrate, and 3.0>«1 fj-3 Kg/m3 of nutrient
flows through the upstream boundary.
Because microbe No. 2 does not exist in
this environment, the initial and boundary
conditions for this component are set to
zero in this simulation. This problem is set
up for a 4 day simulation.
For numerical simulation the region is
divided into 14x16- 224 square elements
of size - 0.1 by 0.1, and 15 x 17 = 255
nodes. A time-step size of 0.05 is used
and an 80 time-step simulation is made.
Simulation Results
To execute this problem, the maximum
control-integers in the MAIN program of
2DFATMIC must be modified as specified
In the manual. Input parameters used in
this example are shown in tabular format in
the manual. The electronic files containing
the model input and output of this example
are attached in the main report. The results
of velocity fields simulated at t=2 days and
4 days are shown in Figure 3.
Rgure 4 contains the simulation results
of the transport and fate of microbes 1 and
3 at t=2 days (see the main report for other
chemicals).
Three-Dimensional
Multicomponent Transport in a
Uniform Flow Field for
3DFATMIC
This example is presented in the
3DFATMIC user's manual to demonstrate
3-D muiticomponent transport behavior.
The kinetic and microbial parameters for
the simulation are shown in the 3DFATMIC
manual. The region is taken as
0
-------
Microbe 3 at Time = 2 Days
Substrate at Time = 2 Days
1.6-
1.4-
1.2-
1.0-
0.8-
0.6
0.4-1
0.2-
-0.0-
0.0003
i r i i i !
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
1.6-
1.4-
1.2-
1.0-
0.8-
0.6-
0.4-
0.2-
-0.0-
0.0025
0.0050.
'0.0075
0/0100 \
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
X
Figure 4. Simulation results of microbes 1 and 3att = 2 days.
input data, the simulation for these four
components is not performed.
Therefore, the initial and boundary
conditions for these four components
are set to zero in the input data file.
Simulation Results
To execute this problem, the
maximum control-integers in MAIN
program of 3DFATMIC must be modified
as specified in the manual. Input
parameters used in this example are
shown in tabular format in the manual.
The electronic files containing the model
input and output of this example are
attached in the main report. Figure 7
shows the simulation results of total
microbial mass distributions at 200 days
on an x-y cross-section.
Unit: meter
/
f"
n '
/-
— —
> t
/ r
; r
» i
o
-^—
ML
' 'r
f
;0
-4
i
d
P
;V'* c
A B
A: (3,2,2)
B: (5,2,2)
C: (5,3,2)
D: (5,3,4)
Figures. The region of interest.
5
-------
6m
S =3 mg/L
>S^2m
2m
Vx = 0.09 m/d >
0^=0.81 m , Of= 0.005m
Dd = 8.05 x 10"5 mVday, Rs = 1 .4
free - exit
boundary
45m
6m
O = 3.5 mg/L
O=l mg/L
\fNjS3~
2m
\ \ \.
Vx = 0.09 m/d
0^= 0.81m , aj= 0.005m
Dd = 8.05 x 105 mVday, RO = 1-0
= 3.5 mg/L
free - exit
boundary
0
The x-z cross-section of region of interest and the associated physical parameters.
45m
Microbe at Time = 100 Days (NXG=NYG=NZG=2)
5.0-
4.5-
4.0-
3.5-
3.0-
2.5-
2.0-
1.5-
1.0-
0.5-
0.0-
I
10
I
15
I
20
I
25
I
30
1
35
40
. i
45
Microbe at Time = 200 Days (NXG=NYG=NZG=2)
5.0-
4.5-
4.0-
3.5-
3.0-
2.5-
2.0-
1.5-
1.0-
0.5-
0.0-
10
15
I
20
!
25
30
35
40
I
45
X
Figure 7. Total microbial mass distributions at t- 100 and t - 200 days on an x-y cross-section. Concentration isolines are in mg/liter of aquifer
materials.
6
-------
-------
Gour-Tsyh (George) Yeh andJing-Ru (Ruth) Cheng are with the Department of
Civil and Environmental Engineering, Pennsylvania State University, University
Park, PA 16802.
Thomas E. Short is author and EPA Project Officer (see below).
The complete report, entitled "2DFATMIC, Two-Dimensional Subsurface Flow, Fate
and Transport of Microbes and Chemicals Model, User's Manual, Version 1.0"
(Order No. PB97-205637; Cost $31.00, subject to change) and "3DFATMIC,
Three-Dimensional Subsurface Flow, Fate and Transport of Microbes and Chemi-
cals Model, User's Manual, Version 1.0" (OrderNo. PB97-205413; Cost $44.00,
subject to change) will be available only from:
National Technical Information Service
5285 Port Royal Road
Springfield, VA 22161
Telephone: 703-487-4650
The EPA Project Officer can be contacted at:
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
U.S. Environmental Protection Agency
Ada, OK 74820
United States
Environmental Protection Agency
Center for Environmental Research Information
Cincinnati, OH 45268
Official Business
Penally for Private Use $300
BULK RATE
POSTAGE & FEES PAID
EPA
PERMIT NO. G-35
EPA/600/SR-97/052
------- |