United States      Office of Research and    EPA/600/R-97/052
          Environmental Protection   Development       August 1997
          Agency         Washington DC 20460
&EPA  2DFATMIC
          Two-Dimensional
          Subsurface Flow, Fate and
          Transport of Microbes and
          Chemicals Model
          User's Manual
          Version 1.0

-------
  2DFATMIC:  User's Manual of a Two-Dimensional
Subsurface Flow, Fate and Transport of Microbes and
                        Chemicals Model
                            Version 1.0
                                  by
                Gour-Tsyh (George) Yeh and Jing-Ru (Ruth) Cheng
                Department of Civil and Environmental Engineering
                      The Pennsylvania State University
                         University Park, PA 16802

                                 and

                             Thomas E. Short
                    U. S. Environmental Protection Agency
                  Robert S. Kerr Environmental Research Center
                 Subsurface Protection and Remediation Division
                           Ada, Oklahoma 74820
                      Cooperative Agreement CR-818322
                             Project Officer
                             Thomas E. Short
                    U. S. Environmental Protection Agency
                  Robert S. Kerr Environmental Research Center
                 Subsurface Protection and Remediation Division
                           Ada, Oklahoma 74820
                 National Risk Management Research Laboratory
                     Office of Research and Development
                     U.S. Environmental Protection Agency
                           Cincinnati, OH 45268

-------

-------

-------
                                    DISCLAIMER

The U.  S. Environmental Protection Agency through its Office of Research and Development
partially funded and collaborated in the research described here under assistance agreement number
CR-818322 to The Pennsylvania State University. It has been subjected to the Agency's peer and
administrative review and has been approved for publication as an EPA document. Mention of trade
names or commercial products does not constitute endorsement or recommendation for use.

When available, the software described  in this  document is  supplied on "as-is" basis without
guarantee or warranty of any kind, express or implied. Neither the United States Government (United
States Environmental Protection Agency, Robert S. Kerr Environmental Research Center), The
Pennsylvania State University, nor any of the authors accept any liability resulting from use of this
software.
                                           11

-------
                                     FOREWORD

The U.S. Environmental Protection Agency is charged by Congress with protecting the Nation's
land, air, and water resources. Under a mandate of national environmental laws, the Agency strives
to formulate and implement actions leading to a compatible balance between human activities and
the ability of natural systems to support and nurture life. To meet these mandates, EPA's research
program is providing data and technical support for solving environmental problems today and
building a science knowledge base necessary to manage our ecological resources wisely, understand
how pollutants affect our health, and prevent or reduce environmental risks in the future.

The National Risk Management Research Laboratory is the Agency's center for investigation of
technological and management approaches for reducing risks from threats to human health and the
environment. The focus of the Laboratory's research program is on methods for the prevention and
control of pollution to air, land, water, and subsurface resources; protection of water quality in public
water systems; remediation of contaminated sites and ground water, and prevention and control of
indoor air pollution. The goal of this research effort  is to catalyze development and implementation
of innovative,  cost-effective environmental technologies;  develop scientific  and engineering
information needed by EPA to support regulatory and policy decisions; and provide technical support
and information transfer to  ensure  effective implementation of environmental regulations and
strategies.

Bioremediation is unique among remediation technologies in  that it degrades  or transforms
contaminants  through the  use, possibly with  manipulative  enhancement, of  indigenous
microorganisms. Bioremediation can be used in many ways - degradation on concentrated organic
contaminants near their sources, as a secondary remediation strategy following physical or chemical
treatment  methods, for sequestration of metals through microbially  mediated transformation
processes, and for remediating large plumes of dilute contaminants that are broadly dispersed in the
environment.   Thus,  bioremediation  has  the  potential to be one of the  most  cost-effective
technologies for dealing with environmental remediation problems.  Yet, realistically quantitative
predictions and assessments of bioremediation technologies appear lacking. In order to meet the
objectives of having a realistic tool for predicting and assessing if a bioremediation technology can
be successfully implemented, the 2DFATMIC model has been developed.  This numerical model
simulates 1) the fate and transport of multiple microbes, electron acceptors, substrates, and nutrients
and density-dependent fluid flow in saturated-unsaturated subsurface media under either steady-state
or transient conditions; 2) multiple distributed and point sources/sinks as well as boundary sources;
and, 3) processes which  degrade and transform contaminants,  cause  the growth  and death of
microbes, and control the fluid flow.
                                         Clinton W. Hall, Director
                                         Subsurface Protection and Remediation Division
                                         National Risk Management Research Laboratory
                                           in

-------
IV

-------
                                 CONTENTS
LIST OF FIGURES  	     v

LIST OF TABLES	    vii

ABSTRACT   	     ix

1. INTRODUCTION	     1

2. DESCRIPTION OF 2DFATMIC MODEL	     3

      2.1 Mathematical Statement of 2DFATMIC  	     3
      2.2 Numerical Approximation 	     13
      2.3 Description of 2DFATMIC Subroutines	    23

3. ADAPTATION OF 2DFATMIC TO SPECIFIC APPLICATIONS  	    58

      3.1 Parameter Specifications	    58
      3.2 Soil Property Function Specifications	    74
      3.3 Input and Output Devices	    76

4. SAMPLE PROBLEM	    77

      4.1 Problem No. 1 : One-Dimensional Column Flow Problem  	    77
      4.2 Problem No. 2 : One-Dimensional Column Transport	    84
      4.3 Problem No. 3 : Two-Dimensional Coupled Flow and
            Multicomponent Transport Problem	    91
APPENDIX A: Data Input Guide 	    105

REFERENCES  	    131
                                     v

-------
                                  LIST OF FIGURES
2.1  The basic structure for coding transport part of 2DFATMIC  	     21

2.2  Program Structure of 2DFATMIC	     25

2.3  Program Structure of 2DFATMIC (Flow Part)	     33

2.4  Program Structure of 2DFATMIC (Transport Part)  	     40

3.1  Finite Element Discretization and Source/Sink Layout of an Illustrated
       Example for Flow	     62

3.2  Finite Element Discretization and Source/Sink Layout of an Illustrated
      Example for Transport	     66

3.3a An Illustrated Example of Global Nodes and Fine Nodes for Initial Condition ...     70

3.3b An Illustrated Example of Forward-tracked Fine Nodal Points	     71

3.3c An Illustrated Example of Fine Nodal Points and Peak/Valley Points at
       the End of the Time Step for Advection	     73

4.1  Problem definition for the one-dimensional transient flow in a soil column  	     78

4.2  Region of interest for Problem 2 	     85

4.3  The concentration profiles of biomass, substrate, oxygen, and nitrate at
       t = 4 days and t = 8 days	     90

4.4  The region of example 3 with 8 different material properties	     92

4.5  The flow field at (a) t = 2 days and (b) t = 4 days  	     100

4.6  The concentration contours of biomass, substrate, oxygen, nitrate, and nutrient at
       t = 2 days and t = 4 days	     101
                                         VI

-------
                                  LIST OF TABLES
4.1  Input data set for Problem 1  	     82




4.2  The list of input parameters for Problem 1  	     83




4.3  Input data set for Problem 2  	     86




4.4  The list of input parameters for Problem 2	     88




4.5  Input data set for Problem 3  	     93




4.6  The list of input parameters for Problem 3	     96
                                          vn

-------
Vlll

-------
                                         ABSTRACT










       This document is the user's manual of 2DFATMIC, a 2-Dimensional Subsurface Flow,  FAte and




Transport of Microbes and Chemicals  Model using a Lagrangian-Eulerian adapted zooming  and peak




capturing (LEZOOMPC) algorithm.  The 2-dimensional 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 by the accuracy requirement of the Eulerian step.  Since this model also includes 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 is used, the better the solution can




be obtained with respect to advection  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.




       The model, 2DFATMIC, is designed to obtain the density-dependent fluid velocity field, and to solve




the advective-dispersive transport equation coupled with biodegradation and microbial biomass production.




The 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. For each




specific application, 68 maximal control-integers must be assigned using PARAMETER statements in the




MAIN program. In addition, if a user uses different analytical forms of boundary conditions, source/sink




strength value functions, and soil property functions from those used in this program, he is instructed to modify




subroutines ESSFCT, WSSFCT, VBVFCT, DBVFCT, NBVFCT, CBVFCT, and SPFUNC, respectively. The




input data to the program include the control indices, properties of the media either in tabular or analytical






                                               ix

-------
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 output includes the spatial distribution of




pressure head, total head, moisture content, Darcy velocity component, concentrations, and material fluxes at




any desired time step. 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.

-------
                                     1. INTRODUCTION










        2DFATMIC (A TWO-Dimensional Subsurface Flow FAte and Transport of Microbes and Chemicals




Model) 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 groundwater environment. 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 is used to discretize the transport equation. This approach




can completely eliminate  spurious oscillation, numerical  dispersion, and peak clipping due to advection




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 time steps is the requirement of accuracy with respect to




dispersion transport which does not pose severe restrictions.




        The purpose of this manual is to provide guidance to users of the computer code for their specific




application.  Section 2.1 lists the governing equations, initial conditions, and boundary conditions for which




the 2DFATMIC is designed to solve.  Section 2.2 describes the numerical procedure used to simulate the




governing equations. Section 2.3 contains the description of all subroutines in 2DFATMIC.  Since occasions




may  arise that requires the user to modify the code, this section should help the user to understand the code




structure so the user can make necessary adjustments for individual purposes.  Section 3.1 contains the control-




parameter specification.  For each application, the user needs to assign 68 maximal control-integers in the




MAIN program.  Section 3.2 describes the required modification of the code so that one might use a different




analytical form of soil property function from the ones used in this report. Section 3.3 describes files required




for the execution of 2DFATMIC. Appendix A contains the data input guide that is essential for any specific




application.







                                                1

-------
        The user may choose whatever consistent units he wants to.  Units of mass (M), length (L), and time




(T) are indicated in the input description.




        The special features of 2DFATMIC are its flexibility and versatility in modeling as wide a range of




problems as possible.  In terms of media, this model  can handle heterogeneous  and anisotropic media




consisting of as many geologic formations as desired.  In terms of source/sink, this model includes both




distributed and point sources/sinks that are spatially and temporally dependent.  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. In terms of composing the matrix, the model provides two options to treat the mass




matrix:  consistent and lumping.  In terms of linearization, the model gives three options for estimating the




nonlinear matrix: 1) exact relaxation, 2) under-relaxation, and 3) over-relaxation.  To properly treat the




advection terms, the model provides three options:  1) the  Galerkin weighting, 2) the upstream weighting, and




3) the Lagrangian approach. To consider the balance between accuracy and efficiency, the model includes two




options in the Lagrangian  step of solving the advective transport: enabling and disabling adapted zooming




scheme.  It also includes two options in the Eulerian step of solving the diffusion problems: enable and disable




of diffusion zooming.  The model can automatically reset the time step size when  boundary conditions or




sources/sinks have changed abruptly.  It checks the mass balance over the entire region for every time step.




Finally, the source program can easily be modified if microbial processes different from those used in this




model  are used.

-------
                         2. DESCRIPTION OF 2DFATMIC MODEL
2.1 Mathematical Statement of 2DFATMIC





       2DFATMIC is 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





              PI  /9    a fa     dS I 5h             _,    p    ,    p*   ,      p   ,
             -^- a'— + p;6 + n — —  = V-[KsKr-(Vh + -£-Vz)] +-£-q(or  --£-q)       (2.1)

             Pw    ne           dh  5t                   Pw       Pw        Pw
in which the saturated hydraulic conductivity Ks is given by



                                        _     _  (P/PW)
                                                                                         (2.2a)
where  6 is the moisture content, a' is the modified compressiblity of the media (1/L), P' is the modified



compressibility of water (1/L), ne is the effective porosity, h is the referenced pressure head defined as p/pwg



in which p is pressure (M/LT2), t is time (T), K^ is the saturated hydraulic conductivity tensor (L/T), K^ is the



relative hydraulic conductivity or relative permeability, z is the potential head (L), q is the source and/or sink



(L3/T/L), and 6 is the moisture content, p and (i are the density (M/L3) and dynamic viscosity (M/LT) at



microbial concentrations Q, C2, C3, and chemical concentrations Cs, C0, Cn, and Cp (M/L3); Ksw, pw and (\, are



the reference saturated hydraulic conductivity tensor, density, and dynamic viscosity, respectively. The values



of reference saturated hydraulic  conductivity are usually taken as the saturated hydraulic conductivity at zero

-------
microbial and chemical concentrations. The strength of source/sink is the discharge or withdraw flow rate q,

and  p* is the density of injected fluid. It should be noted that in this version of 2DFATMIC, a' and P' were

assumed to be zero, i.e., the media were assumed rigid.

        The density and dynamic viscosity of fluid are functions of microbial and chemical concentrations and

are assumed to take the following form (See Appendix B in 3DFATMIC):



                       —  = 1+E  — -- 1C; ;  i =  1, 2,  3, s,  o, n, p               (2.2b)
                       Pw         i  I  Pw   Pi)
                     f-  =  1 + PSCS + P0C0 - PnCn + PpCp + PA + P2C2 + P3C3               (2.2c)



where Cs and ps are dissolved concentration and intrinsic density of the substrate, respectively (M/L3); C0 and

p0 are dissolved concentration and intrinsic density of oxygen (M/L3), respectively; Cn and pn are dissolved

concentration and intrinsic density of nitrate (M/L3), respectively; Cp and ppare dissolved concentration and

intrinsic density of nutrient (M/L3), respectively; Q and pj are dissolved concentration and intrinsic density

of microbe #1 (M/L3), respectively; C2 and p2 are dissolved concentration and intrinsic density of microbe #2

(M/L3), respectively; C3  and p3 are dissolved concentration and intrinsic density of microbe #3 (M/L3),

respectively; PS, PO, Pro Pp, pb P2, and P3 are factors representing the effect of associated species-concentrations

on the viscosity (L2/T). It is assumed that microbe #1 utilizes the substrate under aerobic conditions, microbe

#2 utilizes the substrate under anaerobic conditions, and microbe #3 utilizes the substrate  under both aerobic

and anaerobic conditions.

       The Darcy velocity is calculated as follows.



                                    V  =   -KK-(—Vh+Vz)                               (2.3)
                                                    P

-------
Initial Conditions for Flow Equation


                                    h = hj(x,z)     in R                               (2.4)


where R is the region of interest and h; is the prescribed initial condition which can be obtained by either field

measurement or by solving the steady state version of Eq. (2.1).




Boundary Conditions for Flow Equation

       Dirichlet Conditions:


                                  h = hd(xb>V)     on Bd                            (2.5)

       Neumann Conditions:


                            -n-KsKr-^Vh =  qn(xb,zb,t)     on Bn                     (2.6)
       Cauchy Conditions:
                        -n-KsKr- (-Vh + Vz) =  qc(xb, zb, t)     on Bc                  (2.7)
                                   P
       Variable Conditions - During Precipitation Period:
                  h =  h (xb,zb,t)     iff -n KKr(-^Vh +Vz) > q     on By           (2.8a)
                                                   P
              or
                   -n-KsKr-(-^Vh+Vz) =  qp(xb,zb,t)   iff  h0     on By            (2.8c)


              or

-------
                   h  = hm(xb,zb,t)   ^n-Ks-Kr-(—Vh+Vz)  < qe      on By            (2.8d)


               or


                  -n-KsKr-(-^Vh+Vz)  = qe(xb,zb,t)  iff  h >  hm      on By           (2.8e)
                             P


where (xb, zb) is the spatial coordinate on the boundary; n is an outward unit vector normal to the boundary;

hd, cjj,, and qc are the prescribed Dirichlet functional values, Neumann fluxes, and Cauchy fluxes, respectively;

Bj, Bn, and Bc are the Dirichlet, Neumann, and Cauchy boundaries, respectively; Bv is the variable boundary;

hp is the allowed ponding depth and qp is the throughfall of precipitation on the variable boundary; hm is the

allowed minimum pressure and qe is the allowed maximum evaporation rate on the variable boundary, which

is the potential evaporation.  Only one of Eqs. (2.8a) through (2.8e) is used at any  point on the variable

boundary at any time.

       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.



Governing Equations for Transport

       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 four nonlinear transport and fate equations are

-------
,dC
            = v-eD-vcs-As(e+PbKds)cs+qincs
                                                     Pw   P
                     Y
                       (i)
                            C
               C
                     Y.
                       (2)
                                             C
                                            pn    p
                                                                (2.9)
tf
Y(3)
0
N(3)
n
v(3)
n
cs
K (3) + r
Kso +(-s
c, 1

^ (3) r
Ksn +(-s
co
0 +^0
c 1

n +^n
P
Kpo +^p
cn
P
K (3) + r
pn ^p
            = v-eD-vco-Ao(e+PbKdo)co+qincom
                     C
         C
                                                                (2.10)
                (3)
                o
                     c

c
5C
            = v-eD-vcn-An(e+PbKdn)cn+qincmn
                                               P      Pw   P
                    (2)
                    sn
                              c
         C
                                                                (2.11)
               .(3)
                    (3)
                    sn
       (3)
                                                     C.
r(3)
1

-------
           ,ac
1
ID-VCp - Ap(6 +PbKdp)Cp + qmCpm + ^V-V( JL) - ^qm Cp
r PW r )
\ n, H^
r^ l/r"
1 ° Y'''
f (2)
C ]L(2)
2] " y(2)
I n
u(3)
(~
° Y(3)
0
C
s
Ks(o1}+Cs
r c
s
K (2) r
Ksn +(-s
C,
s
K(3) /-(
L so +(-s
,,(3)r r T
(3) r-n ^s
+ c
" V(3) 17(3) ,r
Y IV ™r I ,
An [x^sn ^s
C
o
Ki1} + C0
c
n
K (2) + r
n ^n
cn
o
K (3) + r
0 ^0
n
R: (3) + r
n ^n
C
p
KP(;)+CP.
I c 1
p
K (2) + r
Kpn +^p
r c i
p
Kpo +^p
p L
K(3) r n
I
I
I
1

r (2.12)



Co)
pn ^pj





The transport equations for the three microbes can be derived based on mass-balance of microbes and are given
as
   (e+p,,^,)—i+v-vc, =v-eD-vc,-A,


VC -A (9+p,K, )C +q. C • +

(i)
*0

cs
(i)
so +^s
co
(1)
0 +^0
\
P n n*
_^VV(— ) - — q
P Pw P my
' cp '
(1)
po p
X">1
A°
J
                                                                              (2.13)


VC2-A2(e+PbKd2)C2 + qinC2m +
*?
[ cs 1
K(2)+C
[ cn 1
Kn(2)-Cn^
/
, P Pw P '"
' Cp '
Kpn2)+Cp.
A
"I
                                                                              (2.14)

-------
                            =  v-eD-vc3-A3(e+PbKd3)c3+qinc3
                                                                               Pw    P
                 (6+PbKd3)C3-
^
cs
Ks(03)+Cs
co
^-(3) +r
0 o.
f CP
po p
                                                                                           (2.15)
where 6 is the moisture content, pb is the bulk density of the medium (M/L3),  t is time, V is the discharge




(L/T), V is the del operator, D is the dispersion coefficient tensor (L2/T), where As, A0, Ap, An, A1; A2, A3 and




Kds,  Kdo, Kjn, Kdp, Kdl, Kd2, K^ are the transformation rate constants and distribution coefficient of the




dissolved substrate, oxygen, nitrate, nutrient, microbe #1, microbe #2, and microbe #3, Qm is the source rate




of water, Csm Com= Cmn Cm Clm C2mi and C3ll^ are the concentrations of substrate, oxygen, nitrogen, nutrient,




microbe #1, microbe #2 and microbe #3 in the source.
        The dispersion coefficient tensor D in Eq. (2.9) to Eq. (2.15) is given by
                           6D  = aT V 6  + (aL-aT)VV/|V  +  6amT6
(2.16)
where |V| is the magnitude of V, 5 is the Kronecker delta tensor, aT is lateral dispersivity, s^ is the longitudinal




dispersivity, a,,, is the molecular diffusion coefficient, and T is the tortuosity. I(C0) = [ 1 + C0/Kc]~l is an




inhibition function  which is under the assumption that denitrifying enzyme inhibition is reversible and




noncompetitive, where I\, is inhibition coefficient (M/L3).  (i0(1), (V^, Mo(3) and ^3) are the maximum specific




oxygen-based growth rates for microbe #1, the maximum specific nitrate-based growth rate for microbe #2,




the maximum specific oxygen-based growth rate for microbe #3,  and the maximum specific nitrate-based




growth rate for microbe #3 (1/T), respectively; Y0(1), Yn(2), Y0(3), and Yn(3) are the yield coefficient for microbe




#1 utilizing oxygen, the yielding coefficient  for microbe #2 utilizing nitrate, the yielding coefficient for

-------
microbe #3 utilizing oxygen and nitrate, in mass of microbe per unit mass of substrate (M/M); Kso(1), Kso(3),




K,,n(2), Ksn(3), Kp0(1), Kp0(3), Kpn(2), Kpn(3) are the retarded substrate saturation constants under aerobic conditions




with respect to microbe #1, #3, the retarded substrate saturation constants under anaerobic conditions with




respect to microbe #2, #3, the retarded nutrient saturation constants under aerobic conditions with respect to




microbe #1, microbe #3, and the retarded nutrient saturation constants under anaerobic conditions with respect




to microbe #2, microbe #3, respectively; K0(1), K0(3), K^, K^ are the retarded oxygen saturation constants




under aerobic conditions with respect to microbe #1, microbe #3, and the retarded nitrate saturation constant




under anaerobic conditions with respect to microbe #2 and microbe #3 (M/L3), respectively.  A0(1), A0(3), An(2),




and An(3) are the microbial  decay constants of aerobic respiration of microbe #1 and microbe #3, and the




microbial decay constants of anaerobic respiration of microbe #2 and microbe #3 (1/T), respectively; y0(1), y0(3),




yn(2),  and yn(3) are the oxygen-use or nitrate-use  for syntheses by microbe #1, microbe #2, or microbe #3,




respectively; ce0(1), ce0(3), cen(2), and cen(3) are the oxygen-use or nitrate-use coefficient for energy by microbe #1,




microbe #2, or microbe #3, respectively; F0(1), F0(3), Fn(2), and Fn(3) are the oxygen or nitrate saturation constants




for decay with respect to microbe #1, microbe #2,  or microbe #3 (M/L3), respectively; and e0(1), e0(3), en(2), and




en(3) are the nutrient-use coefficients for the production of microbe #1, microbe #2, or microbe #3 with respect




to aerobic or anaerobic respiration, respectively.
                                                   10

-------
Initial Conditions for Transport

                                    Cs=CS1(X>Z)
                                    C0=C01(x,z)
                                    Cn=Cni(x,z)
                                    Cp=cPi(x'z)     in  R                             (2.17)
                                    CI=CH(X,Z)
                                    C2=C2l(x,z)
                                    C3=C3i(x,z)
Prescribed Concentration (Dirichlet) Boundary Conditions
                                  Cs=Csd(xb,zb,t)
                                  C0=Cod(xb,zb,t)
                                  Cn=Cnd(xb,zb,t)
                                  Cp=CPd(xb'zb't)     on Bd                           (2.18)
                                  Ci=Cid(xb.V)
                                  C2=C2d(Xb'Zb't)
                                  C3=C3d(Xb'Zb't)
                                              11

-------
Variable Boundary Conditions





                    n-(VCs-6D-VCs)=n-VCsv(xb,zb,t)




                    n-(VC0-6D-VC0)=n-VCov(xb,zb,t)




                    n-(VCn-6D-VCn)=n-VCnv(xb,zb,t)




                    n-(VCp-6D-VCp)=n-VCpv(xb,zb,t)      // n-V^O             (2.19a)




                    n-(VC1-0D-VC1)=n-VClv(xb,zb,t)




                    n-(VC2-6D-VC2)=n-VC2v(xb,zb,t)




                    n-(VC3-6D-VC3)=n-VC3v(xb,zb,t)




                            n-(-6D-VCs)=0




                            n-(-6D-VCp)=0




                            n-(-6D-VCn)=0




                            n-(-6D-VCp)=0     //n-V>0                     (2.19b)




                            n-(-6D-VC1)=0




                            n-(-6D-VC2)=0




                            n-(-6D-VC3)=0



Cauchy Boundary Conditions





                       n-(VCs-6D-VCs)=qsc(xb,zb,t)




                       n-(VC0-6D-VC0)=qoc(xb,zb,t)




                       n-(VCn-6D-VCn)=qnc(xb,zb,t)




                       n-(VCp-6D-VCp)=qpc(xb,zb,t)     on  BC                 (2.20)




                       n-(VC1-6D-VC1)=qlc(xb,zb,t)




                       n-(VC2-6D-VC2)=q2c(xb,zb,t)




                       n-(VC3-6D-VC3)=q3c(xb,zb,t)






                                        12

-------
Neumann Boundary Conditions
                             n-(-6D-VCs)=qsn(xb,zb,t)
                             n-(-6D-VC0)=qon(xb,zb,t)
                             n-(-6D-VCn)=qnn(xb,zb,t)
                             n-(-6D-VCp)=qpn(xb,zb,t)     on  Bn                       (2.21)
                             n-(-6D-VC1)=qln(xb,zb,t)
                             n-(-6D-VC2)=q2n(xb,zb,t)
                             n-(-6D-VC3)=q3n(xb,zb,t)
where CS1, C01, Cm, Cpl, CH, C2l, and C3l, are the initial concentrations of substrate, oxygen, nitrogen, nutrient,
microbe #1, #2, and #3; R is the region of interest; (xb,zb) is the spatial coordinate on the boundary; n is an
outward unit vector normal to the boundary; Csd, Cod, Cnd, Cpd, Cld, C2d, C3d, and Csv, Cov, Cnv, Cpv, Clv, C2v, C3v,
are the prescribed concentrations of substrate, oxygen, nitrogen, nutrient, microbe #1, #2, and #3, on the
Dirichlet boundary and the specified concentrations of water through the variable boundary, respectively; Bd,
and Bv are the Dirichlet and variable boundaries respectively; q,,, qoc, q^, q^, qlc, cbo, q3C and q^, qon, q^, qpn,
qln, q2n, q3n, are the prescribed total flux and gradient flux of substrate, oxygen, nitrogen, nutrient, microbe #1,
#2, and #3 through  the Cauchy and Neumann boundaries Bc and Bm respectively.

2.2 Numerical Approximation

Flow Equation
       The pressure head can be approximated with Eq. (2.22) by the finite element method.
                                                                                         (2.22)
                                               13

-------
where N is the total number of nodes in the region and Nj and hj are the basis function and the amplitude of


h respectively, at nodal point j.  Substituting Eq. (2.22) into Eq. (2.1) and choosing Galerkin finite element


method, the governing flow equation can be approximated to the following.
                              Pw
FNjdR
                                        dtv
                                        dt
                                              /'(VNi)-K(VNj)dR
                                                                h.
                  I W
                        --2-)qdR- f(VNi)-K-(-P-Vz)dR+ fn-K-(Vh + -H-Vz)N;dB,

                          I W        D          I  W         R           I W
                                where i = l,2,...,N.      F =  —
                                                            dt
                                                                                        (2.23)
Equation (2.23) written in matrix form is:


                                  dh

                                   dt
                                                                                        (2.24)
where {dh/dt} and {h} are the column vectors containing the values of dh/dt and h respectively at all nodes;


[M] is the mass matrix resulting from the storage term;  [S] is stiffness matrix resulting from the action of


conductivity; {G}, {Q} and {B} are the load vectors from the gravity force, internal source/sink, and boundary


conditions, respectively.  The matrices, [M] and [S] are given by
                                          E  K-£-
                                          eMe J     Pw
                                       E /v(NaJK-v(N^R
where
       R, = the region of element e,


       Me = the set of elements that have a local side ce~P coinciding with the global side i-j,


       N»e = the a-th local basis function of element e.
                                                                                        (2.25)
                                                                                        (2.26)
                                              14

-------
Similarly the load vectors {G}, {Q} and {B} are given by
                                           f (VNa>K-

                                           J
                       p                               (2.27)
                       rw
                                                                                       ^    *
Qi =  E  fNae-Pl(or--^)qdR
      "** i    PW      PW
                                                                                       (2.29)
where



       Be = the length of boundary segment e,



       Nse = the set of boundary segments that have a local node a coinciding with the global node i.



The following steps illustrate the incorporation of boundary conditions into matrix equation by finite element



method.



       For the Cauchy boundary condition given by Eq. (2.7), we simply substitute Eq. (2.7)  into Eq. (2.28)



to yield a boundary-element column vector {Bce}  for a Cauchy segment:




                                         1BC6}  =  {qce}                                   (2.30)




where {qce} is the Cauchy boundary flux vector given by
                                                       a =  1  or 2                     (2 31)
The Cauchy boundary flux vector represents the normal fluxes through the two nodal points of the segment



Be on Bc.
                                              15

-------
       For the Neumann boundary condition given by Eq. (2.6), we substitute Eq. (2.6) into Eq. (2.28) to



yield a boundary-element column vector {Bne} for a Neumann segment:





                                         !Bn6} = {qn6}                                  (2.32)



where {qnae} is the Neumann boundary flux vector given by:



                      qnea =  f N>K-^Vz-Naeqn dB ;   a =  1  or 2               (2 33)


                             Bel        Pw          )




which is independent of the pressure head.



       The implementation of the variable-type boundary condition is more involved. During the iteration



of boundary conditions on the variable boundary, one of Eqs.(2.8a) through (2.8e) is used at a node. If either



Eq.(2.8b) or (2.8e) is used, we substitute it into Eq. (2.28) to yield a boundary element column vector {Bve}



for a variable boundary segment:




                                         {Bve} = {qve}                                  (2.34)





where {qve} is the variable boundary flux given by:





                                                      e-^
ae-^qpdB,  or qya =  -JN(Xe-^qedB ;   a  = 1 or 2
                                                                                      (2 35)
Assembling over all Neumann, Cauchy, and variable boundary segments, we obtain the global boundary



column vector {B} as:




                                          IB} =  {q}                                   (2.36)


in which



                           (q)  =  E  {qne) + E (qce} + E (0                    n 3?)
                                  eeN        eeN^        eeN,,                          V '   '
                                              16

-------
where Nne, Nce, and Nve are the number of Neumann boundary segments, Cauchy boundary segments, and


variable boundary segments with flux conditions imposed on them, respectively. The boundary flux {B} given


by Eqs. (2.36) and (2.37) should be added to the right hand side of Eq. (2.24).


       At nodes where Dirichlet boundary conditions are applied, an identity equation is generated for each


node and included in the matrices of Eq. (2.24).  The Dirichlet nodes include the nodes on the Dirichlet


boundary and the nodes on the variable boundary to which either Eq.(2.8a), (2.8c), or (2.8d) is applied.





Transport Equation


       To simplify the notation, the subscript s, o, n, p, 1, 2, or 3 in Eq. (2.9) to (2.15) will be dropped for


the development of numerical algorithm in this section. Since the hybrid Lagrangian-Eulerian approach is used


to simulate Eq. (2.9) to (2.15), it is written in the Lagrangian-Eulerian form as
             (6+pbKd)—=V-(6D-VC)-A(eC+pbKd)+QCin--2-QC+—V-V(-2-)C
                       Dt                                   p       p      Pw         (238)


                     - f(C1,C2,C3,Cs,C0,Cn,Cp)C + g(C1,C2,C3,Cs,C0,Cn,Cp)C
                                                                                       (2.39)
where f(Q, C2, C3, Cs, C0, Cn, Cp) is a microbial-chemical interaction function and g(Q, C2, C3, Cs, C0, Cn, Cp)


is amicrobial growth function . Applying the Galerkin finite element method to Eq. (2.9) through Eq. (2.15),


we obtain
                                       [D] + [K] + [B]){C}  = {Q} + {B}                (2.40)
where {C} is a vector whose components are the concentrations at all nodes, {DC/Dt} is the time derivative


of {C} with respect to time, [M] is the mass matrix associated with the time derivative term, [A] is the stiffness
                                              17

-------
matrix associated with the velocity term which is only computed when steady-state is considered, [D] is the


stifrhess matrix associated with the dispersion term, [K] is the stiffness matrix associated with the decay term


and microbial-chemical interaction, [B] is the stifrhess matrix resulting from boundary conditions, {Q} is the


load vector associated with all source/sink terms,  and {B} is the load vector associated with boundary


conditions.  These matrices  and vectors are given as follows.
                                      eeM
                                           fNae(8+PbKd)NpedR
                                  A,i =  E  fW>-VNRedR
                                                                            (2.4 la)
                                                                                      (2.41b)
                        Dr  =  E
                          1J    eeM,
 EfNj
eeMe J
                                  fN>-V)NpedB + E  fN>-V)NpedB
                                                    R *  J
                 Q, =  E  fNaeQinCmdR+E  fNaeg(C1,C2,C3,Cs,C0,Cn,Cp)dR
                      eefvf  ^               seMa £
                                                                                      (2.4 Ic)
                   A(8 + PbKd) + f(C1,C2,C3,Cs,C0,Cn,Cp) + ^Qm - ^ V-V-P- kedR     (2.4 Id)

                                                          i        i      «
                                                                            (2.41e)
                                                                             (2.4 If)
                B.
                        DD+
                        BeeBv Bp
                    jN>-V)CydB
                                                                                      (2.4 Ig)
where Bv+ is  that part of variable boundary for which the flow is directed into the region, Cv is the


concentration of the incoming fluid through the variable boundary segment Bv+, and Bc, Bn are the Cauchy and


Neumann boundary segments.
                                              18

-------
       The numerical algorithm to solving this partial differential equation is a modified Lagrangian-Eulerian

method with adapted zooming and peak capturing (LEZOOMPC).  Basically, LEZOOMPC is a modified

method of the zoomable hidden fine-mesh approach (LEZOOM) (Yeh, 1990) and exact peak capturing and

oscillation-free scheme (EPCOF) (Yeh et al., 1992) to solve advection-dispersion transport equations.  The

detail algorithm of the LEZOOMPC can be found elsewhere (Cheng, et. al., 1995).  Two of the most important

terminologies used in the LEZOOMPC algorithm are the rough element and smooth element. Two terms need

to be defined, which shall used throughout this document. A smooth element is defined as an element within

which any physical quantity at all points can be interpolated with its node values to within error tolerance. A

rough element is defined as an element within which there exists at least one point for which the physical

quantity cannot be interpolated with its node values to within error tolerance.

       Figure 2.1 is the basic concept structure of solving transport of 2DFATMIC. It contains two main

steps, Lagrangian and Eulerian steps.  First, the concentrations at last time step f1, Cn's, are known quantities

at the beginning of this time step. Second (GNTRAK module), we compute the Lagrangian concentration Q*'s

at global nodes using  the backward node tracking as
                              Q*  =  E  Cj-N/x^Zi*), i = l,2,.,N                        (2.42)
                                      j=i
in which
                                                                                        (2.43)
                                       zVzdt
where N^x/*) is the base function associated with node x^ evaluated at x/*; Vx and Vz are the velocities along

x- and y- direction, respectively.
                                              19

-------
        Third (HPTRAK modules), all the activated fine grid nodes and the global nodes are forwardly tracked

to obtain the Lagrangian concentration Cjf by the following equations:

                   Cjf = C(x/,z/,tn+1)  = C(Wn)  =Cjn     j=l,2,	,N+Nn            (2.44)

in which
                                        tn+i                    Vi
                             Xjf = Xj + I Vxdt  ,   Zjf =  Zj + I Vzdt                       (2.45)
                                        tn                     tn

It should be noted that C/'s are exact if C^'s are exact and Eq.(2.45) is integrated exactly.
                                                20

-------
Lag rang ian Step
    Eulerian Step

>
f
GNTRAK
>
f
HPTRAK
>
f
SFDET
>
f
FGDET
>
f
ISEHIL




v
(cs )

>
1
DFPREP
>
1
ASEMBL
>
f
SOLVE






•w
              Figure 2.1 The basic structure for coding transport part of 2DFATMIC
                                          21

-------
       Fourth (SFDET and FGDET modules), we determine if an element is a rough element in SFDET

module (Yeh, 1990) based on the prescribed error tolerance. The criterion is shown as follows:
                                        •> a _ p f\,„ f
                                        "j   UJ /^J

-------
which is based on Eqs. (2.9) to (2.15), Where [Ae] is the element coefficient matrix, {Ce} is the unknown




vector of the concentration, and {Re} is the element load vector. Element e can be a global element or a




subelement generated in a rough region by DFPREP.  Then, this module assembles all the element matrix




equations to a global matrix equation, which will be solved by a diffusion solver.




       The module SOLVE solves the assembled global matrix equation by a direct solver or pointwise




iteration method.  If the diffusion zoomed approach is activated in the Eulerian step, the latter is selected




forcefully.




       At the very end of this time step, i.e. at tn+1, the concentrations at all fine nodal points (including




regularly refined nodes and peak/valley nodes) generated in Lagrangian step, are obtained by finite element




interpolations as follows:
                                                    k e A11 Fine  Nodal  Points
2.3 Description of 2DFATMIC Subroutines







       2DFATMIC consists of a MAIN program and 77 subroutines. The MAIN is utilized to specify the




sizes of all arrays. The control and coordinate activity are performed by the subroutine HTMICH.  Figure 2.2




shows the structure of the program. The functions of these subroutines are described below.




Subroutine HTMICH




       The subroutine HTMICH controls the entire sequence of operations, a function generally performed




by the MAIN program. It is, however, preferable to keep a short MAIN and several subroutines with variable




short allocation. This makes it possible to place most of the FORTRAN deck on a permanent file and to deal




with a specific problem without making changes in array  dimensions throughout all subroutines.




       The flow of data input for the model is also anchored by the HTMICH. The subroutine RDATIO is




called to read the geometric and material data.  HTMICH then calls subroutine IBE2D to generate IBE array





                                              23

-------
which stores the index of boundary elements.  The source/sink data and boundary conditions for flow and




transport are read by subroutines FDATIO and TDATIO, respectively.  The microbial-chemical reaction




parameters are also read in subroutine TDATIO. Figure 2.2 shows the flow chart of this subroutine.




       Depending on the combinations of the parameters KSSh, KSSt, NTI, and IMOD, the subroutine




HTMICH will perform either the steady-state flow and/or transport computations only, or the transient state




flow and/or transport computations using the flow and/or transport steady-state solution as the  initial




conditions, or the transient flow and/or transport computation using user-supplied initial conditions.




       HTMICH calls the subroutines ESSFCT, WSSFCT, CBVFCT, NBVFCT,  VBVFCT, and DBVFCT




to obtain sources/sinks and boundary values; subroutine SPROP to obtain the relative hydraulic conductivity,




water capacity, and moisture content from the pressure head;  subroutine VELT to compute Darcy's velocity;




subroutine FSFLOW to calculate flux through all types of boundaries and water accumulated in the media;




subroutine FPRINT to print out the results; and subroutine FSTORE to store the flow variables for plotting;




subroutine HMCHYD to perform the  flow computations;  subroutine FLUX to compute  material flux;




subroutine AFABTA to obtain upstream weighting factors based on velocity and dispersivity; subroutine




TSFLOW to calculate material fluxes through all types of boundaries and the amount of water accumulated




in the media; subroutine TPRINT to print out the transport computation results; subroutine TSTORE to store




the transport computation results for plotting;  subroutine THNODE to compute  the value of the moisture




content plus bulk density times distribution  coefficient; subroutine  DISPC to compute the dispersion




coefficients; subroutine HMCTRN to perform the transport computations.
                                              24

-------
Figure 2.2  Program structure of 2DFATMIC (MAIN)
                     25

-------
Subroutine RDATIO




       Subroutine RDATIO reads data set 2 to data set 6 described by Appendix A. It also prints all the input




information, and calls subroutine READR and READN, respectively, to read and/or automatically generate




real and integer numbers; subroutine LNDGEN to generate two pointer arrays (one is the node stencil and the




other is the  connecting elements to any node); subroutine SURF to identify the boundary segments and




boundary nodes.




Subroutine READR




This subroutine is called by the subroutines RDATIO, FDATIO, TDATIO, and HTMICH to read some real




numbers and automatically generate other real numbers if required. Automatic generation of regular patterned




data is built into the subroutine (see Appendix A).




Subroutine READN




This subroutine is also called by the subroutines RDATIO, FDATIO, and TDATIO to generate read some




integers and automatically generate other integers if required (see Appendix A).




Subroutine LNDGEN




This subroutine is called by the controlling subroutine HMCTRN to preprocess the two pointer arrays.  One




is the global node connectively (stencil) LRN(J,N), where LRN(J,N) is the global node number of J-th node




connected to the global node N. This pointer array serves two purposes: (1) to assemble the global matrix in




compressed form when the pointwise iteration method is used, (2) to speed up locating fictitious particles in




the Lagrangian step.  The other pointer array is the global elements connected to nodes LRL(K,N), where




LRL(K,N) is the global element number of the K-th elements connected to the global node N.  These two




pointer arrays are generated based on the element connectivity IE(M,J). Here, IE(M,J) is the global node




number of J-th node of element M. These arrays are used in particle tracking, i.e., subroutine GNTRAK and




HPTRAK.
                                              26

-------
Subroutine SURF




       Subroutine SURF identifies the boundary sides, sequences the boundary nodes, and computes the




directional cosine of the surface sides.  The mapping from boundary nodes to global nodes are stored in




NPBB(I) (where NPBB(I) is the global node number of i-th boundary node). The element  number associated




with the boundary sides are stored in ISB array. The length and directional cosines for each side are stored




in DLCSB array, respectively.  The local and global nodal number of two nodes of each side are also stored




in ISB. The information contained in NPBB, ISB, and DLCSB along with the number of boundary nodes and




the number of boundary sides are returned to subroutine RDATIO for other users.




Subroutine IBE2D




       The  subroutine IBE2D is used to generate the index of boundary element stored in IBE array. If




IBE(M) = 0, it means no boundary element side in the M-th element.  If IBE(M)=12, the element side




connected by local node number 1 and 2 is the boundary element side of the M-th element.




Subroutine FDATIO




       The  subroutine FDATIO is called by the program HTMICH and to control the  input of initial




condition, sources/sinks, and boundary conditions, in time and space, assigned to each specified node/element




for flow simulations. Users need to give the initial values, source/sink profiles, and boundary profiles, to




specify the global node/element numbers of the boundary,  and to assign boundary profile type to each




node/element for flow simulations. The source/sink type for each node/element  is also assigned in  this




subroutine according to the data given by the users.




Subroutine TDATIO




       The  subroutine TDATIO is called by the program HTMICH and to control the  input of initial




condition, sources/sinks, and boundary conditions, in time and space, assigned to each specified node/element




for transport simulations.  Users need to give the initial values, source/sink profiles, and boundary profiles, to




specify the global node/element numbers of the boundary, and to assign boundary profile type to each
                                              27

-------
node/element for transport simulations. The source/sink type for each node/element is also assigned in this




subroutine according to the data given by the users.




Subroutine ESSFCT




        This subroutine is called by the subroutine HTMICH to compute the elemental source strength. It uses




linear interpolation of the tabular data or it computes the value with an analytical function. If the latter option




is used, the user must supply the function to this subroutine.




Subroutine WSSFCT




        This subroutine is called by the subroutine HTMICH to compute the well source strength. It uses the




linear interpolation of the tabular data or it computes the value with analytical function. If the latter option is




used, the user must supply the function  into this subroutine.




Subroutine VBVFCT




        This subroutine is called by the subroutine HTMICH to compute the variable boundary values. It uses




linear interpolation of the tabular data or it computes the value with an analytical function. If the latter option




is used, the user must supply the function to this subroutine.




Subroutine DBVFCT




        This subroutine is called by the subroutine HTMICH to compute the Dirichlet boundary values.  It




uses linear interpolation of the tabular data or it computes the value with an analytical function.  If the latter




option is used, the user must supply the function to this subroutine.




Subroutine CBVFCT




        This subroutine is called by the subroutine HTMICH to compute the Cauchy fluxes.  It uses linear




interpolation of the tabular data or it computes the value with an analytical function. If the latter option is used,




the user must supply the function to this subroutine.




Subroutine NBVFCT




        This subroutine is called by the subroutine HTMICH to compute the Neumann fluxes.  It uses linear
                                               28

-------
interpolation of the tabular data or it computes the value with an analytical function. If the latter option is used,

the user must supply the function to this subroutine.

Subroutine FSFLOW

       This subroutine is used to compute the fluxes through various types of boundaries and the increasing

rate of water content in the region of interest.  The function of FRATE(7) is to store the flux through the whole

boundary enclosing the region of interest. It is given by
                              FRATE(7) = J(Vxnx  +  Vznz)dB  ,                       (2.49)
where B is the global boundary of the region of interest; Vx and Vz are Darcy's velocity components; and nx

and nz are the directional cosines of the outward unit vector normal to the boundary B. FRATE(l) through

FRATE(5) store the flux through Dirichlet boundary BD, Cauchy boundary Bc, Neumann boundary BN, the

seepage/evapotranspiration boundary Bs, and infiltration boundary Bn respectively, and are given by
                             FRATE(l) = J(Vxnx  + Vznz)dB ,                      (2.50a)
                             FRATE(2) = f(V^x  + Vznz)dB ,                      (2.50b)
                                           B.,
                             FRATE(3) = J(Vxnx  + Vznz)dB ,                      (2.50c)
                             FRATE(4) = J(Vxnx  + Vznz)dB ,                      (2.50d)
                              FRATE(5) = J(Vxnx  + Vznz)dB ,                      (2.5oe)
                                             29

-------
       FRATE(6), which is related to the numerical loss, is given by
                                                     5

                         FRATE(6)  =  FRATE(7) - £ FRATE(I)                   (2.51)

                                                    1=1
       FRATE(8) and FRATE(9) are used to store the source/sink and increased rate of water within the



media, respectively:
                                 FRATE(8) = -/-C-qdR ,                           (2.52)




and



                                   „-,_,»>.     e p d9 dh ,_
                               FRATE(9)  =  /-£-       dR                          (2.53)

                                            JRpwdhat





where R is the region of interest. If there is no numerical error in the computation, the following equation



should be satisfied:




                         FRATE(9) =  -[FRATE(7) + FRATE(8)]                   (2.54)




and FRATE(6) should be equal to zero. Equation (2.54) simply states that the negative rate of water going



out from the region through the entire boundary and due to  a source/sink is equal to the rate of water



accumulated in the region.



Subroutine 04TH and 03TH



       These subroutines are used to compute the contribution of the increasing rate of the water content from



an element e
                                 QTHP = J-^__dR,                           (2.55)


                                          T? f"W
                                            30

-------
The computation of the above integration is straightforward.  Q4TH is used to evaluate the integral for




quadrilateral elements and Q3TH for triangular elements.




Subroutine FPRINT




       This subroutine is used to line-print the flow variables. These include the fluxes through variable




boundary surfaces, the pressure head, total head, moisture content, and Darcy's velocity components.




Subroutine FSTORE




       This subroutine is used to store the flow variables on Logical Unit  11. It is intended to be used for




plotting. The information stored includes region geometry, subregion data,  and hydrological variables such




as pressure head, total head, moisture content, and Darcy's velocity components.




Subroutine TPRINT




       This subroutine is used to line-print the simulation results of the contaminant transport. These include




the material flux components and the concentration at each node.




Subroutine TSTORE




       This subroutine is used to store the simulation results on Logical Unit 12. It is intended for plotting




purposes.  The information stored includes region geometry, concentrations, and material flux components at




all nodes for any desired time step.




Subroutine ADVWRK




       This subroutine is called by HTMICH to prepare all the working arrays to be further used in ELETRK.




This subroutine is called only when the transient simulation is required.




Subroutine HMCHYD




       HMCHYD calls subroutine SPROP to obtain the relative hydraulic conductivity, water capacity, and




moisture content from the pressure head; subroutine VELT to compute Darcy's velocity; subroutine BCPREP




to determine if a change of boundary conditions is required; subroutine FASEMB to assemble the element
                                               31

-------
matrices over all elements; subroutine FBC to implement the boundary conditions; subroutine BANSOL or




Subroutine PISS to solve the linear matrix equations. Figure 2.3 shows the flow chart of this subroutine.
                                             32

-------
HMCHYD
           SPROP
           BCPREP
           FASEMB
             FBC
           BANSOL
       OR
             PISS
            VELT
SPFUNC
                          FQ3
                          FQ4
F2CNVB
                          FQ3D
 FQ4D
                        BANSOL
             BASE
BASE
          Figure 2.3 Program Structure of 2DFATMIC (Flow Part)
                       33

-------
Subroutine SPROP




       This subroutine calculates the values of moisture content, relative hydraulic conductivity, and the water




capacity.   This subroutine calls subroutine SPFUNC to calculate the soil property function by either tabular




input or analytical functions.




Subroutine SPFUNC




       This subroutine calculates the soil property function by either tabular input or analytical functions.




When analytical functions are used, the users must supply the functional form and modify this subroutine.




Subroutine BCPREP




       This subroutine is called by HMCHYD to prepare the infiltration-seepage boundary conditions during




a rainfall period or the seepage-evapotranspiration boundary conditions during non-rainfall periods.  It decides




the number of nodal points on the variable boundary to be considered as Dirichlet or Cauchy points.  It




computes  the number of points that change boundary conditions from ponding depth (Dirichlet types) to




infiltration (Cauchy types), or from infiltration to ponding depth, or from minimum pressure (Dirichlet types)




to infiltration during rainfall periods. It also computes the number of points that change boundary conditions




from potential evapotranspiration (Cauchy types) to minimum pressure, or from ponding depth to potential




evapotranspiration, or from minimum pressure to potential evapotranspiration during non-rainfall periods.




Upon completion, this subroutine returns the Darcy flux, infiltration/potential evapotranspiration rate, the




ponding depth nodal index, the flux-type nodal index, the minimum pressure nodal index, which are stored




in RSVAB array, and the number of nodal points (NCHG) that have changed boundary conditions.




Subroutine FASEMB




       This subroutine calls FQ4 or FQ3 to evaluate the element matrices.  It then sums over all element




matrices to form a global matrix equation governing the referenced pressure head at all nodes.




Subroutine FQ4 and FQ3




       These subroutines are called by the subroutine FASEMB to compute the element matrix given by
                                               34

-------
                                  QA(I,J)  =  NFN/dR ,                           (2.56a)
                              QB(I,J)  = J(VN1e)-KKr-(VN/)dR ,                      (2.56b)
where F is the soil property function.  Subroutine FQ4 and FQ3 also calculate the element load vector given

by
                   RQ(I) = J[(VN1e)-KsKr-(-E-Vz)  - Nle-(or--P-)       ,            (2.56c)
                                             PW           PW      PW
where q is the source/sink. Subroutine FQ4 is to compute Eq.(2.56a) to Eq.(2.56c) for quadrilateral elements,

while FQ3 is to compute those for triangular elements.

Subroutine BASE

        This subroutine is called by  subroutines FQ4D, Q4TH and FQ4 to evaluate the value of the base

function at a Gaussian point.  The computation is straightforward.

Subroutine FBC

        This subroutine incorporates Dirichlet, Cauchy, Neumann, and variable boundary conditions.  For a

Dirichlet boundary condition, an identity algebraic equation is generated for each Dirichlet nodal point. Any

other equation having this nodal variable is modified accordingly to simplify the computation.  For a Cauchy

boundary, the integration of Eq. (2.31) is obtained in subroutine FBC, and the result is added to the load

vector.  For a Neumann boundary, the integrations of both the Neumann and gravity fluxes (Eq.(2.33)) are

added to the load vector. The subroutine FBC also implements the variable boundary conditions.  First, it

checks all infiltration-evapotranspiration-seepage points, identifying any of them that are Dirichlet points. If

there are Dirichlet points, the method of incorporating Dirichlet boundary conditions mentioned above is used.

If a given point is not the Dirichlet point, the point is bypassed. Second, it checks  all rainfall-evaporation-

                                              35

-------
seepage points again to see if any of them is a Cauchy point.  If it is a Cauchy point, then the computed flux

by infiltration or potential evapotranspiration (Eq.(2.35)) is added to the load vector. If a given point is not

a Cauchy point, it is bypassed. Because the infiltration-evaporation-seepage points are either Dirichlet or

Cauchy points, all points are taken care of in this manner.

Subroutine F2CNVB

       This subroutine is called by the subroutines FBC to compute the surface node flux of the type

                                    RQ(I)  = -  NXqvdB
                                                                                        (2
                                               B     w
                                    RQ(I)  = -   N
or
                             RQ(I)  =
                                                                                        (2.57c)
where qv, qc, and qn are the rainfall/evaporation, Cauchy, or Neumann fluxes.

Subroutine BANSOL

       This subroutine is called by the subroutine VELT and HMCHYD to solve for the matrix equation of

the type:

                                         [C]{x}  =  {y}                                  (2.58)


where [C] is the coefficient matrix and {x} and {y} are two vectors, {x} is the unknown to be solved and {y}

is the known load vector.  The computer returns the solution {y} and stores it in {y}. The computation is a

standard banded Gaussian direct elimination procedure.
                                               36

-------
Subroutine PISS

       This subroutine is called by either HMCHYD or HMCTRN which is determined by identification

integer, IFT, to solve for the matrix equation of the type

                                          [C](x} =  {y}                                   (2.59)

where [C] is the coefficient matrix and {x} and {y} are two vectors, {x} is the unknown to be solved, and {y}

is the known load vector. The computer returns the solution {x} and stores it in {x}.  The computation is a

pointwise iteration method. Three options are provided for the point iteration solution strategy, the under-

relaxation, Gauss-Seidel, and over-relaxation.

Subroutine VELT

       This subroutine calls FQ4D and FQ3D to evaluate the element matrices and the  derivatives of the total

head. It then sums over all element matrices to form a matrix equation governing the velocity components at

all nodal points. To save computational time, the matrix is  diagonalized by lumping. The velocity components

can thus be solved point by point.  The computed velocity field is then returned to HTMICH or HMCHYD

through the argument. This velocity field is also passed  to subroutine BCPREP to evaluate the Darcy flux

across  the seepage-infiltration-evapotranspiration surfaces, and to  subroutine  HMCTRN to perform the

transport computation.

Subroutine F04D and F03D

       Subroutine FQ4D and FQ3D are called by the subroutine VELT to compute the element matrices

given by


                                    QB(I,J) = jN^/dR                             (2.60a)
                                               Be

where N;e and Nje are the base functions for nodal point i and j of element e, respectively. Subroutine FQ4D

and FQ3D also evaluate the element load vector:

                                              37

-------
                  QRX(I) =  -N^-K^.-CVN^hjdR  -  N^-K^-VzdR          (2.60b)
                               "            D               "
                               R           P               R
                  QRZ(I)  = -N^k-KKCVNj^LdR-N^k-KK^VzdR           (2.60c)
                               Re            P              Re

where

      Hj = the total head at nodal point j,

      i = the unit vector along the x-coordinate,

      k = the unit vector along the z-coordinate,

      Ks = the saturated hydraulic conductivity tensor,

      Kj = the relative hydraulic  conductivity.

FQ4D is used to evaluate the integration for quadrilateral elements and FQ3D for triangular elements.

Subroutine HMCTRN

       The subroutine HMCTRN controls the entire sequence of transport computations. HMCTRN calls

subroutine THNODE to compute the value of moisture content plus bulk density times distribution coefficient

in the case of linear isotherm; subroutine DISPC to calculate the dispersion coefficient (Eq. (2.16)) associated

with each Gaussian point in every element; subroutine AFABTA to obtain upstream weighting factors based

on velocity and dispersivity; subroutine TASEMB to assemble  the element matrices over all elements;

subroutine TBC to implement the boundary conditions; subroutine TBC1 to apply intra-boundary conditions

with the implementation  of the slave point concept to circumvent the incompatibility problem; subroutine

SOLVE or subroutine PISS to solve the resulting matrix equations with the direct Gaussian elimination method

or pointwise iteration methods; subroutine FLUX to compute material flux; subroutine TSFLOW to calculate

the fluxes through all types of boundaries and accumulated in the media; subroutine TPRINT to print out the

results; subroutine GNTRAK to perform particle tracking of global nodes; subroutine HPTRAK to perform

particle tracking of fine nodal points, and subroutine ADVBC to implement boundary conditions in the


                                              38

-------
Lagrangian step; subroutine  SFDET to determine sharp front elements;  subroutine FGDET to  imbed




subelements into every sharp front element; subroutine ISEFflL to prepare ISE array which stores the indices




of subelements and to determine the activation of the points with the highest or lowest concentrations in each




subelement; subroutine DFPREP to prepare the fine mesh nodes and elements for diffusion zooming.  Figure




2.4 shows the flow chart of this subroutine.




Subroutine THNODE




       This subroutine is called by HTMICH and HMCTRN to compute (6 +pbKd) for the linear isotherm




model.




Subroutine DISPC




       Subroutine DISPC calculates the dispersion coefficient associated with each Gaussian point in every




element.  The three components of the dispersion coefficient tensor D are calculated using Eq. (2.16) and




stored in AKDC array.




Subroutine AFABTA




       This subroutine calculates the values of the upstream weighting factors along 4 or 3 sides of all




elements.




Subroutine TASEMB




       This subroutine calls TQ4 and TQ3 to evaluate the element matrices. It then sums over all element




matrices to form a global matrix equation governing the concentration distribution at all nodes.




Subroutine TQ4 and Subroutine TQ3




       These subroutines are called by the subroutine TASEMB to compute the element matrix given by
                                    QAaJ)=/Nie0NjedR ,                            (2.61a)
                                              39

-------
WARMSG! CENTER!  GRID  11 BASEI  IELENO
    WARMSq | CKTRI  |NEWTRl| [ONLINE
      Figure 2.4 Program structure of 2DFATMIC (Transport Part 1 of 2)




                                 40

-------
ALGBDY—  FIXCHK

REPLAS

ONLINE

FCOS
LOCQ2D

FCOS

LOCQ2D

>

REPLAS

VALBDL

MMLOC

CKCNEL

WRKARY

TRACK2

IBWCHK

CKCOIN

TRACK 1

ONLINE

CKCNEL
—

—
               Figure 2.4  Program structure of 2DFATMIC (Transport Part 2 of 2)
                                          41

-------
                                QAAaJ)=/NiepbKdNjedR ,                         (2.61b)
                                          Re

                              QB(I,J)=J(VN1e)-6D-(VNJe)dR ,                       (2.61c)
                                      Re

                             QBB(I,J) =
                                QVaJ)=/NieV-(VNje)dR ,                         (2.61e)
                          qc(i,j)= fN^--Q --v-v-^N/dR ,                   (2.6if)
                                  Re     '      '      ' w

                  QDftJ) = ±|Nj (6 +PbKd) f(C1,C2,C3,Cs,C0,Cn,Cp) NjdR            g
                             Re

where f(Cb C2, C3, Cs, C0, Cn, Cp) should be evaluated using those seven dissolved concentrations at previous

iteration. Subroutines FQ3 and FQ4 also calculate the element load vector given by:
                                   QR(I)=jN1eQCmdR ,                            (2.61h)
Subroutine SHAPE

       This subroutine is called by subroutines TQ4D, TQ4 and Q4Rto evaluate the value of the base and

weighting functions and their derivatives at a Gaussian point. The computation is straightforward.

Subroutine BASE1

       This subroutine is called by INTERP, INTERPK, ADVRX, SFDET, TASEMB, and HMCTRN to

compute the base function values associated with a specified point based on the given two-dimensional global

coordinates. For the cases of quadrilateral elements, it calls XSI2D to calculate the local coordinates, and

                                            42

-------
computes base functions with these determined local coordinates.  For the cases of triangular elements, the




base functions can be analytically determined based on the given global coordinates.




Subroutine XSI2D




       This subroutine is called by BASE1 to compute the local coordinate of an element given the global




coordinate within that element. With the local coordinate, the base functions can then easily be computed by




formula.




Subroutine TBC




       This  subroutine incorporates Dirichlet, variable boundary, Cauchy,  and Neumann boundary




conditions. For a Dirichlet boundary condition, an identity is generated for each Dirichlet nodal point.  Any




other equations having this nodal variable are modified accordingly to simplify the computation.  For a variable




surface, the integration of the normal velocity times the incoming concentration is added to the load vector and




the integration of normal velocity is added to the matrix. For the Cauchy boundaries, the integration of Cauchy




flux is added to the load vector and the integration of normal velocity is added to the matrix.  For the Neumann




boundary, the integration of gradient flux is added to the load vector.




Subroutine Q2CNVB




       These subroutines are called by the subroutines TBC to compute the surface node flux of the type
                                       RQO)=|NieqdB ,                                (2.62a)
where q is either the Cauchy flux, Neumann flux, or n VCV. It also computes the boundary element matrices
                                    BQaJ)=/NieVNjedR ,                            (2.62b)
                                               43

-------
Subroutine TBC1




       This subroutine is called whenever the total number of nodes for composing the matrix is greater than




the total number of global nodes, i.e., the diffusion zooming scheme is employed.  The "slave point" concept




takes care of the incompatibility for the intraboundary points between rough and smooth regions.  This




subroutine implements the concept so that the entries for the intraboundary points of the matrix equation can




be modified. If there are diffusion fine grids falling on the global boundaries, the  "slave point" concept also




resolves the problems of implementation of boundary conditions for these fine grids. Subroutine SOLVE




       This subroutine  is called by HMCTRN to solve a matrix equation of the type





                                         [C]{x}={y}                                  (2.63)








where [C] is the coefficient matrix and {x} and {y} are two vectors, {x} is the unknown to be solved and {y}




is the known load vector. The computer returns the solution {y} and stores it in {y}. SOLVE uses a standard




banded Gaussian direct-elimination procedure.




Subroutine FLUX




       This subroutine calls TQ3D and TQ4D to evaluate the element matrices and the derivatives of




concentrations.  It  then sums over all  element matrices to form a matrix equation governing the flux




components at all nodal points.  To save computational time, the matrix is diagonalized by lumping. The flux




components due to dispersion can thus be solved point by point. The flux due to the velocity is then added




to the computed flux due to dispersion.  The computed total  flux field is then returned to HTMICH and




HMCTRN through the argument.




Subroutine TQ4D and TQ3D




       Subroutine  TQ4D and TQ3D  are called by the  subroutine FLUX to compute the element matrices




given by
                                             44

-------
                                    QBftJ)=|NieNjedR ,                             (2.64a)
where N;e and Nje are the basis functions for nodal point i and j of element  e, respectively which are




computed by subroutine SHAPE in TQ4D. Subroutine TQ4D also evaluates the element load vector:
                                         |Niei-eD-(VNje)CjdR ,                       (2.64b)
                                        |Niek-0D-(VNje)CjdR ,                       (2.64c)
where Cj is the concentration at nodal point j, i is the unit vector along the x-direction, k is the unit vector




along the z-coordinate, 6 is the moisture content, and D is the dispersion coefficient tensor. Subroutine TQ4D




is used for quadrilateral elements and TQ3D is used for triangular elements.




Subroutine TSFLOW




       This subroutine is used to compute the fluxes  through various types of  boundaries and the




accumulation rate of material in the region of interest.  FRATE(7) is to store the flux through the whole




boundary
                                FRATE(7)=J(Fxnx+Fznz)dB ,                          (2.65)
where B is the global boundary of the region of interest; Fx, and Fz are the flux components; and nx, and nz




are the directional cosines of the outward unit vector normal to the boundary B.  FRATE(l) stores the fluxes




through Dirichlet boundary Bd. FRATE(2) and FRATE(3) store the fluxes through Cauchy and Neumann
                                             45

-------
boundaries, respectively. FRATE(4) and FRATE(5) store incoming flux and outgoing flux rates, respectively,


through the variable boundaries Bv~ and Bv+, as given by
                              FRATE(l)=J(Fxnx+Fznz)dB  ,                       (2.66a)
                              FRATE(2)=J(Fxnx+Fznz)dB  ,                       (2.66b)
                              FRATE(3)=J(Fxnx+Fznz)dB  ,                       (2.66c)
                              FRATE(4)=J(Fxnx+Fznz)dB,                       (2 66d)
                                         B _
                              FRATE(5)=J(Fxnx+Fznz)dB,                       (2 66e)

                                         B +
where Bv~ and Bv+ are that part of variable boundary where the fluxes are directed into the region and out from


the region, respectively. FRATE(6) stores the fluxes through unspecified boundaries as


                                                  5

                           FRATE(6)=FRATE(7)-^FRATE(I)                    (2.67)
                                                 1=1

FRATE(8) and FRATE(9) store the accumulation rates in dissolved and adsorbed phases, respectively, as


given by
                                                                                  (2.68)
                                             R  5t
                                 FRATE(9)=f-^dR,                           (2.69)
                                             J   at
   at
R
                                           46

-------
FRATE(IO) stores the loss rate due to decay and FRATE(1 1) through FRATE(13) are set to zero as given by


                              FRATE(10)=pl(eC+pbS)dR ,                        (2.70)
                                          R

                        FRATE(11)=FRATE(12)=FRATE(13)=0  ,                  (2.71)

FRATE(14) is used to store the source/sink rate as

                  FRATE(14)=J[QCffi(l +sign(Q))+QC(l -sign(Q))]/2dR ,            (2.72)
                               R
If there is no numerical error in the computation, the following equation should be satisfied:

                                     14
                                                  O                              (2.73)
                                    1=7
and FRATE(6) should be equal to zero.  The calculation of the contributions to FRATE(8), FRATE(9),

FRATE(1 1), and FRATE(14) are in subroutine Q3R or Q4R.

Subroutine 04R or 03R

       These subroutines are used to compute the contributions to FRATE(8), FRATE(9), FRATE(1 1), and

FRATE(14) by quadrilateral and triangular elements, respectively:


                                    QRM=j6CdR ,                             (2.74a)
                                           R


                                     QDM=JSdR ,                              (2.74b)
                                           R


                    SOSM=J[QCm(l +sign(Q))+QC(l -sign(Q))]/2dR ,              (2.74c)
                            R

The computation of the above integrals is straightforward.
                                           47

-------
Subroutine ADVBC




       This  subroutine is called by HMCTRN to implement the boundary conditions.  For Dirichlet




boundaries, the Lagrangian concentration is specified. For variable boundaries, if the flow is directed out of




the region, the fictitious particle associated with the boundary node must come from the interior nodes. Hence




the Lagrangian concentration for the boundary node has already been computed from subroutine GNTRAK




and the implementation for such a boundary segment is bypassed. Thus, care should be taken to ensure that




the subroutine GNTRAK is called before the subroutine ADVBC.  If the flow is directed into the region, the




concentration of  incoming  fluid is  specified  for variable boundaries.   An  intermediate Lagrangian




concentration is calculated according to
                                                                                       (2?5a)
where C" is the concentration due to the boundary source at the boundary node i, Vn is the normal vertically




integrated Darcy's velocity, and Cm is the concentration of incoming fluid.




       Cauchy boundary conditions are normally applied to the boundary where flow is directed into the




region, and the material flux of incoming fluid is specified.  The intermediate Lagrangian concentration is thus




calculated according to
                                   Ci"=|NiqcdB/|NiVndB ,                           (2.75b)
where Q* is the concentration due to boundary source at the boundary node i, Vn is the normal component of




the Darcy's velocity, and qc is the Cauchy flux of the incoming fluid.




       The final Lagrangian concentration on either the variable or Cauchy boundary is obtained by using




the value C" and C", which is the concentration at the previous time step, as follows.
                                              48

-------
                                                                                         (2.75c)
Subroutine GNTRAK




        This subroutine is called by HMCTRN to control the process of backward or forward particle tracking




starting from global nodes.  It is designed to get the Lagrangian concentrations of all the fictitious particles




arriving at the global nodes at the current time step. In the subroutine, each particle is tracked one element by




one element until either the tracking time  is completely consumed or the particle encounters a specified




boundary side. During the particle tracking, this subroutine calls ELETRK to track a particle in the element




being considered. In order to make the particle tracking complete and remedy the given velocity field error




on the unspecified boundaries, subroutine FIXCHK calls subroutine ALGBDY to continue tracking particles




along the unspecified/Neumann boundaries.  At the end of each particle tracking, this subroutine calls INTERP




to calculate base functions such that the Lagrangian concentration can be computed by interpolation with those




base function values.




Subroutine HPTRAK




        This subroutine is called by HMCTRN to compute the locations and concentrations of all forward-




tracked fine-mesh nodes.  Basically, the algorithm of this subroutine is the same as that of subroutine




GNTRAK.




Subroutine ADVRX




        This subroutine is used to the following seven nonlinear simultaneous ordinary differential equations
                                               49

-------
       ^dsl
          Dt
                          21
                               C
        C
                                      K
                                        (3)
                                                  C
                C.
                                                 po
                                                 C
                                                       fc.
                                                                      (2.76a)
                           DC
                SO   ^s
                          c
                             c
                                                   O   ~0
                                                                      (2.76b)
-[(e+p*
                                         (3) + r
                                         po +Lp.
       y(2V(2)
       in Hn
                sn    s
Kf+C[
                                                    c
                                                   n   ^n
                                                                      (2.76c)
-[(0+p/,^
                                         P"   ^P
                                50

-------
           /         ^DCn
           (0+pbKdPhr =
uK,
(D ^
° Y0(1)
cs
K»)+C«.
co
Ko1) + C0
f Cp
po p
u(2)
(2) ^n
" V(2)
1n
cs
Ks(2)^
f C" 1
-rr \^-) f~\
11 + n.
f CP
pn + p
                                        (3)
                                    (3)
                                      .(3)
                                + e:
                                     Y
      K
                                                      c
                                  uc
                                                                                    (2.76d)
                      DC
                                           .(!)
                                C
                                                                                    (2.76e)
cs
v- (2) + r
Ksn +CS
f C" 1
K (2) r
n n.
cp
K (2) r
Nn +Lp.
                                                                            i
                                                                            A
                                                                                    (2.76f)
                      DC
                         (3)
                                           .(3)
                                              KS(03)+CS
                    K;
                                C
                              C,
    C,
                                                pn
                                                                                    (2.76g)
This subroutine is called right after the Lagrangian concentrations have been obtained.
                                             51

-------
Subroutine RXRATE




       This subroutine is called by subroutine ADVRX, and TASEMB during steady state simulations.




Basically the subroutine calculates the removal rates of 7 chemical components which are represented as the




terms within the braces on the right hand side of Eqs. (2.76a) through (2.76g).  The values of each bracket




within the braces are returned to the calling subroutines for each component.




Subroutine SFDET




       This subroutine determines if an element is a rough element based on the prescribed error tolerance




criteria shown in Eq.(2.46). If the M-th element is a rough element, the array IE(M,7) is activated to 1.




Subroutine FGDET




       This subroutine generates regular fine grids within each rough element based on the information of




IE(M,7) resulting from  subroutine SFDET.   It calls  subroutine HPTRAK to obtain the Lagrangian




concentrations of each activated fine grid.




Subroutine ISEHIL




       This subroutine removes all the forward-tracked nodes in a smooth element and stores the indices of




subelements into ISE array. In addition to regular refinement, subroutine ISEHIL also captures all the highest




and lowest concentrations within each subelement. The subelements, where the high-low points are located,




are determined by subroutine KGLOC. Once these high-low points are activated, subroutine TRIANG is




called to triangulate this subelement and the indices of each triangular are also stored in ISE array.




Subroutine TRIANG




       This subroutine calls subroutine CENTER and GRID to triangulate all the activated nodes in each




subelement according to the Delauney triangulation method.




Subroutine DFPREP




       This subroutine prepares all the needed information for assembling the fine grid elemental matrices.
                                              52

-------
It calls subroutines INTERPK to obtain the Lagrangian concentrations at each fine nodal points of the diffusion




step. This subroutine also prepares element indices for each subelement in the Eulerian step for composing




the matrix equation and stores the arrays for the intra-boundary points between rough and smooth regions to




circumvent the incompatibility by implementing the "slave point" concept.




Subroutine WARMSG




       The arguments passed to this subroutine are N, MAXN, SUBNAM, VARNAM, and NO. The stop




statement is activated whenever N is greater than MAXN, and a message is written in the output file to indicate




which variable is overflow in subroutine SUBNAM.




Subroutine INTERPK




       This subroutine has the same capability as subroutine INTERP.  Moreover, this subroutine once called




would interpolate the concentrations of all components instead of one component as in subroutine INTERP.




Subroutine KGLOC




       This subroutine is called by subroutine INTERP, INTERPK, and ISEHIL to obtain the subelement




on which the peak/valley point falls if the global element is a rough element.




Subroutine ELENOD




       This subroutine determines the number of nodes of element M by using the IE(M,4) information.




Subroutine VALBDL




       This subroutine calculates two interpolated values by the passed working arrays and basis functions.




Subroutine CENTER




       This subroutine is called by TRIANG to calculate the geometric center of each triangle.




Subroutine GRID




       This subroutine is called by subroutine TRIANG to triangulate all the nodes in each subelement




according to the Delaunay triangulation method. This  triangular information is used by GNTRAK for




backward-tracked interpolated concentration of next time step.
                                              53

-------
Subroutine NEWTRI




       This subroutine is called by subroutine GRID to produce new triangles in each subelement generated




in rough global elements.




Subroutine CKTRI




       This subroutine is called by subroutine GRID to check the existence of the triangles produced by




subroutine NEWTRI; i.e., to keep the triangles or not.




Subroutine ONLINE




       This subroutine adjusts the particle coordinates to be on the same line with the other two points.




Subroutine REPLAS




       This subroutine replaces the last four arguments with the first four arguments.




Subroutine WRKARY




       This subroutine prepares four real and one integer working arrays for later usage.




Subroutine REP3RI




       This subroutine replaces the last three real and one integer arguments with the first four arguments.




Subroutine MOVCHK




       This subroutine determines the concentrations and travel time of a fixed particle.




Subroutine ELETRK




       This subroutine is called by GNTRAK and HPTRAK to implement particle tracking in one element.




In order to increase the accuracy of particle tracking, regular refinement is considered to make the element be




composed of as many subelements as wished, and the tracking is executed one subelement by one subelement.




This subroutine is designed to control the process of particle tracking among the subelements. During the




particle tracking, this subroutine calls (1) TRACK 1 to track a particle in the subelement being considered if




that particle is right standing on a node of the subelement,  and (2) TRACK2 to track a particle if that particle
                                              54

-------
is not on any nodes of the subelement. The tracking will not be stopped until either the tracking time is




completely consumed or the particle encounters a side of the element.




Subroutine FIXCHK




        This is a control panel to check the ongoing process when a particle hits the boundary of the region




of interest. The backward tracked concentrations are obtained by interpolation if the boundary is  specified




including  Dirichlet, Cauchy, and Variable types.   Otherwise, the particle tracking continues along the




unspecified boundary till either the specified boundaries are encountered or tracking time is consumed.




Subroutine ALGBDY




        This subroutine is called by FIXCHK to control the process of backward or forward particle tracking




along the unspecified boundaries.  In the subroutine, the particle tracking is executed one boundary side by




one boundary side based on the nodal velocity component along the side being considered. The tracking will




not be  stopped until either the tracking time is completely consumed or the particle encounters a specified




boundary side. For accuracy, using the average velocity of both the source point and the target point to locate




the target point is firstly considered in the subroutine. However, if this average velocity approach is not able




to deal with very complex velocity fields, the single velocity of the source point is used to determine the




location of the target point.




Subroutine CKCNEL




        This subroutine checks the elements connected to a specific side segment.




Subroutine CKCOIN




        This subroutine checks to see if a specific point coincides with a global node.




Subroutine CKSIDE




        This subroutine checks to see if a specific point is on a side segment.




Subroutine INTERP




        This subroutine is called by GNTRAK and HPTRAK to obtain the Lagrangian concentrations by finite
                                               55

-------
element interpolation of the concentrations of the previous time step. A backward-tracked fictitious particle




can either land in a smooth element or a rough element.  If it lands in a rough element, a determination must




be made as to which subelement of the rough element it would land? This determination is made by calling




subroutine KGLOC.




Subroutine MMLOC




       This subroutine is called by ELETRK to locate the particle associated with a specific subelement for




subsequent elemental tracking. If this particle coincides with the nodes of a subelement, ICODE=0 is returned.




Subroutine TRACK 1




       This subroutine is called by ELETRK to compute the particle tracking in a specified subelement when




the source point coincides with a global node of the element. This subroutine determines (1) whether the




particle would move backwards or forwards into the element or not,  (2) which side of the element the particle




would head onto if the particle does move backwards or forwards into this element, and (3) where the target




point would be. After determining onto which side the particle is going to  move, this subroutine computes the




exact location of the target point on the side analytically.  For accuracy, using the average velocity of both the




source point and the target point to locate the target point is considered first in the subroutine.  However, if




this average velocity approach is not able to deal with very complex velocity fields, the single velocity of the




source point is used to determine the location of the target point.




 Subroutine TRACK2




       This subroutine is called by ELETRK to compute the particle tracking in a specified subelement when




the source point does not coincide with a global node of the element. This subroutine determines (1) whether




the particle would move backwards or forwards into the element,  (2) which side of the element the particle




would head onto if the particle does move backwards or forwards into this element, and (3) where the target




point would be. After determining which side the particle is going to move onto, this subroutine computes the




exact location of the target point on the side analytically.  For accuracy, using the average velocity of both the
                                               56

-------
source point and the target point to locate the target point is considered first in the subroutine.  However, if




this average velocity approach is not able to deal with very complex velocity fields, the single velocity of the




source point is used to determine the location of the target point.




Subroutine LOCQ2D




        This subroutine locates the target point of a particle tracking on a line segment in a specified element.




All the computations are made according to the average velocity approach and the single velocity approach,




as the index parameter IJUDGE is 1 and 2, respectively.  When the average velocity approach is considered,




the Newton-Ralphson method is used to solve a nonlinear algebraic equation such that the local coordinates




of the target point on the pre-determined element segment can be determined. With these local coordinates,




the location of the target point can be easily determined based on both the velocity of the source point and the




geometrical relationship between the source point and the pre-determined element side.




Function FCOS




        This function computes the cross of the vector of a given line segment with a specified vector whose




starting point stands on the line. The result helps to determine where  the endpoint of the specified vector is




located.




Subroutine IBWCHK




        This subroutine is called whenever the "in-element tracking" hits the global element side to determine




on which element side the target point is located.
                                                57

-------
       3.  ADAPTATION OF 2DFATMIC TO SITE SPECIFIC APPLICATIONS



       The following describes what one should do to the 2DFATMIC code for each site-specific

application and what data files should be prepared.


3.1    Control-integer Specifications


       For each site-specific problem, the users only need to specify the size of the problem by assigning

68 maximum control-integers with the PARAMETER statement in the MAIN program. The list and

definitions of the maximum control-integers required for both flow and transport simulations are given

below:

Maximum Control-Integers for the Spatial Domain

   MAXNPK = maximum no. of nodal points;
   MAXELK = maximum no. of elements;
   MXBESK = maximum no. of boundary-element surfaces;
   MXBNPK = maximum no. of boundary nodal points;
   MAXBWK = the maximum number of global nodes and/or fine nodes connected to any nodal point
              if the matrix equation is solved by PISS or maximum no. of bandwidth which is defined
              as 2*(maximum difference of global nodal number in an element)+l if the matrix
              equation is solved by the direct Gaussian elimination;
   MXJBDK = maximum no. of nodal points connected to any nodal point resulting from all global and
              all subelements surrounding that point;
   MXKBDK  = maximum no. of elements connected to any nodal point;
   MXADNK  = maximum no. of points used to solve matrix equations for transport;

Maximum Control-Integers for the Time Domain

   MXNTIK = maximum no. of time steps;
   MXNDTK = maximum no. of DELT changes;

Maximum Control-Integers for Material and Soil Properties

   MXMATK = maximum no. of material types;
   MXSPPK = maximum no. of soil parameters per material to describe soil characteristic curves;
   MXMPPK = maximum no. of material properties per material;
                                           58

-------
   The maximum control-integers for flow simulations and their definitions are given as the following:

Maximum Control-Integers for Source/sinks, flow

   MXSELF = maximum no. of source elements;
   MXSPRF = maximum no. of source profiles;
   MXSDPF = maximum no. of data points on each element source/sink profile;
   MXWNPF = maximum no. of well nodal points;
   MXWPRF = maximum no. of well source/sink profiles;
   MXWDPF = maximum no. of data points on each well source/sink profile;

Maximum Control-Integers for Cauchy Boundary Conditions, flow

   MXCNPF = maximum no. of Cauchy nodal points;
   MXCESF = maximum no. of Cauchy element surfaces;
   MXCPRF = maximum no. of Cauchy-flux profiles;
   MXCDPF = maximum no. of data points on each Cauchy-flux profile;

Maximum Control-Integers for Neumann Boundary Conditions, flow

   MXNNPF = maximum no. of Neumann nodal points;
   MXNESF = maximum no. of Neumann element surfaces;
   MXNPRF = maximum no. of Neumann-flux profiles;
   MXNDPF = maximum no. of data points on each Neumann-flux profile;

Maximum Control-Integers for Rainfall-Seepage Boundary Conditions, flow

   MXRNPF = maximum no. of variable nodal points;
   MXRESF = maximum no. of variable element surfaces;
   MXRPRF = maximum no. of rainfall profiles;
   MXRDPF = maximum no. of data point on each rainfall profile;

Maximum Control-Integers for Dirichlet Boundary Conditions, flow

   MXDNPF = maximum no. of Dirichlet nodal points;
   MXDPRF = maximum no. of Dirichlet total head profiles;
   MXDDPF = maximum no. of data points on each Dirichlet profile;
   The maximum control-integers for transport simulations and their definitions are given as:

Maximum Control-Integers for Source/sinks, transport

   MXSELT = maximum no. of source elements;
   MXSPRT = maximum no. of source profiles;
   MXSDPT = maximum no. of data points on each element source/sink profile;
   MXWNPT = maximum no. of well nodal points;
   MXWPRT = maximum no. of well source/sink profiles;

                                          59

-------
   MXWDPT = maximum no. of data points on each well source/sink profile;

Maximum Control-Integers for Cauchy Boundary Conditions, transport

   MXCNPT = maximum no. of Cauchy nodal points;
   MXCEST = maximum no. of Cauchy element surfaces;
   MXCPRT = maximum no. of Cauchy-flux profiles;
   MXCDPT = maximum no. of data points on each Cauchy-flux profile;

Maximum Control-Integers for Neumann Boundary Conditions, transport

   MXNNPT = maximum no. of Neumann nodal points;
   MXNEST = maximum no. of Neumann element surfaces;
   MXNPRT = maximum no. of Neumann-flux profiles;
   MXNDPT = maximum no. of data points on each Neumann-flux profile;

Maximum Control-Integers for Flowin-Flowout Boundary Conditions, transport

   MXVNPT = maximum no. of variable nodal points;
   MXVEST = maximum no. of variable element surfaces;
   MXVPRT = maximum no. of rainfall profiles;
   MXVDPT = maximum no. of data point on each rainfall profile;

Maximum Control-Integers for Dirichlet Boundary Conditions, transport

   MXDNPT = maximum no. of Dirichlet nodal points;
   MXDPRT = maximum no. of Dirichlet total head profiles;
   MXDDPT = maximum no. of data points on each Dirichlet profile;

Control-Integers for Number of Components in the system

   MXNCCK = maximum no. of components in this system;

Maximum Control-Integers for Refined System

   MXKGLDK = maximum no. of subelements in the Eulerian step;
   MXLSVK = maximum no. of subelement sides located on the intra-boundaries between extended
              rough and smooth regions;
   MXMSVK = maximum no. of global element sides located on the intra-boundaries between extended
              rough and smooth regions;
   MXNDBK = maximum no. of diffusion fine nodal-points located on the global boundary;
   MXNEPK = maximum no. of all forward tracked nodal points in the region of interest when the exact
              peak capture and oscillation free (EPCOF) numerical scheme is used. When EPCOF is
              not used, set MXNEPK = 1;
   MXEPWK = maximum no. of forward tracked nodal points in any rough element when the exact
              peak capture and oscillation free (EPCOF) numerical scheme is used. When EPCOF is
              not used, set MXEPWK = 1;
                                          60

-------
    MXNPWK = maximum no. of fine nodal-points in any global element for particle tracking;
    MXELWK = maximum no. of subelement elements in any global element for particle tracking;
    MXNPWS = maximum no.  of fine nodal-points in  any global element which surrounds point
              sources/sinks  for   obtaining  more  accurate  Lagrangian  concentrations  with
              injection/extraction wells in the region of interest;
    MXELWS = maximum no. of subelements in any global element which surrounds any point wih a
              source/sink  for   obtaining  more   accurate   Lagrangian  concentrations  with
              injection/extraction wells in the region of interest.
    MXNPFGK = maximum no. of forward tracked nodal points over the region of interest or maximum
              no. of fine nodal points plus peak/valley nodal points;
    MXKGLK = maximum no. of subelements in the  Lagrangian step;
    For flow simulations only, we demonstrate how to specify the above maximum control-integers with

PARAMETER statement in the MAIN in the following by an example.



    Assume that a region of interest is discretized by 11x11 nodes and 10x10 elements (Fig. 3.1).  In

other words, the region is discretized with 11  nodes along the longitudinal or x-direction and 11 nodes

along the vertical or z-direction.  Since the region has a total of nodes  11x11 = 121, the maximum

number of nodes is MAXNPK =121. The total number of elements is IQx 10 = 100; i.e. MAXELK = 100.

For this simple discretization problem, the maximum number of nodes connecting to any of the 121 nodes

in the region of interest is 9; i.e., MXJBDK = 9; and the maximum number of elements connecting to any

of the 121 nodes is 4; i.e. MXKBDK = 4. There will be 10 element segments each on the bottom and top

faces of the region and 10 element segments each on the left and right faces of the region.  Thus, there

will be a total of 40 element segments; i.e., MXBESK = 40 .  Similarly, we can compute the boundary

nodes to be 40, i.e., MXBNPK = 40. Assume that the direct Gaussian elimination method is used to solve

the matrix equation. Then MAXBWK is set equal to the half bandwidth plus 1.  For this discretization,

the  computed half band width plus 1 is: 2*(13-1)+1 = 25; thus, MAXBWK = 25.  Since this is a flow

problem, no local grid refinement is involved; thus MXADNK = MAXNPK+0.
                                            61

-------
        11

        10


        9


        8

        7
        6
  Cauchy
  Boundary
        5
©
0
©
           0
           ©
           ©
           ©
                                                Variable Boundary
                12      23       34      45      56      67
                                       Neumann Boundary
                                                                     89
                                                                                       21
                                                                           >120
\116
iDirichlet
I Boundary
'
                                                                                      m
                                                                                      112
                                                                   1001T1
   Figure 3.1 Finite Element Discretization and Source/Sink Layout of an Illustrated Example for Flow



    In this model, six material properties (six saturated hydraulic conductivity components) per material

were included for flow problems. Assume that the whole region of interest is made of three different

kinds of materials. The characteristic curves of each material are assumed  to be described by four

parameters. We then have MXMATK = 3. MXMPPK = 6. and MXSPPK = 4. If we assume that we will

make a 500 time step simulation and we  will reinitiate the time step size for 20 times during  our

simulation, then we have MXNTIK = 500 and MXNDTK = 20.
                                              62

-------
    We will assume that there will be a maximum of 11 elements that have the distributed sources/sinks




(i.e., MXSELF =11) and a maximum of 10 nodal points that can be considered as well sources/sinks (i.e.,




MXWNPF= 10). We will also assume that there will be two different distributed source/sink profiles and




four distinct point source/sink profiles.  Then we will have MXSPRF = 2 and MXWPRF = 4.  Let us




further assume that four data points are needed to  describe the  distributed source/sink profiles as  a




function of time and that 8 data points are required to describe point source/sink profiles (i.e., MXSDPF




=_! and MXWDPF = 8).









    To specify maximum control-integers for the boundary conditions, let us assume that the top face is




a variable boundary (i.e., on the air-soil interface, either ponding, infiltration, or evapotranspiration may




take place). On the left face, fluxes from the adjacent aquifer are known.  On the right face, the total head




is assumed known.  On the bottom face, natural drainage is assumed to occur (i.e., the gradient of the




pressure head can be assumed zero). There are 11 nodes on the left face and 10 element segments; thus




MXCNPF =11 and MXCESF = 10. It is further assumed that there are two different fluxes going into




the region through the left face and that each flux can be described by four data points as a function of




time (i.e., MXCPRF = 2, and MXCDPF = 4).  On the bottom surface, there are 11 nodes and 10 element




segments. Since the gradient of pressure head on the bottom surface is zero, there is only one Neumann




flux profile, and two data points, one at zero time and the other at infinite time, are sufficient to describe




the constant value of zero.   Hence, we have MXNNPF =  11. MXNESF = 10.  MXNPRF = 1. and




MXNDPF = 2. On the top face, there will be 11 nodes and 10 element segments. Let us assume that there




are three  different rainfall intensities that might be falling on the air-soil interface, and that each rainfall




intensity  is a function of time and can be described by 24 data points. With these descriptions, we have




MXRNPF= 11. MXRESF = 10.  MXRPRF = 3. and MXRDPF = 24.  On the right face, there are 11




nodes. Let us assume that there are eleven different values of the total head, one each on a vertical line
                                             63

-------
of the right face. We further assume that each of these twenty total heads can be described by 8 data

points as a function of time. We then have MXDNPF = 11. MXDPRF = 11. and MXDDPF = 8.



   From the above discussion, the following PARAMETER statements can be used to specify the

maximum control-integers in the MAIN for the problem at hand:
   PARAMETER(MAXELK=100,MAXNPK=121,MXBNPK=40,MXBESK=40,MAXBWK=25)
   PARAMETER(MXJBDK=9,MXKBDK=4,MXADNK=MAXNPK+0)
   PARAMETER(MXMATK=3,MXSPPK=4,MXMPPK=6)
   PARAMETER(MXNTIK=500,MXNDTK=20)

   PARAMETER(MXSELF=11,MXSPRF=2,MXSDPF=4,MXWNPF=10,MXWPRF=4,MXWDPF=8)
   PARAMETER(MXCNPF= 11 ,MXCESF=10,MXCPRF=2,MXCDPF=4)
   PARAMETER(MXNNPF=11 ,MXNESF=10,MXNPRF= 1 ,MXNDPF=2)
   PARAMETER(MXRNPF=11,MXRESF=10,MXRPRF=3,MXRDPF=24)
   PARAMETER(MXDNPF=11,MXDPRF=20,MXDDPF=8)

   PARAMETER(MXSELT=1,MXSPRT=1,MXSDPT=1,MXWNPT=1,MXWPRT=1,MXWDPT=1)
   PARAMETER(MXCNPT=1,MXCEST=1,MXCPRT=1,MXCDPT=1)
   PARAMETER(MXNNPT= 1 ,MXNEST= 1 ,MXNPRT= 1 ,MXNDPT= 1)
   PARAMETER(MXVNPT= 1 ,MXVEST= 1 ,MXVPRT= 1 ,MXVDPT= 1)
   PARAMETER(MXDNPT= 1 ,MXDPRT= 1 ,MXDDPT= 1)

   P ARAMETER(MXNCCK= 1)

   PARAMETER(MXKGLDK= 1 ,MXLS VK= 1 ,MXMSVK= 1 ,MXNDBK= 1)
   PARAMETER(MXNEPK= 1 ,MXEPWK= 1)
   P ARAMETER(MXNPWK= 1 ,MXELWK= 1 ,MXNPWS=1 ,MXELWS=1)
   P ARAMETER(MXNPFGK= 1 ,MXKGLK= 1)
   In the following, for transport simulations only, we demonstrate how to specify the maximum control-

integers with PARAMETER statements in the MAIN with an example.



   Let us assume that a region of interest is discretized by 11x11 nodes and 10x10 elements (Fig. 3.2).

In other words, we are discretizing the region with 11 nodes along the longitudinal or x-direction and 11

nodes along the vertical or z-direction.  Since we have a total of nodes 11x11 = 121, the maximum

                                      64

-------
number of nodes is MAXNPK =121. The total number of elements is IQx 10 = 100; i.e. MAXELK = 100.




For this simple discretization problem, the maximum number of nodes connecting to any of the 121 nodes




in the region of interest is 17; i.e., MXJBDK =17. (e.g., there are four global elements connected to




global node 25, thus there are 9 global nodes connected to  node 25.   In addition, there are  four




subelements surrounding node 25, thus 8 additional fine nodal points are connected to node 25.  Hence




the total number of global nodes and fine nodal points connected to node 25  is 17.)  The maximum




number of elements connecting to any of the 121 nodes is 4; i.e. MXKBDK = 4. There will be 10 element




segments each on the bottom and top faces of the region and 10 element segments each on the left and




right faces of the region.  Thus, there will be  a total of 40 element segments; i.e.,  MXBESK = 40 .




Similarly, the number of boundary nodes is 40; thus MXBNPK = 40 (Fig. 3.2).
                                            65

-------
           Dirichlet Boundary
                                              Cauchy Boundary
       101 P
       91 F
      8


      7


Variable 6
Boundary
      5


      4


      3
                 12      23       34      45      56      67      78
                                       Neumann Boundary
                                                                                      121
                                                                                      120
                                                                                      119
                                                                                      118
                                                                                      117
                                                                                      116
                                                                                      Variable
                                                                                      Boundary
                                                                                      115
                                                                                      114
                                                                                      113
                                                                                      112
                                                                           100
                                                                                    111
 Figure  3.2 Finite Element Discretization and Source/Sink Layout of an Illustrated Example for Transport

    We assume that there are 25 rough elements as shown in Figure 3.2. Each rough element is refined

by 4 (2 x 2 = 4) subelements in the diffusion step computation.  It is seen (Fig. 3.2) that the maximum

number of nodal points (including both global and fine nodal points) connecting to any nodal points is

9. We assume that adaptive zooming will be enforced in the diffusion step of solving the transport

equation, which necessitates the use of point iteration method  to solve the matrix equation.  Hence,

MAXBWK is set to be the maximum number of nodal points connecting to any node; i.e. MAXBWK =

iL This is different from the flow case described previously, where the direct Gaussian elimination

method was used to solve the matrix equation; thus MAXBWK was set to be the maximum half

bandwidth plus 1.
                                              66

-------
    Five additional diffusion points for each rough element will be added to the region of interest with




a 2 x 2 subelement refinement. When some of the rough elements are connected, some of these diffusion




points will be duplicated.  For example, there are 101 additional diffusion points, not 5 X 25 = 125 points




because some of the rough elements are connected.  Thus, MXADNK = MAXNPK +101. If all 25 rough




elements were disconnected,  MXADNK would be specified as MXADNK = MAXNPK  + 5 x 25 =




MAXNPK +125.  For sake of safety, sometime we may want to assign the latter value.




    In this model, we have included eight material properties per material for transport problems. We will




assume that the whole region of interest is composed of three different kinds of materials. We then have




MXMATK = 3 and MXMPPK = 8. If we assume that we will make a 500 time step simulation and we




will reinitiate the change on the time step size for 20 times during  our simulation, then we have MXNTIK




= 500 and MXNDTK = 20.




    Because there are 25 rough elements and each rough element  is refined by 2 x 2 = 4 subelements, the




maximum number of subelements for assembling  in the diffusion step is MXKGLDK = 25 x 4 = 100.




Figure 3.2 shows that the number of global element sides located  on the  intra-boundaries between rough




and smooth regions is 51,  and the number of subelement sides on these intra-boundaries is 102.  Hence,




MXMSVK = 51. and MXLSVK= 102. The number of diffusion fine nodal points on the global boundary




is MXNDBK =12 (Fig. 3.2).  The values of MXMSVK naturally depends on where and how the rough




elements are located in  the region of interest.  In the worst  case  when  all  25 rough elements are




disconnected, the MXMSVK would be  equal to 25 x 4 = 100 because each rough element would be




surrounded by four smooth elements.  The value of  MXLSVK would depend on how each rough element




is refined. For this particular case, MXLSVK =  100 x 2 = 200 because every global element side is




refined by 2 subelement sides.
                                            67

-------
   We will assume that there will be a maximum of 11 elements that have the distributed sources/sinks




(i.e., MXSELT= 11) and a maximum of 10 nodal points that can be considered as well sources/sinks (i.e.,




MXWNPT =10) as shown in Figure 3.2. We will also assume that there will be two different distributed




source/sink profiles and four distinct point source/sink profiles (Fig. 3.2). Then we will have MXSPRT




= 2 and MXWPRT = 4. Let us further assume that four data points are needed to describe the distributed




source/sink profiles as a function of time and that 8 data points are required to describe point source/sink




profiles (i.e., MXSDPT = 4 and MXWDPT=8).




   To specify maximum control-integers for boundary conditions, let us assume that nodes 9, 10, 11, 22,




and 33 (Fig. 3.2) are the Dirichlet boundary nodes.  The remaining top face is a Cauchy boundary. The




remaining left face and the whole right face are the variable boundaries (i.e., either Cauchy condition with




known incoming concentration or Neumann condition with zero gradient of concentration).  On the




bottom face, the Neumann condition is assumed to occur with zero concentration gradient.  For the five




Dirichlet nodes, let us assume that there are two profiles to be applied to these five nodes. We further




assume that each profile can be described by 8 data points for the Dirichlet boundary values. From this




description, we have MXDNPT = 5, MXDPRT = 2, and MXDDPT = 8.  There are 8 element segments




which consist of 9 nodes on the remainder of the top  face; thus MXCNPT =  9 and MXCEST = 8. It is




further assumed that there are two different fluxes going into the region through the top face and that each




flux can be described by four data points as a  function of time (i.e., MXCPRT = 2, and MXCDPT = 4).




On the bottom face, there are  10 element segments which consist of 11 nodes. Since the gradient of




concentration on the bottom surface is zero, there is only one Neumann flux profile, and two data points,




one at zero time and the other at infinite time, are sufficient to  describe the constant value  of  zero.




Hence, we have MXNNPT=11. MXNEST=10. MXNPRT=  1. and MXNDPT = 2. On part of the left




face and the whole right face, there are a total of 18 element segments, which consist of 20 nodal points.




Let us assume that there are two concentration profiles imposed on the left face  and one concentration
                                            68

-------
imposed on the right face, that is a total of 3 concentration profiles are used for the variable boundary




condition.  We further assume that each concentration profile is a function of time and can be described




by 24 data points. With these descriptions, we have MXVNPT = 20. MXVEST= 18. MXVPRT = 3. and




MXVDPT = 24.




    There are seven components, say microbe #1, microbe #2, microbe #3, substrate, oxygen, nitrate, and




nutrient, involved in this system; i.e., MXNCCK = 7.




    The numerical schemes for solving transport equations are LEZOOMPC plus keeping EPCOF points




in the Lagrangian step.  For practical problems, EPCOF points will not  be kept; thus, MXNEPK = 1.




MXEPWK = 1. In the Lagrangian step, each element is assumed to be refined by 4 subelements (NXW




= 2 and NZW = 2) for accurate tracking. With this assumption, we have MXNPWK = (NXW+1)  x




(NZW+1) = 9. MXELWK = (NXW x NZW) = 4.  For each element connected to the sources/sinks, we




assume to refine it with 3x3 elements for accurate computation of Lagrangian concentrations to yield




MXNPWS = 16 and MXELWS = 9.




    The specification of MXNPFGK and MXKGLK is much more involved. These two control integers




depend on many things: (1) how all the nodal points (including global  nodes and fine nodal points) at the




beginning of a time-step simulation are forwardly tracked, (2) how many elements are rough at the end




of the time-step computation,  (3) how each rough element is refined, and (4) how many peak/valley




points are kept.  Figure 3.3 helps to illustrate how  these two control integers are determined.  Let us




assume that at the beginning of a time-step computation, there are 24 global elements and 35 global nodes




(Fig. 3.a).  We further assume that four elements are rough and each rough element is refined by 2 x 2  =




4 subelements. The total number of nodes in the region of interest is  51: nodes 1 through 35 are global




nodes and nodes 36  through 51 are fine nodal points. At this moment,  the number of subelements would




be 4x4 = 16, and the number of fine nodal points would be 4 x (3 x 3) = 36 (in computing the number
                                            69

-------
of fine nodal points we have included the global nodes in rough elements) even though there are only 26




fine nodal points.
           29
30
31
32
33
34
35
22
4$
15
37}
Q
c
Y.. 1
V
( >
s
}
k\
\
W
( \
4tf '
f
4CP
f s
3> )
f
'36°
V
23 }
( \
47° ;
•s
16 )
^ s
39° ;
•s
51°
< \
49^ '
t
44"
^s
;
/
9 74f

24
5CP
17
43°
10

25
18
11

26
19
12

27
20
13
23456
28
21
14
7
               X



       Figure  3.3a An Illustrated Example of Global Nodes and Fine Nodes for Initial Condition
    According to the LEZOOMPC algorithm, we need to perform forward node tracking.  Since we have




seven components, we need to perform a forward node tracking with respect to each component.  Let us




assume, the forward tracking velocity with respect to Component 1  is given in Figure 3.3(b). We need




to forward track all 51 nodal points. After forward tracking, some of these 51 points will remain in the
                                             70

-------
region of interest and some of them will go out of the region.  For this particular case, 32 nodal points
remain in the region of interest (Fig. 3.3b).  These 32 nodal points are referred to as the intermediate-fme-
nodal points with respect to Component 1.  Similarly, we have to perform forward tracking of 51 nodal
points with respect to Components 2 through 7.
 Y*
           29
    Flow direction:
30          31          32
33
34
          22
           15
                         23
                         16
                                            ^2*
             24
                                            x
                                            ^5*


 26
                                                11
 27
                                    19
                                                              x
                                                              ^42*
                                    12
             20
             13
35
                                                            28
                                                            21
                                                                                    14
              X
              Figure 3.3b  An Illustrated Example of Forward-tracked Fine Nodal Points

    At the end of forward tracking, the number of fine nodal points would be the sum of all intermediate-
fine-nodal points that still remain in the region of interest.  In general, each chemical component has a
tracking velocity of its own; thus there will be a set of intermediate-fine-nodal points corresponding to
                                            71

-------
each component.  Let us assume that the numbers of intermediate-fine-nodal points corresponding to




Components 2 through 7 are 30, 34, 35, 36, 33, and 40,  respectively. The number of fine nodal points




would then be 32 + 30 + 34 + 35 + 36 + 33 + 40 = 240.




    The next step in LEZOOMPC algorithm is the determination of rough elements.  Let's assume that




we  have determined 4 elements to be rough based on all the forward tracked nodes, and these four




elements happen to be the elements formed by global nodes 11-12-19-18, 12-13-20-19, 18-19-26-25, and




19-20-27-26, respectively (Fig. 3.3c).  The final step in the LEZOOMPC is the regular refinement and




peak/valley capturing at the end of the time-step simulation.  These four rough elements are refined with




4 x  (2 x 2) = 16 subelements and 4 x (3 x 3) = 36 fine nodal points.  Let us further assume that there is




a peak/valley nodal point in each of 9 subelements out of the 16 subelements (Fig. 3.3c).  The number




of fine  nodal points at the end of the time-step computation would be equal to 36  + 9 = 45 (Fig. 3.3c).




For those 9 subelements that have  peak/valley nodal points, each has to be triangularized. The number




of triangles  in each subelement depends on the number of peak/valley nodal points. The  number of




triangles can be computed with 2Np + 2 (where Np is the number of peak/valley nodal points).  For this




particular case, each subelement with a peak/valley nodal point can  be triangularized with 4 triangular




elements.  Thus, the  number of subelements after regular refinement and triangulation is  equal to  9




(subelements with one peak/valley nodal points) x  4 + 7 (subelements without peak/valley nodal points)




= 43 (Fig. 3.3c).
                                            72

-------
           29
30
31
32
33
34
35
              X
           Figure 3.3c An Illustrated Example of Fine Nodal Points and Peak/Valley Points
                           at the End of the Time Step for Advection
                                                                                 14
   The maximum number of fine-nodal points for the present time-step would be MXNPFGK = max

(36,240, 45) = 240. The maximum number of subelements would be MXKGLK = max (16,43) = 43. The

specification of MXNPFGK and MXKGLK should be large enough for all time steps, not just for a

particular time step.



   From the above discussion, the following PARAMETER statements can be used to specify the

maximum control-integers in the MAIN for the problem at hand:
                                           73

-------
PARAMETER(MAXNPK=121,MAXELK=100,MXBNPK=40,MXBESK=40,MAXBWK=9)
PARAMETER(MXJBDK=9,MXKBDK=4,MXADNK=MAXNPK+125)
PARAMETER(MXMATK=3,MXSPPK=6,MXMPPK=8)
PARAMETER(MXNTIK=500,MXNDTK=20)

PARAMETER(MXSELF= 1 ,MXSPRF=1 ,MXSDPF= 1 ,MXWNPF=1 ,MXWPRF=1 ,MXWDPF= 1)
PARAMETER(MXCNPF=1 ,MXCESF= 1 ,MXCPRF=1 ,MXCDPF= 1)
PARAMETER(MXNNPF= 1 ,MXNESF= 1 ,MXNPRF= 1 ,MXNDPF= 1)
PARAMETER(MXRNPF=1 ,MXRESF= 1 ,MXRPRF=1 ,MXRDPF= 1)
PARAMETER(MXDNPF= 1 ,MXDPRF=1 ,MXDDPF= 1)

PARAMETER(MXSELT=11,MXSPRT=2,MXSDPT=4,MXWNPT=10,MXWPRT=4,MXWDPT=8)
PARAMETER(MXCNPT=9,MXCEST=8,MXCPRT=2,MXCDPT=4)
PARAMETER(MXNNPT= 11 ,MXNEST= 10,MXNPRT= 1 ,MXNDPT=2)
PARAMETER(MXVNPT=20,MXVEST= 18,MXVPRT=3 ,MXVDPT=24)
PARAMETER(MXDNPT=5,MXDPRT=2,MXDDPT=8)

PARAMETER(MXNCCK=7)

PARAMETER(MXKGLDK= 100,MXLS VK= 102,MXMSVK=51 ,MXNDBK= 12)
PARAMETER(MXNEPK= 1 ,MXEPWK= 1)
PARAMETER(MXNPWK=9,MXELWK=4,MXNPWS=16,MXELWS=9)
PARAMETER(MXNPFGK=240,MXKGLK=43)
3.2 Soil Property Function Specifications


   Analytical functions are used to describe the relationships of water content, water capacity, and

relative hydraulic conductivity with the pressure head.  Therefore, the user must supply three functions

to compute the water content, water capacity, and relative hydraulic conductivity based on the current

value of the pressure head. The parameters needed to specify the functional form are read and stored in

SPP. One example is shown in the subroutine SPROP in the source code. In this example, the water

content, water capacity, and relative hydraulic conductivity are given by (van Genuchten 1980):

                                          e -e
                               6  =  e
                                    r
                                    1
                                       74

-------
                                = a(n-i)[i-f(e)]m[f(6)](es-er)                      (3.2)
in which
                                 f(6) = [6-er]/[es-er]1/m                           (3.4)

and

                                       m =1  - -                                 (3.5)
    To further demonstrate how we should modify the subroutine SPROP in Appendix A to accommodate

the material property functions that are different from those given by Eqs. (3.4) through (3.5), let us

assume that the following Fermi types of functions are used to represent the unsaturated hydraulic

properties (Yeh 1987):

                           e = er  + (es-er)/{i+exP[-a(h-he)]}                     (3.6)

                   d6/dh = a(es-er)exp[-oc(h-he)]/{l+exp[-a(h-he)]}2 ,             (3.7)

and
                          log10(Kr) = €/{l+exp[-p(h-hk)]} - e ,                    (3.8)


where 6S, 6r, a,  and he are the parameters for computing the water content and water capacity; and P, e,

and hk are the parameters for computing the relative hydraulic conductivity. The source code, subroutine

SPFUNC, must be changed, for this example, to the following form for computing the moisture content

and water capacity
     WCR=SPP(1,MTYP,1)
     WCS=SPP(2,MTYP,1)
     ALPHA=SPP(3,MTYP, 1)
     HTHETA=SPP(4,MTYP, 1)
     EPS=SPP(1,MTYP,2)


                                           75

-------
     BETA=SPP(2,MTYP,2)
     HSUBK=SPP(3,MTYP,2)
C
C	SATURATED CONDITION
C
     IF(HNP.LE.O.O) THEN
      TH=WCS
      IF(ISP.EQ.l) GOTO 900
      DTH=O.ODO

      USKFCT=1.0DO
     ELSE
C
C	UNSATURATED CASE
C
     EXPAH=DEXP(-ALPHA*(HNP-HTHETA))
     TH=WCR+(WCS-WCR)/(1 .ODO+EXPAH)
     IF(ISP.EQ.l) GOTO 900
     DTH=ALPHA*(WCS-WCR)*EXPAH/(1.0DO+EXPAH)**2
     AKRLOG=EPS/(1.0DO+DEXP(-BETA*(HNP-HSUBK))) - EPS
     UKFCT=10.0DO**AKRLOG
3.3 Input and Output Devices


   Five logical units are needed to execute 2DFATMIC. Units 15 and 16 are standard card input and

line printer devices, respectively.  Unit 11 must be specified to store the flow simulation results, which

can be used for plotting purposes. Unit 12 must be specified to store the transport simulation results,

which can be used for plotting purposes. Unit 13 is used to store the boundary arrays for later uses, if

these arrays are computed for the present job. Unit 14 is used to store pointer arrays for later uses, if these

arrays are generated for the present job. For large problems, our experience has indicated that it would

take too much time to process the boundary arrays and to generate pointer arrays. Hence, it is advisable

that for multijob executions, these boundary and pointer arrays should be computed only once and written

on units 13  and 14, respectively. Once they are stored on units 13 and 14, the IGEOM described in

Appendix A should be  properly  identified for the new job so they can be read via units 13 and  14,

respectively.
                                          76

-------
                               4. SAMPLE PROBLEMS






       To verify 2DFATMIC, three illustrative examples are used. The first one is a one-dimensional


flow problem through a column. The second one is a one-dimensional transport problem through a


column. The third one is a two-dimensional biodegradation problem.
4.1 Problem No. 1: One-Dimensional Column Flow Problem



Problem Description


       This example is selected to represent the simulation of a one-dimensional problem with the


revised FEMWATER. The column is 200 cm long and 50 cm wide (Figure 4.1). The column is assumed


to contain the soil with a saturated hydraulic conductivity of 10 cm/day, a porosity of 0.45 and a field


capacity of 0.1.  The unsaturated characteristic hydraulic properties of the soil in the column are given


as:



                                                / h -  h '
                               e  = es - (es - ej  ——i-                          (4.1)



and

                                            e  - er

                                      Kr = y-rf                                (4.2)
                                             s     r



where hb and ha are the parameters used to  compute the water  content  and the relative hydraulic


conductivity.
                                           77

-------
                      Variable B.C.
zero flux
                        h = 0.0
                                         200cm
                                                  zero flux
                        50 cm





Figure 4.1  Problem definition for the one-dimensional transient flow in a soil column



                                  78

-------
       The initial conditions are assumed: a pressure head of-90 cm is imposed on the top surface of the

column, 0 cm on the bottom surface of the column, and -97 cm elsewhere.  The boundary conditions are

given as: no flux is imposed on the left, front, right, and back surfaces of the column; pressure head is

held at 0 cm on the bottom surface; and variable condition is used on the top surface of the column with

a ponding depth of zero, minimum pressure of -90 cm, and a rainfall of 5 cm/day for the first ten days and

a potential evaporation of 5 cm/day for the second ten days.

       The region of interest, i.e. the whole column, will be discretized with 1 x 40 = 40 elements with

element size = 50 x 5 cm, resulting in 2 x 41  = 82 node points.  A variable time-step size is used. The

initial time step size is 0.05 days and each subsequent time-step size is increased by 0.2 times with a

maximum time-step size not greater than 1.0 day. Because there is an abrupt change in the flux value

from 5 cm/day (infiltration) to -5 cm/day (evaporation) imposed on the top surface at day 10, the time-step

size is automatically reset to 0.05 day on tenth day.  A 20 day simulation will be made with the revised

FEMWATER.  With the time-step size described above, 44 time-steps will be needed. The pressure head

tolerance for nonlinear iteration is 2 10"2 cm, and the relaxation factors is set equal to 0.5.

       To execute the problem, the maximum control-integers in the MAIN program should be specified
as:
   PARAMETER(MAXELK=40,MAXNPK=82,MXBESK=82,MXBNPK=82,MAXBWK=7)
   PARAMETER(MXJBDK=9,MXKBDK=4,MXADNK=MAXNPK+1)
   PARAMETER(MXMATK= 1 ,MXSPPK=4,MXMPPK= 10)
   PARAMETER(MXNTIK=44,MXNDTK=3)

   PARAMETER(MXSELF=1 ,MXSPRF=1 ,MXSDPF=1 ,MXWNPF=1 ,MXWPRF=1 ,MXWDPF= 1)
   PARAMETER(MXCNPF=1 ,MXCESF= 1 ,MXCPRF=1 ,MXCDPF= 1)
   PARAMETER(MXNNPF= 1 ,MXNESF= 1 ,MXNPRF= 1 ,MXNDPF=1)
   PARAMETER(MXRESF=1 ,MXRNPF=2,MXRPRF=1 ,MXRDPF=4)
   PARAMETER(MXDNPF=2,MXDPRF=1 ,MXDDPF=2)

   PARAMETER(MXSELT= 1 ,MXSPRT= 1 ,MXSDPT= 1 ,MXWNPT= 1 ,MXWPRT= 1 ,MXWDPT= 1)
   PARAMETER(MXCNPT= 1 ,MXCEST= 1 ,MXCPRT= 1 ,MXCDPT= 1)
   PARAMETER(MXNNPT= 1 ,MXNEST= 1 ,MXNPRT= 1 ,MXNDPT= 1)
   PARAMETER(MXVEST= 1 ,MXVNPT= 1 ,MXVPRT= 1 ,MXVDPT= 1)
   PARAMETER(MXDNPT= 1 ,MXDPRT= 1 ,MXDDPT= 1)

                                         79

-------
   PARAMETER(MXNCCK= 1)

   PARAMETER(MXKGLDK= 1 ,MXLSVK= 1 ,MXMS VK= 1 ,MXNDBK= 1)
   PARAMETER(MXNEPK= 1 ,MXEPWK= 1)
   PARAMETER(MXNPWK= 1 ,MXELWK= 1 ,MXNPWS= 1 ,MXELWS= 1)
   PARAMETER(MXNPFGK= 1 ,MXKGLK= 1)
      To reflect the soil property function given by Eqs. (4.1) and (4.2), we have to modify the

subroutine SPFUNC. A segment of the code in the subroutine SPFUNC must be modified as:

C %%%% EXAMPLE 1
C
   WCR=SPP(1,MTYP,1)
   WCS=SPP(2,MTYP,1)
   HAA=SPP(3,MTYP,1)
   HAB=SPP(4,MTYP,1)
CC
CC	SATURATED CONDITION
CC
C
C %%%% EXAMPLE 1
C
   IF(HNP.GT.O.O) GO TO 700
C
    TH=WCS
    IF(ISP.EQ.1)GOTO900
    DTH=O.ODO
   AKHCX=SATKX*RHOM/AMU
   AKHCZ=SATKZ*RHOM/AMU
   AKHCXZ=S ATKXZ * RHOM/AMU
   GOTO 900
CC
CC	UNSATURATED CASE
CC
  700   CONTINUE
C
C %%%% EXAMPLE 1
C

   THMKG=WCs-(WCS-WCR)*(-hnp-haa)/(hab-haa)
   TH=THMKG
   IF(ISP.EQ.l) GOTO 900
   rk=(thmkg-wcr)/(wcs-wcr)
   DTH=-(wcs-wcr)/(hab-haa)
CC
   IF(RK.LT.CNSTKR)RK=CNSTKR

                                     80

-------
    AKHCX=S ATKX* rk* RHOM/AMU
    AKHCZ=SATKZ*rk*RHOM/AMU
    AKHCXZ=S ATKXZ * rk* RHOM/AMU
   ENDIF
C
 900 RETURN
   END
Input and Output for Problem No. 1

       With the above descriptions, the input data can be prepared according to the instructions given

in Appendix A. The input parameters are shown in Table 4.1 and the input data are given in Table 4.2.

The simulation results are identical to those obtained with FEMWATER. Thus, the flow module of

2DFATMIC is verifiied.  To save space, output is available in electronic form.
                                           81

-------
Table 4.1  The list of input parameters for Problem 1
Parameters
number of points
AX
AZ
KS77
er
es
ha
hb
initial time-step size
time-step size increment
percentage
maximum time-step size
no. of times to reset
time-step size
time to reset time-step
size
total simulation time
no. of time-steps
tolerance for nonlinear
iteration
relaxation factor for
nonlinear iteration
Pw
Hw
g
Notation in the
data input guide
NNP
XAD
ZAD
PROPf(l,2)
SPP(1,1,1)
SPP(2,U)
SPP(3,1,1)
SPP(4,1,1)
DELT
CHNG
DELMAX
NDTCHG
TDTCH(l)
TMAX
NTI
TOLBh
OMEh
PROP(I,4)
PROP(I,9)
PROP(I,5)
Value
82
50
5
10
0.15
0.45
0
-100
0.05
0.2
1
1
10
22
44
2xlQ-2
0.5
1.0
9483.264
7.316xl012
Unit
dimensionless
cm
cm
cm/day
dimensionless
dimensionless
cm
cm
day
dimensionless
day
dimensionless
day
day
dimensionless
cm
dimensionless
g/cm3
g/cm/day
cm/day2
Data set
3. A.
3.B.
3.B.
6. B.
7. B.
7. B.
7. B.
7. B.
8.B.
8.B.
8.B.
8. A.
8.E.
8.B.
8. A.
9. B.
9. B.
6.B.
6.B.
6.B.
                       82

-------
                          Table 4.2 Input Data Set for Problem 1
    1    Test of 2DFATMIC, flow only,  FEMWATER.EX1,  cm,g,day
    1  1  1  1  10 100  0.5  1  iitr,ichng,inter,ibuf,imod,nitrht,omeht,igeom
c ***** data set 2
00  0  0  0  10  0120 IMID,IWET,ILUMP,IOPTIM,IPNTS,IVML,NSTRh,nstrt,iquar, idis
C ***** data set 3
82                        NNP
1  40  2   0.0  0.0  0.0      X
2  40  2   5.0  0.0  0.0      COORDINATE
0   0  0   0.0  0.0  0.0
1  40  2  0.0  5.0  0.0      Z
2  40  2  0.0  5.0  0.0      COORDINATE
0  0   0  0.0  0.0   0.0
c ***** data set 4
40  1  0                  NEL  NMAT  KCP
1124311  40        IE
C ****** data set  5
0                         NCM
c ***** data set 6
10  1  1                  NMPPM NCC  IRXN
0.00 10.0 0.0  l.OdO  7.316D12   0.36    0.36   3.6D-5   948.OdO 1.2    PROP
0.5                                    DKD
0 . 5                                    TRANC
1.0                                DINTS
0 . 0                                RHOMU
C ***** data set 7
0  4                                   KSP  NSPPM
0.15  0.45  0.0  -1.0d2                THKCAH
0.0   0.0   0.0   0.0  0.0             AKPROP
C ***** data set 8
44 3                                   NTI,NDTCHG
0.05  0.2  l.ODO  20.OdO               DELT,CHNG,DELMAX,TMAX
555050500050005005000055505050005000500500005  KPR
111010100010001001000011101010001000100100001  KDSK
1.Odl  2.Odl  1.OD3 8                   TDTCH
c ***** data set 9
20  50  300  Oil                   NCYL,NITER,NPITER,KSTR,KSS,KGRAV
2.0D-2  2.0D-2  1.0   0.5  0.5  0.0      TOLAH,TOLBH,WF,OMEH,OMIH,cnstkr
c ***** data set 10
11    1   0.00.00.0              I.C.
3  77   1  -97.0  0.0  0.0
81  1   1  -90.0  0.0  0.0
0  0   0    0.0   0.0  0.0
c ***** data set 11
0000                            NSELH,NSPRH,NSDPH,KSAIH
c ***** data set 12
0000                            NWNPH,NWPRH,NWDPH,KWAIH
c ***** data set 13
1   2   140                       NRES,NRNP,NRPR,NRDP,KRAIH
0.0  5.0  10.0  5.0   10.001  -5.0  1.0D38   -5.0    TRF,RF OF PROFILE 1
11   1   81   1                    IRTYP(1,2)
0000    0
111    10                    IRTYP(1,1)
000    00
11   10.ODO  0.ODO  0.0
000  O.ODO  O.ODO  0.0           RSVAB(1,1)  	 HCON
111  -90.ODO  O.ODO  0.0
000  O.ODO  O.ODO  0.0           RSVAB(1,2)  	 HMIN
10  40  120   Oil      MI,NSEQ,M,IS1,IS2,MIAD,MAD,IS1AD,IS2AD
00   0000   000
c ***** data set 14
2120                         NDNPH,NDPRH,NDDPH,KDAIH
0.0  0.0  1.0D38  0.0                 THDBF,HDBF
1    1111                      IDTYP(1,2)
0    0000
1    1110                      IDTYP(1,1)
0    0000



                                         83

-------
c ****** data  set 15
00000                           NCES,NCNP,NCPR,NCDP,KCAIH
C ****** data  set 16
00000                           NNES,NNNP,NNPR,NNDP,KNAIH
0
4.2 Problem No. 2:  One-Dimensional Column Transport


Statement of Problem No. 2

       A simple problem is presented here to illustrate the application of this document.  This is a one-

dimensional solute transport and biotransfbrmation between x = 0 and x = 50.0 cm (Figure 4.2).  Initially,

the concentrations of microbe #1 (Q) and microbe #2  (C2) are zero, and those  of microbe #3  (C },

substrate (Cs), oxygen (C0), nitrate  (Cn) and nutrient (Cp) are 3.64xlQ-3, l.OxlO'3, 2.0xlQ-3, S.OxlO'3 and

1.0 xl02mg/cm3 throughout the region of interest. This case is to demonstrate that multiple electron

acceptor respiration can occur and the simulation is done with excess nutrient.  The concentration at x =

0.0 is maintained at Q = C2 = 0, C3 = 3.64x 10'3, Cs = 2.0xlQ-2, C0 = 2.0xlQ-3, Cn  = 5.0x 10'3, and Cp =

l.Ox 102 mg/cm3.  The natural condition of zero gradient flux is imposed at x = 50.0.

       A bulk density of 1.2, a dispersivity of 0.5, an effective porosity of 0.4 (not  used in the program)

are assumed. A specific discharge (Darcy velocity) of 10.0 is assumed and a moisture content of 1.0 is

used. The retardation factors of microbe #1, microbe #2, microbe #3, substrate, oxygen, nitrate, and

nutrient are 105, 105,  105, 1.1, 1.0, 1.0 and 2.2, respectively.

       For numerical simulation the region is divided into 100 elements of equal size of 0.5. A time- step

size of 0.1 is used and 80 time-step simulation is made.  For this discretization, mesh Peclet number is

Pe = 1 and Courant number Cr = 2.
                                              84

-------
246 200 2(
T
1 cm

CO

(V)



(s>
135 199 2(



                  Figure 4.2 Region of Interest for Example 2
      To execute the problem, the maximum control-integers in the MAIN must be specified as

   PARAMETER(MAXELK=100,MAXNPK=202,MXBESK=202,MXBNPK=202,MAXBWK=11)
   PARAMETER(MXJBDK= 11 ,MXKBDK=4,MXADNK=MAXNPK+3 000)
   PARAMETER(MXMATK= 1 ,MXSPPK=4,MXMPPK= 10)
   PARAMETER(MXNTIK= 80,MXNDTK=3)

   PARAMETER(MXSELF=1 ,MXSPRF=1 ,MXSDPF=1 ,MXWNPF=1 ,MXWPRF=1 ,MXWDPF= 1)
   PARAMETER(MXCNPF=1 ,MXCESF= 1 ,MXCPRF=1 ,MXCDPF=1)
   PARAMETER(MXNNPF= 1 ,MXNESF= 1 ,MXNPRF= 1 ,MXNDPF=1)
   PARAMETER(MXRESF=1 ,MXRNPF= 1 ,MXRPRF=1 ,MXRDPF=1)
   PARAMETER(MXDNPF= 1 ,MXDPRF=1 ,MXDDPF= 1)

   PARAMETER(MXSELT= 1 ,MXSPRT= 1 ,MXSDPT= 1 ,MXWNPT= 1 ,MXWPRT= 1 ,MXWDPT= 1)
   PARAMETER(MXCNPT= 1 ,MXCEST= 1 ,MXCPRT= 1 ,MXCDPT= 1)
   PARAMETER(MXNNPT= 1 ,MXNEST= 1 ,MXNPRT= 1 ,MXNDPT= 1)
   PARAMETER(MXVEST= 1 ,MXVNPT=2,MXVPRT= 1 ,MXVDPT=2)
   PARAMETER(MXDNPT=2,MXDPRT=6,MXDDPT=2)

   PARAMETER(MXNCCK=7)

   PARAMETER(MXKGLDK=3 000,MXLS VK= 1 ,MXMSVK= 1 ,MXNDBK= 1200)
   PARAMETER(MXNEPK= 1 ,MXEPWK= 1)
   PARAMETER(MXNPWK=3 6,MXELWK=25 ,MXNPWS=3 6,MXELWS=25)
   PARAMETER(MXNPFGK=35000,MXKGLK=14000)
Input and Ouput for Problem No. 2

      Table 4.3 lists the input parameters and Table 4.4 shows the input data set for this example. To

save space, the output isavailable in electronic form.
                                     85

-------
Table 4.3  The list of input parameters for Problem 2
Parameters
number of points
AX
AZ
KA1
KA2
Kd3
Kds
Kdo
Kdn
KdD
Pb
OL
Mo0)
Un®
Mo(3)
Un*3'
Y 0)
Y (2)
A n
Y (3)
1 0
Y (3)
-*- n
Kso(1)
K (2)
Avsn
Kso(3)
Ksn(3)
K0(1)
IC2)
Notation in the
data input guide
NNP
XAD
ZAD
RKD(l)
RKD(2)
RKD(3)
RKD(4)
RKD(5)
RKD(6)
RKD(7)
PROPt(l,10)
PROPt(l,6)
GRATE(l)
GRATE(2)
GRATE(3)
GRATE(4)
YCOEFF(l)
YCOEFF(2)
YCOEFF(3)
YCOEFF(4)
RTARDS(l)
RTARDS(2)
RTARDS(3)
RTARDS(4)
RTARDO(l)
RTARDO(2)
Value
202
0.5
1
105
105
105
0.1
0
0
1.2
1.0
0.5
0.0
0.0
0.62
0.58
0.45
0.5
0.45
0.5
0.04
0.04
0.04
0.04
7.7xlQ-4
2.6xlQ-3
Unit
dimensionless
cm
cm
cmVmg
cmVmg
cmVmg
cm3/mg
cmVmg
cm3/mg
cm3/mg
mg/cm3
cm
I/day
I/day
I/day
I/day
mg/mg
mg/mg
mg/mg
mg/mg
mg/cm3
mg/cm3
mg/cm3
mg/cm3
mg/cm3
mg/cm3
Data set
3. A.
3.B.
3.B.
6. C.
6. C.
6. C.
6. C.
6. C.
6. C.
6. C.
6. B.
6. B.
19. A.
19. A.
19. A.
19. A.
19. B.
19. B.
19. B.
19. B.
19. C.
19. C.
19. C.
19. C.
19. D.
19. D.
                       86

-------
K0(3)
K(3)
Avn
K (1)
DO
K (2)
-*-*-Dn
K (3)
DO
K (3)
A^-on
v V
1 0
v ^
I n
v (3)
I 0
v (3)
In
a0(1)
CCn(2)
a0(3)
an(3)
i (i)
Ao
V2)
V3)
1 (3)
An
-P (1)
-1 o
-P (2)
x n
Y (3)
r (3)
-1 n
e0(1)
en(2)
e0(3)
en(3)
Kc
maximum time-step size
no. of times to reset
time-step size
RTARDO(3)
RTARDO(4)
RTARDN(l)
RTARDN(2)
RTARDN(3)
RTARDN(4)
SCOEFF(l)
SCOEFF(2)
SCOEFF(3)
SCOEFF(4)
ECOEFF(l)
ECOEFF(2)
ECOEFF(3)
ECOEFF(4)
DCOEFF(l)
DCOEFF(2)
DCOEFF(3)
DCOEFF(4)
SATURC(l)
SATURC(2)
SATURC(3)
SATURC(4)
PCOEFF(l)
PCOEFF(2)
PCOEFF(3)
PCOEFF(4)
COFK
DELMAX
NDTCHG
7.7xlO-4
2.6xlQ-3
i.oxio-3
i.oxio-3
i.oxio-3
i.oxio-3
1.4
2.2
1.4
2.2
0.0402
0.1
0.0402
0.1
0.02
0.02
0.004
0.004
v.vxio-4
2.6xlQ-3
7.7xlQ-4
2.6xlQ-3
0.122
0.122
0.122
0.122
l.lxlQ-4
0.1
0
mg/cm3
mg/cm3
mg/cm3
mg/cm3
mg/cm3
mg/cm3
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
I/day
I/day
I/day
I/day
mg/cm3
mg/cm3
mg/cm3
mg/cm3
dimensionless
dimensionless
dimensionless
dimensionless
mg/cm3
day
dimensionless
19. D.
19. D.
19. E.
19. E.
19. E.
19. E.
19. F.
19. F.
19. F.
19. F.
19. G.
19. G.
19. G.
19. G.
19. H.
19. H.
19. H.
19. H.
19.1.
19.1.
19.1.
19.1.
19. J.
19. J.
19. J.
19. J.
19. K.
8.B.
8. A.
87

-------
total simulation time
no. of time-steps
tolerance for nonlinear
iteration
relaxation factor for
nonlinear iteration
TMAX
NTI
TOLBt
OMEt
8
80
IxlO'4
1.0
day
dimensionless
dimensionless
dimensionless
8.B.
8. A.
17. B.
17. B.
                             Table 4.4 Input Data Set for Problem 2
    1 adv-dis, Multi-electron  acceptor,  excess nutrient, VB, cm, mg,  day
    11111 100  0.5 1  iitr,ichng,inter,ibuf,imod,nitrht,omeht,igeom
c ***** DATA SET 2  : OPTION  PARAMETERS FOR BOTH FLOW AND TRANSPORT
00  1  1  0  000 12 0 IMID,IWET,ILUMP,IOPTIM,IPNTS,IVML,NSTRH,NSTRT,iquar,idis
c ***** DATA SET
202
                     NODAL  COORDINATE
                          NNP
                          0 .ODO
                          0 . ODO
                         0 .  ODO
                         0.  ODO
                         0.  ODO
1 100  2   O.ODO
2 100  2   O.ODO
000   O.ODO
1 100  2   O.ODO
2 100  2   l.ODO
0  0   0  0.0  0.0
c ***** DATA SET 4
100 1  0
113   4   21
c ***** DATA SET 5
0
C ***** DATA SET 6
10  7  1
0.00 10.0 0.0
1.0d5  1.0d5  1.0d5
0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  0.0  0.0
c ****** DATA SET 7
0  4
0.15  0.45  0.0  -1.0d2
0.0   0.0   0.0   0.0  0.0
c ***** DATA SET 8
80 1
0.1DO  0.0  9.0D1
 0.5DO  O.ODO      X
 0.5DO  O.ODO      X
0 . ODO
0 .ODO
O.ODO  O.ODO      COORDINATE
   0 . 0
   ELEMENT INCIDENCES
        NEL   NMAT  KCP
  1  100       IE
   MATERIAL  TYPE  CORRECTION
        NCM
   MATERIAL  PROPERTIES
        NMPPM  NCC  IRXN
0 . 5
0 . 0
.0
1. 0
0 . 0
0.0 0.0
1 . 2
TRANC
DINTS
RHOMU
9 .48d-5
DKD



                                                        1 . 0  PROP
555550 0 0 000 00 00
5
1 0000 0 0 000 00 00
1
              1.OD38
c ***** DATA SET 17  :
50   999 11   -11
l.Od-3  l.Od-4  1.0
1105151
0.0005  0.0005
c ***** DATA SET 18  :
               1.0d3 7.316D12
                      0.1   0.0
                      0  0.0  C
                      .0  1.0
                     0.0  0.0
                      SOIL  PROPERTIES
                                       KSP   NSPPM
                                       SPP(J,I,1)
                                       SPP(J,I,2)
                     TIME CONTROL  PARAMETERS
                                       NTI,NDTCHG
                     9O.ODO             DELT,CHNG,DELMAX,TMAX
                                        50000
   201
    0
   201
    0
   201
       1  O.ODO   0.0
       0  0.0  0.0  0.0
        1  0.OdO    0.0
       0  0.0  0.0  0.0
        1 3.64d-8 0.0  0.0
                      10000

                     TDTCH
    BASIC REAL AND  INTEGER PARAMETERS FOR TRANSPROT
              nitert,npitert,kstrt,ksst,kvit,miconf
   0.0  l.OdO  1.OdO   1.0     tolat,tolbt,wt,wv,omet,omit,aphag
    112   IZOOM,IDZOOM,IEPC,NXA,NZA,NXD,NZD,NXW,NZW,IDETQ
                    ADPEPS, ADPARM
    INITIAL/PRE-INITIAL CONDITIONS FOR TRANSPORT
     0.0               1C FOR Microbe #1
       0 . 0
  1C FOR microbe #2

1C FOR Microbe #3
    0  0  0.0  0.0  0.0

-------
1  201  1  1 .OD-3  0.0  0.0
0   0  0  0.0   0.0  0.0
1  201  1  2.OD-3  0.0  0.0
0   0  0  0.0   0.0  0.0
1  201  1  5.Od-3  0.0  0.0
0   0  0  0.0   0.0  0.0
1  201  1  1.0D2  0.0  0.0
0   0  0  0.0   0.0  0.0
c ***** DATA SET 19
 0.0  0.0  0.620.58
 0.45  0.50.450.5
4.0D-2  4.0D-2  4.0D-2
7.7D-4  2.6D-3  7.7D-4
l.Od-3  l.OD-3  l.OD-3
1.4  2.2    1.4   2.2
4.02D-2 l.OD-1 4.02D-2
                                          1C  FOR Substrate
                    1C  FOR  Oxygen
                    1C  FOR  Nitrate
                   1C  FOR  Nutrient
0.02 0.02 0.004
7.7D-4  2.6D-3
1.22D-1 1.22D-1
  1.ld-4
MICROBE-CHEMICAL INTERACTION  CONSTANTS
                  GRATE
                YCOEFF
  4.0D-2        RTARDS
  2.6D-3        RTARDO
  l.OD-3        RTARDN
                SCOEFF
                ECOEFF
                DCOEFF
              SATURC
       1.OD-1
  0 . 004
 7.7D-4 2.6D-3
1.22D-1 1.22D-1
  Kso,  Ksn
  Ko,  Kn
  Kpo,  Kpn
 gammao,  gamman
  alphao, alphan
  lambdao, lambdan
GAMMAo,  GAMMAn
c
0
c
0
c
1
0.
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
c
2
0 .
0 .
0.
0.
0 .
0 .
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
* * *
0
** *
0
* * *
2
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
** *
6
0
0
0
0
0
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
** DATA
0 0

** DATA
0 0

** DATA
1 2
0 .0
1
0
1
0
1
0
1
0
1
0
1
0
1
0 0
1 201
0 0
100 1
0 0
0
1
1
0
1
0
1
0
1
0
1
0
1
0
1
0


2
0
** DATA
2 0
2 . OD-2
2 . OD-3
3 .64d-6
1 .002
5 . Od-3
0 . 0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0

1
1

1

1
6
0
6
0
3
0
1
0
2
0
5
0
4
0
1
0
SET 20 : ELEMENT SOURCE

SET 21 : WELL SOURCE/SI

SET 22 : VARIABLE B.C.

.Od38 O.OdO
0
0
0
0
0
0
0
0
0
0
0
0
0

1
0
11 11
0000
SET 23 : DIRICHLET B.C.

.OD38 2. OD-2
.OD38 2. OD-3
1.0D38 3.64d-8
.0038 1.0D2
1.0d38 5. Od-3
.Od38 O.OdO
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
                PCOEFF   Epsilon
               COFK
               /SINK
               NSEL,NSPR,NSDP,ksai
               \TK
               NWNP,NWPR,NSDP,kwai

               NVES,NVNP,NVPR,NVDP,kdai

               IvTYP FOR Microbe #1

               IVTYP for Microbe #2

               IVTYP for Microbe #3

               IVTYP for Substrate

               IVTYP for Oxygen

               IVTYP for Nitrate

               IVTYP for Nutrient

               NPVB

               ISV
                                     DIRICHLET BOUNDARY CONDITION
                                     TCDBF  CDBF    Profile 1
                                                    Profile  2
                                                    Profile  3
                                                    Profile  4
                                                    Profile  5
                                                    Profile  6
                                     IDTYP for Microbe #1

                                     IDTYP for Microbe #2

                                     IDTYP for Microbe #3

                                     IDTYP for substrate

                                     IDTYP for oxygen

                                     IDTYP for nitrate

                                     IDTYP for nutrient

                                                   NPDB
                                            89

-------
                        CAUCHY B.C.
                        NEUMANN B.C.
c *****  DATA SET  24
0   0000
c *****  DATA SET  25
0   0000
c *****  DATA SET  26  :  HYDROLOGICAL VARIABLES
1 201  1  10.0  0.0   0.0                 VELOC(NNP,1)
0   0  0   0.0  0.0   0.0
1 201  1    0.0   0.0   0.0                VELOC(NNP,2)
0  0   0   0.0   0.0   0.0
1 99   1   l.OdO   0.0                     TH(1..4,NEL)
0  0   0   0.0   0.0
0
nces,ncnp,ncpr,ncdp,kcai

nnes,nnnp,nnpr,nndp,knai
Simulation Results

The results at t = 4 days and t = 8 days are shown as follows. These results are similar to those obtained

in Widdowson, et. al. (1988). Thus, the transport module of 2DFATMIC is indirectly verified.
                                                                t = 8 days
                     20      30
                       Distance
         Figure 4.3 The concentration profiles of biomass, substrate, oxygen, and nitrate at
                                   t = 4 days and t = 8 days
                                             90

-------
4.3  Problem No. 3:  Two-Dimensional Coupled Flow and Multicomponent Transport





       Problem







Statement of Problem No. 3




       This problem is presented in "Denitrification in nonhomogeneous laboratory scale aquifers: 5:




user's manual  for the mathematical model LT3VSI" by G.A. Bachelor et al. reported in 1990.  The




example aquifer used for this problem is 1.4 meters long in the X direction, 1.6 meters thick in the Z




direction, and shown in Figure 4.4.  This aquifer has 8 different materials, one injection well at (0.1,0.1)




and one extraction well at  (1.3,1.5).  The hydrological and microbial dynamical data are all from this




report. The initial condition of the flow field is obtained by simulating steady state of flow field without




sources and sinks.  Then the flow field and concentration distribution are updated at each time step.




       There are two types of microbes included in the system, say microbe #1 and microbe #3 with




1.7?xlO~4  Kg/m3 initially.  The initial concentrations of the chemicals  are 5><10"3 Kg/m3 of substrate,




5 x 10"3 Kg/m3 of oxygen, 5 x 10"3 Kg/m3  of nitrate, and 3 x 10"3 Kg/m3 of nutrient. The daily inj ection and




withdrawal rates of water are 3.75xlO~3 and 3.75xlO~3 m3, respectively.  The total hydraulic head is 1.0




m at the upstream boundary AB (Figure 4.4) and 0.0 m at the downstream boundary CD (Figure 4.4). For




transport simulation, variable boundary condition is implemented at the downstream boundary CD (Figure




4.4) and l.TTxlO'4 Kg/m3 of microbe #1, l.TTxlO'4 Kg/m3 of  microbe #3, l.SxlO'2 Kg/m3 of substrate,




5.Ox 10"3 Kg/m3 of oxygen, 5.Ox 10"3 Kg/m3 of nitrate, and 3.Ox 10"3 Kg/m3 of nutrient influents through




the upstream boundary.  Because microbe #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 4 day




simulation.
                                             91

-------
               No flux for both flow and transport
B l.o
1.2
1.0
h= 1.0
0.6
0.4
Ann
(1)
(6)
J3)
(0.1,0.1)
(5)
(4)
(7)
(6)
(5)
(7)
(8)

•
(1.3,1.5)
(3)
(4)
(2)
                                                            D
                                                           h  = 0.0 for flow
                                                           Variable
                                                           boundary for
                                                           transport
        0.0
0.3    0.5
 0.9    1.1       1.4
No flux for both flow and transport
         (1) : material type 1

             Figure 4.4 The region of Problem 3 with 8 different material properties



      For numerical simulation the region is divided into 14 x 16 = 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 80 time-step simulation is made.

      To execute the problem, the maximum control-integers in the MAIN must be specified as

PARAMETER(MAXELK=224,MAXNPK=255,MXBESK=60,MXBNPK=60,MAXBWK=33)
PARAMETER(MXJBDK= 17,MXKBDK=4,MXADNK=MAXNPK+4200)
PARAMETER(MXMATK=8,MXSPPK=2,MXMPPK= 10)
PARAMETER(MXNTIK=8 0,MXNDTK=3)
                                       92

-------
PARAMETER(MXSELF= 1 ,MXSPRF=1 ,MXSDPF= 1 ,MXWNPF=2,MXWPRF=2,MXWDPF=3)
PARAMETER(MXCNPF=1 ,MXCESF= 1 ,MXCPRF=1 ,MXCDPF= 1)
PARAMETER(MXNNPF= 1 ,MXNESF= 1 ,MXNPRF= 1 ,MXNDPF= 1)
PARAMETER(MXRESF= 1 ,MXRNPF= 1 ,MXRPRF=1 ,MXRDPF= 1)
PARAMETER(MXDNPF=34,MXDPRF=2,MXDDPF=2)

PARAMETER(MXSELT= 1 ,MXSPRT= 1 ,MXSDPT= 1 ,MXWNPT=2,MXWPRT=2,MXWDPT=3)
PARAMETER(MXCNPT= 1 ,MXCEST= 1 ,MXCPRT= 1 ,MXCDPT= 1)
PARAMETER(MXNNPT= 1 ,MXNEST= 1 ,MXNPRT= 1 ,MXNDPT= 1)
PARAMETER(MXVEST= 16,MXVNPT= 17,MXVPRT= 1 ,MXVDPT=2)
PARAMETER(MXDNPT= 17,MXDPRT=5 ,MXDDPT=2)

PARAMETER(MXNCCK=7)

PARAMETER(MXKGLDK=4200,MXLSVK=700,MXMSVK=150,MXNDBK=500)
PARAMETER(MXNEPK= 1 ,MXEPWK= 1)
PARAMETER(MXNPWK=3 6,MXELWK=25 ,MXNPWS=3 6,MXELWS=25)
PARAMETER(MXNPFGK=7100,MXKGLK= 175 00)
Input and Ouput for Problem No. 3

      Table 4.5 shows the input parameters used for this example. Table 4.6 shows the input data set

for the transient simulation of the coupled problem. The output results are given isavailable in electronic

form.
                  Table 4.5 The list of input parameters for Problem 3
Parameters
number of points
AX
AZ
no. of materials
Mo0)
^(2)
Mo(3)
u/3)
Notation in the
data input guide
NNP
XAD
ZAD
NMAT
GRATE(l)
GRATE(2)
GRATE(3)
GRATE(4)
Value
510
0.1
0.1
8
4.0
0.0
4.0
2.5
Unit
dimensionless
m
m
dimensionless
I/day
I/day
I/day
1/dav
Data set
3. A.
3.B.
3.B.
4. A.
19. A.
19. A.
19. A.
19. A.
                                     93

-------
Y 0)
Y (2)
A n
Y (3)
1 0
Y (3)
-*- n
Kso(1)
K (2)
Avsn
Kso(3)
Ksn(3)
K0(1)
K(2)
Avn
K(3)
0
K(3)
n
K (1)
DO
K (2)
A^-on
K (3)
DO
K (3)
-"-^DII
v (1)
I 0
v (2)
In
v &
1 0
v <3>
I n
a0(1)
an(2)
a0(3)
an(3)
AO(I)
1 (2)
An
1 (3)
Ao
1 (3)
An
J1 (1)
YCOEFF(l)
YCOEFF(2)
YCOEFF(3)
YCOEFF(4)
RTARDS(l)
RTARDS(2)
RTARDS(3)
RTARDS(4)
RTARDO(l)
RTARDO(2)
RTARDO(3)
RTARDO(4)
RTARDN(l)
RTARDN(2)
RTARDN(3)
RTARDN(4)
SCOEFF(l)
SCOEFF(2)
SCOEFF(3)
SCOEFF(4)
ECOEFF(l)
ECOEFF(2)
ECOEFF(3)
ECOEFF(4)
DCOEFF(l)
DCOEFF(2)
DCOEFF(3)
DCOEFF(4)
SATURC(l)
0.4
0.17
0.4
0.17
0.018
0.018
0.018
0.018
3.0xlQ-5
2.0xlQ-5
3.0xlQ-5
2.0xlQ-5
3.0xlQ-4
3.0xlQ-4
3.0xlQ-4
3.0xlQ-4
1.0
0.375
1.0
0.375
0.004
0.002
0.004
0.002
0.02
0.02
0.02
0.02
3.0xlQ-5
Kg/Kg
Kg/Kg
Kg/Kg
Kg/Kg
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
Kg/m3
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
I/day
I/day
I/day
I/day
Kg/m3
19. B.
19. B.
19. B.
19. B.
19. C.
19. C.
19. C.
19. C.
19. D.
19. D.
19. D.
19. D.
19. E.
19. E.
19. E.
19. E.
19. F.
19. F.
19. F.
19. F.
19. G.
19. G.
19. G.
19. G.
19. H.
19. H.
19. H.
19. H.
19.1.
94

-------
rn(2)
-P (3)
-1 o
-P (3)
-1 n
e0(1)
en(2)
e0(3)
e (3)
^n
Kc
no. of elements
no. of materials to be
corrected
steady-state for flow
transient-state for
transport
initial time-step size
time-step size increment
percentage
maximum time-step size
no. of times to reset
time-step size
total simulation time
no. of time-steps
tolerance for flow
nonlinear iteration
relaxation factor for
flow nonlinear iteration
tolerance for transport
nonlinear iteration
relaxation factor for
transport nonlinear
iteration
Pw
U,»,
SATURC(2)
SATURC(3)
SATURC(4)
PCOEFF(l)
PCOEFF(2)
PCOEFF(3)
PCOEFF(4)
COFK
NEL
NCM
KSSh
KSSt
DELT
CHNG
DELMAX
NDTCHG
TMAX
NTI
TOLAh
OMEh
TOLBt
OMEt
PROP(I,4)
PROP(I,9)
2.0xlQ-5
s.oxio-5
2.0xlO-5
0.05
0.021
0.05
0.021
l.lxlO'4
224
134
0
1
0.05
0
0.05
0
4
80
IxlQ-2
1.0
IxlO'4
1.0
1000.0
948.3264
Kg/m3
Kg/m3
Kg/m3
dimensionless
dimensionless
dimensionless
dimensionless
Kg/m3
dimensionless
dimensionless
dimensionless
dimensionless
day
dimensionless
day
dimensionless
day
dimensionless
m
dimensionless
dimensionless
dimensionless
Kg/m3
Kg/m/day
19.1.
19.1.
19.1.
19. J.
19. J.
19. J.
19. J.
19. K.
4. A.
5. A.
9. A.
17. A.
8.B.
8.B.
8.B.
8. A.
8.B.
8. A.
9. B.
9. B
17. B.
17. B.
6. B.
6. B.
95

-------
1 g
PROP(I,5)
7.316xl010
m/day2
6.B. |
                        Table 4.6 Input Data Set for Problem 3
  1 advection-dispersion, flow & transport,multi-materials,  Kg,  m,  day
  1  1  1  1  11   100  0.5 0 iitr,ichng,inter,ibuf,imod,nitrht,omeht,igeom
***** DATA SET 8: CONTROL PARAMETER
c ***** DATA SET 2
255
1 16 15 0.0 0.0
2 16 15 0.1 0.0
3 16 15 0.2 0.0
4 16 15 0.3 0.0
5 16 15 0.4 0.0
6 16 15 0.5 0.0
7 16 15 0.6 0 .
8 16 15 0.7 0 .
9 16 15 0.8 0
10 16 15 0.9 0.
11 16 15 1.0 0 .
12 16 15 1.1 0 .
13 16 15 1.2 0.
14 16 15 1.3 0.
15 16 15 1.4 0 .
0 0 0 0.0 0.0
1 16 15 0.0 0.1
2 16 15 0.0 0.1
3 16 15 0.0 0.1
4 16 15 0.0 0.1
5 16 15 0.0 0.1
6 16 15 0.0 0.1
7 16 15 0.0 0 .
8 16 15 0.0 0 .
9 16 15 0.0 0
10 16 15 0.0 0.
11 16 15 0.0 0 .
12 16 15 0.0 0 .
13 16 15 0.0 0.
14 16 15 0.0 0.
15 16 15 0.0 0 .
0 0 0 0.0 0.0
c ***** DATA SET 3 :
224 8 0
1 1 2 17 16 1
c ***** DATA SET 4 :
134
1 3 14 3 0
2 3 14 3 0
3 3 14 3 0
68120
20 8 1 2 0
34 8 12 0
48 8 12 0
68 5 14 2 0
69 5 14 2 0
70 5 14 2 0
180 3 14 3 0
181 3 14 3 0
182 3 14 3 0
152 2140
166 2140
: NODAL COORDINATE
NNP
0. 0
0. 0
0 . 0
0 . 0
0. 0
0. 0
0 0.0
0 0.0
. 0 0.0
0 0.0
0 0.0
0 0.0
0 0.0
0 0.0
0 0.0
0 . 0
0. 0
0. 0
0 . 0
0 . 0
0. 0
0. 0
1 0.0
1 0.0
. 1 0.0
1 0.0
1 0.0
1 0.0
1 0.0
1 0.0
1 0.0
0 . 0

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z

ELEMENT INDICES
NEL
14 16
NMAT KCP

MATERIAL CORRECTION
NCM































                                       96

-------
4
5
57
71
60
74
62
76
66
80
94
95
153
167
181
182
0
3
3
2
2
1
1
3
3
1
1
3
3
1
1
3
3
0
14
14
1
1
1
1
1
1
1
1
14
14
1
1
14
14
0
4
4
6
6
5
5
7
7
8
8
7
7
5
5
6
6
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
c *****
10   7
l.OD-2 1
0.0  0.0
1.Od-2
l.ODO  1
0.0  0.0
0 . ODO
l.OD-1 1
0.0  0.0
1.Od-4
3.16D-1
0.0  0.0
0. OdO
5.62D-2
0.0  0.0
0. OdO
3.16D-2
0.0  0.0
1.Od-3
l.OD-1 1
0.0  0.0
0 . OdO
3.16D-1
0.0  0.0
0. OdO
   1 . Od3
   0 . 0
p *****
DATA SET 5: MATERIAL  PROPERTIES
  1               NMPPM  NCC  IRXN
.OD-2 0.0  1.0d3  7.316dlO  O.Od-2  O.Od-2 O.OOd-5 0.3 1414.!
  0.0  0.0  0.0   0.0  0.0      DKD
PROP  MAT 1
                                                  TRANC
                                                0.3  1901.9

                                                  TRANC
                                                 0.3  1689. 5

                                                  TRANC
PROP  MAT 2
 PROP  MAT 3
l.Od-2  l.Od-2  l.Od-2   l.Od-2   l.Od-2   l.Od-2
.ODO  0.0  1.0d3 7.316dlO  0.OdO   0.OdO   0.ODO
  0.0  0.0  0.0  0.0  0.0      DKD
O.OdO   O.OdO   O.OdO    0.OdO    0.OdO    0.OdO
.OD-1 0.0 1.0d3 7.316dlO O.OD-3  O.OD-3  O.OOD-5
  0.0  0.0  0.0  0.0  0.0      DKD
l.Od-4  l.Od-4  l.Od-4   l.Od-4   l.Od-4   l.Od-4
3.16D-1 0.0 1.0D3  7.316dlO O.OOD-3  O.OOD-3 O.OOD-4 0.3 1755.8  PROP  MAT  4
  0.0  0.0  0.0  0.0  0.0      DKD
O.OdO   O.OdO   O.OdO    O.OdO    O.OdO    O.OdO    TRANC
5.62D-2 0.0 1.0D3  7.316dlO O.OOD-3  O.OOD-3 O.OOD-5 0.3 1472.8  PROP  MAT  5
  0.0  0.0  0.0  0.0  0.0      DKD
O.OdO   O.OdO   O.OdO    O.OdO    O.OdO    O.OdO    TRANC
3.16D-2 0.0 1.0D3  7.316dlO O.OOD-3  O.OOD-3 O.OOD-5 0.3 1515.8  PROP  MAT  6
  0.0  0.0  0.0  0.0  0.0      DKD
l.Od-3  l.Od-3  l.Od-3   l.Od-3   l.Od-3   l.Od-3   TRANC
.OD-1 0.0 1.0D3 7.316dlO O.OOD-3 O.OOD-3  O.OOD-5 0.3 1512.4 PROP MAT 7
  0.0  0.0  0.0  0.0  0.0      DKD
O.OdO   O.OdO   O.OdO    O.OdO    O.OdO    O.OdO    TRANC
3.16D-1 0.0 1.0D3  7.316dlO O.OOD-3  O.OOD-3 O.OOD-4 0.3 1706.1  PROP  MAT  8
  0.0  0.0  0.0  0.0  0.0      DKD
O.OdO   O.OdO   O.OdO    O.OdO    O.OdO    O.OdO    TRANC
  1.0d3  1.Od3  1.0d3  1.Od3  1.Od3  1.Od3  DINTS
0.0  0.0  0.0  0.0  0.0  0.0  RHOMU
DATA SET 6:  MATERIAL PROPERTIES
                              KSP  NSPPM
1 2
-1000 . 0
-1000 . 0
-1000 .0
-1000 .0
-1000 . 0
-1000 . 0
-1000 .0
-1000 .0
0.465
0.285
0.365
0.323
0.387
0 . 412
0.364
0.322
1 . 0
1 . 0
1. 0
1. 0
1 . 0
1 . 0

1000 . 0
1000 . 0
1000 . 0
1000 . 0
1000 . 0
1000 . 0
1000 . 0
1000 . 0
0.465
0.285
0.365
0.323
0.387
0 . 412
0.364
0.322
1 . 0
1 . 0
1. 0
1. 0
1 . 0
1 . 0
                                  97

-------
  1.0     1.0
  1.0     1.0
  0.0     0.0
  0.0     0.0
  0.0     0.0
  0.0     0.0
  0.0     0.0
  0.0     0.0
  0.0     0.0
  0.0     0.0
c ***** DATA SET 7: TIME  CONTROL PARAMETERS
80  1                                   NTI,NDTCHG
0.05  0.0  l.ODO 4.0               DELT,CHNG,DELMAX,TMAX
52 0 0 0 05 00 00 0050    0 0 0500055000 5000      555
55555555
00 0 0 0 00 00 00 0000    0 0 0000000000 0000      000
00000000
1.Od5  2.Od4  1.OD3 8                   TDTCH
c ***** DATA SET 9: CONTROL  PARAMETERS FOR FLOW
20  100  500  101                  NCYL,NITERH,NPITERH,KSTRH,KSSH,KGRAV
l.OD-3  l.OD-3  1.0 1.0   1.0   0.0         TOLAH,TOLBH,WF,OMEH,OMIH,cnstkr
c ***** DATA SET 10:  I.C.  FOR  FLOW
1  254  10.00.00.0               I.C. OF H
0   0   0   0.0  0.0  0.0
c ***** DATA SET 11:  ELEMENTAL SOURCES FOR FLOW
0000                             NSELH,NSPRH,NSDPH,KSAIH
c ***** DATA SET 12:  WELL SOURCES  FOR FLOW
2  2
0 . 0
0. 0
 3
0 . 0
0 .0
 0
 0
 0
 1
 0
                                       NWNPH,NWPRH,NWDPH,KWAIH
                                   7.5d-3
                                   -4.Od-3
                                       IWTYPH (1,2)
                                       IWTYPH (2,2!
                    TSOSFH,SOSFH
000
c *****
34  2  2
0.0  1.0
0.0  0.0
1161
18  16
0   0
1161
18 16  1
000
0
   0  0
  ** ** *
100  500
1.Od-3
1  105
0 . 0005
C *****
   254 1
    0  0
   254
    0  0
   254
    0  0
   254
    0  0
   254
    0  0
   254
    0  0
 0
  l.Od-6  7.5d-3  1.0D38
  l.Od-6 -4.Od-3  1.0D38
 17  0
239  0
 0   0
 1   1
 0   0
DATA SET 13: VARAIABLE  B.C.
 0   0
DATA SET 14: DIRICHLET  B.C.
  0
  1.OD38  1.0
  1.OD38  0.0
  1  15
1  15 15
00     0
  1  0
  2  0
  0  0
DATA SET 15: CAUCHY  B.C.
 0  0
DATA SET 16: NEUMANN B.C.
 0  0                          NNESH,NNNPH,NNPRH,NNDPH,KNAIH
DATA SET 17: CONTROL PARAMETERS FOR TRANSPORT
 1120                     nitert,npitert,kstrt,ksst,kvi,miconf
l.Od-3  1.0  0.0  l.OdO   1.OdO  1.0    tolat,tolbt,wt,wv,omet,omit,alphg
 555221        izoom,idzoom,iepc,nxa,nya,nxd,nzd,nxw,nyw,idetq
0.0005   1.0                  ADPEPS, ADPARM
DATA SET 18: I.C. FOR TRANSPORT
  1.77D-4 0.00.0             1C FOR microbe #1
  0.0  0.0  0.0
10.00.00.0               1C FOR microbe #2
  0.0  0.0  0.0
1  1.77D-4  0.00.0           1C FOR microbe #3
  0.0  0.0  0.0
1  5.OD-3  0.0  0.0            1C FOR substrate
  0.0  0.0  0.0
1  5.OD-3  0.00.0            1C FOR oxygen
  0.0  0.0  0.0
1  5.OD-3  0.00.0            1C FOR nitrate
  0.0  0.0  0.0
IWTYPH(1-.2,1)


NRES,NRNP,NRPR,NRDP,KRAIH

NDNPH,NDPRH,NDDPH,KDAIH
THDBF,HDBF

IDTYPH(1..22,2)


IDTYPH(1..22,1)



NCESH,NCNPH,NCPRH,NCDPH,KCAIH
                                          98

-------
                        1.8D-2
                        2.OD-5
                        3.OD-4
1  254  1  3.OD-3  0.0  0.0
0   0  0  0.0   0.0  0.0
c ***** DATA SET 19: MICROBE-CHEMICAL
 4.0  0.0  4.0  2.5
 0.4  0.4  0.4  0.17
1.8D-2  1.8d-2  1.8D-2
3.OD-5  3.0d-5  3.OD-5
3.0d-4  3.0d-4  3.OD-4
1.0  0.375  1.0   0.375
0.004  0.002  0.004  0.002
0.02 0.02 0.02 0.02
3.OD-5  2.0d-5   3.OD-5 2.OD-5
0.05 0.021 0.05 0.021
  1.ld-4
c ***** DATA SET 20:
0000
                                       1C  FOR  nutrient

                                      INTERACTION PARAMETERS
                                     GRATE
                                 YCOEFF
                                  RTARDS
                                  RTARDO
                                  RTARDN
                                 SCOEFF
                                 ECOEFF
                                 DCOEFF
                                  SATURC
                                 PCOEFF
                                 COFK
                 ELEMENTAL SOURCES FOR TRANSPORT
                                 NSEL,NSPR,NSDP,ksai
                                       Kso,  Ksn
                                       Ko,  Kn
                                       Kpo,  Kpn
                                       gammao,  gamman
                                       alphao,  alphan
                                       lambdao,  lambdan
                                       GAMMAo,  GAMMAn
                                       Epsilon
2  2
0. 0
0 . 0
  17
   1
   0
   1
   0
   1
   0
   1
   0
   1
   0
   1
   0
   1
      3
     0 .0
     0 . 0
      239
      1  1
      0  0
    DATA SET 21:
     0
           1.Od-6
      0. 0
      0 . 0
             WELL SOURCES FOR TRANSPORT
                              NWNPH,NWPRH,NWDPH,KWAIH
               7.5d-3  0.0  1.0D38  7.5d-3  0.0
         Od-6 -4.0d-3 0.0  1.0D38  -4.0d-3  0.0
 16
 0 .0
 1
 0
 1
 0
 1
 0
 1
 0
 1
 0
 1
 0
 1
 0
 1
 0
 1
 0
c *•>
17
0. 0
0 . 0
0 . 0
0. 0
0. 0
1
0
1
0
1
0
 17
  0. 0
 15  1
 0
 15
 0
 15
 0
 15
 0
 15
 0
 15
 0
 15
 0
16
 0
15
DATA SET ;
 120
    .Od38
                     VARIABLE B.C.
           0. OdO
0
    1
    0
    1
    0
    1
    0
    1
    0
    1
    0
    1
    0
    1
    0
    15
    0
    1
    0
               0
               0
               0
               0
               0
               0
               0
               0
               0
               0
               0
               0
               0
               0
               15
               0
               2
               0
               FOR TRANSPORT
                 NVES,NVNP,NVPR,NVDP,kdai
                  Profile
                 Profile type for microbe

                 Profile type for microbe
*** DATA SET 23:
520
 1.77D-4  1.0D38
 0 . 0
 1.5D-2
 5.OD-3
 3.OD-3
    1
      1.OD38
      1.OD38
    1.OD38
    1.OD38
                          #1

                          #2
 14   1   1
  000
 DIRICHLET B.C.

  1.77D-4
  0 . 0
  1.5D-2
5.OD-3
3.OD-3
                 Profile type for microbe #3

                 Profile type for substrate

                 Profile type for oxygen

                 Profile type for nitrate

                 Profile type for nutrient
 NPDB
 VB element side

FOR TRANSPORT
 DIRICHLET BOUNDARY CONDITION
 TCDBF  CDBF    FOR microbe  #1

               FOR substrate
               FOR oxygen
               FOR nutrient
 IDTYP for microbe #1

 IDTYP for microbe #2

 IDTYP for microbe #3
                                         99

-------
1
0
1
0
1
0
1
0
1
0
16
0
16
0
16
0
16
0
16
0
1
0
1
0
1
0
1
0
1
0
3
0
4
0
4
0
5
0
1
0
0
0
0
0
0
0
0
0










15
0
                                        IDTYP  for substrate

                                        IDTYP  for oxygen

                                        IDTYP  for nitrate

                                        IDTYP  for nutrient
                                                       NPDB
c *****  DATA SET 24:  CAUCHY B.C.  FOR TRANSPORT
0
    0000
nces,ncnp,ncpr,ncdp,kcai
c *****  DATA SET 25:  NEUMANN B.C.  FOR TRANSPORT
0   0000
0
nnes,nnnp,nnpr,nndp,knai
Simulation Results

       The results of velocity fields and concentration contours simulated at t = 2 days and 4 days are

shown in Figure 4.5 and Figure 4.6.
                      TIME = 2 Days
                 TIME = 4 Days
  N
1.6-

1.4-
1.2-
1 0


ox-


.6

04-


0.2-

00-


















                                               N
               -+-—-+—-+—-••+-——h--^ * * t- -
           0.0  0.2  0.4  0.6  0.8  1.0  1.2   1.4
                           X
1.6-
1.4-
1.2-
1.0-
08-


.6 —

04-


0.2-

00-














      0.0  0.2  0.4  0.6  0.8   1.0  1.2   1.4
                      X
             •=0.30
      --=0.30
                   Figure 4.5 The flow field at (a) t = 2 days and (b) t = 4 days
                                            100

-------
The following three pages are Figure 4.6. They contain the simulation results of transport and fate of

chemicals and microbes.
                                                                Microbe 1 at Time  =  4  Days
                Microbe  1  at Time = 2 Days
        1.6


        1.4


        1.2


        1.0


        0.8


        0.6


        0.4


        0.2


       -0.0
0.0003'
                       0.00061
                          X
                                                        1.4 —
1.2-
                                     1.0 —
                                     0.2 —
                                                       -0.0-
             i    i     i    i     i    i     i     r
             0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4
                                          1    I     I     I    I     I    I     T
                                          0.0  0.2  0.4  0.6  0.8  1.0  1.2 1.4
                Microbe  3  at Time = 2 Days
                                              Microbe 3 at Time = 4 Days
        1.4 —
        1.2-
        1.0 —
        0.2 —
       -0.0-
                                     1.6 —


                                     1.4 —


                                     1.2-


                                     1.0 —


                                     0.8 —


                                     0.6-


                                     0.4-


                                     0.2 —


                                    -0.0-
             1    I     I    I     I    I     I     T
             0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4
                                          1    I     I     I    I     I    I     T
                                          0.0  0.2  0.4  0.6  0.8  1.0  1.2 1.4
                                              101

-------
        Substrate at Time  =  2  Days
    Substrate at Time = 4 Days
1.2 —
0.4 —
0.2-
     i    i     i    i     i    i     i    r
    0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4
                                               1.2 —
                                               0.4 —
                                               0.2-
i     i    i     i    i     i    i     r
0.0 0.2  0.4 0.6  0.8 1.0  1.2 1.4
          Oxygen at Time = 2 Days
     Oxygen  at  Time  =  4  Days
1.2 —
0.4 —
0.2-
     i    i     i    i     i    i     i    r
    0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4
                                               1.2 —
                                               0.4 —
                                               0.2-
i     i    i     i    i     i    i     r
0.0 0.2  0.4 0.6  0.8 1.0  1.2 1.4
                                      102

-------
         Nitrate at Time  =  2  Days
     Nitrate at Time = 4 Days
1.2 —
0.4 —
0.2-
     i    i     i    i     i    i     i    r
    0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4
                                               1.2 —
                                               0.4 —
                                               0.2-
                                                      /\
i     i    i     i    i     i    i     r
0.0 0.2  0.4 0.6  0.8 1.0  1.2 1.4
         Nutrient at Time = 2 Days
    Nutrient  at  Time =  4  Days
1.2 —
0.4 —
0.2-
     i    i     i    i     i    i     i    r
    0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4
                                               1.2 —
                                               0.4 —
                                               0.2-
i     i    i     i    i     i    i     r
0.0 0.2  0.4 0.6  0.8 1.0  1.2 1.4
                                      103

-------
104

-------
                             APPENDIX A: Data Input Guide



                   **** Data sets 2 through 26 must be preceded by a card *****

                    ******** containing the description of the data set ********
1.      PROBLEM IDENTIFICATION AND DESCRIPTION
       Two Records per Problem: FORMAT(I5,A72)

       Record one contains the following variables.
       (1) NPROB = Problem number
       (2) TITLE = Array for the title of the problem. It may contain up to 72 characters from column 6 to
              column 77.

       Record two contains the following 8 variables (FREE FORMAT)
       (1) IITR = Is non-linear convergence table information to be printed? 0= No, 1 = Yes.
       (2) ICHNG = Is cyclic change of rainfall-seepage nodes to be printed? 0 = No, 1 = Yes.
       (3) INTER = Is non-linear iterate information to be printed? 0 = No, 1 = Yes.
       (4) IBUG = Is debug information to be printed? 0 = No, 1 = Yes.
       (5) IMOD = Indicator of simulated problems:
              0  = Read input data only,
              1  = Only transport is simulated,
              10 = Only flow is simulated,
              11= Both flow and transport are simulated.
       (6) NITRHT = No. of iterations allowed for solving steady-state flow-transport iteration.
       (7) OMEHT = Iteration parameter of flow for coupled flow-transport iteration loop:
              0.0 ~ 1.0 = Under-relaxation,
              1.0 = Exact relaxation,
              1.0 ~ 2.0 = Over-relaxation.
       (8)  IGEOM = Integer indicating if (1) the geometry, boundary and pointer arrays are to be printed;
              (2) the boundary and pointer arrays are to be computed or read via logical units. If to be
              computed, they should be written on logical units. If IGEOM is even number, (1) will not be
              printed.  If IGEOM is odd number, (1) will be printed.  If IGEOM is less than or equal to 1,
              boundary arrays will be computed and written on logical unit, but if IGEOM is greater than
              1, boundary  arrays will be read via logical unit 13. If IGEOM is less than or equal to 3,
              pointer arrays will be will be computed  and written on logical unit  14, but if IGEOM is
              greater than 3, pointer arrays will be read via logical unit 14.
2.      OPTION PARAMETERS FOR BOTH FLOW AND TRANSPORT
       One record  contains the following 10 variables (FREE FORMAT)
       (1) IMID = Mid-difference control: 0 =No mid-difference, 1 = Mid-difference is used. (W=1.0)
       (2) IWET = Weighting factor control: 0 = Galerkin weighting, 1 = Upstream weighting.

                                            105

-------
       (3) ILUMP = Mass matrix lumping control: 0 = No lump, 1 = Lump.
       (4) IOPTIM = Optimization factor for computing indicator: 0 = Optimization factor is set to be 1.0,
              0.0, -1.0, or APHAG by the program depending on the velocity and dispersion coefficient,
              1 = Optimization factor is to be computed.
       (5) IPNTS = Is pointwise iteration method to be used to solve the matrix equations? 0 = No,
               1 =Yes.
       (6) IVML = Is velocity matrix to be lumped? 0 = No, 1 = Yes.
       (7) NSTRH = No. of logical records to be read via logical unit 11 for restarting calculation. 0 = No
              restart.
       (8) NSTRT = No. of logical records to be read via logical unit 12 for restarting calculation. 0 = No
              restart.
       (9) IQUAR = Index of using quadrature for numerical integration;
              11= Nodal quadrature for surface integration, Nodal quadrature for element integration,
              12 = Nodal quadrature for surface integration, Gaussian quadrature for element integration,
              21 = Gaussian quadrature for surface integration, Nodal quadrature for element integration.
              22 = Gaussian  quadrature for  surface  integration,  Gaussian quadrature  for  element
                      integration,
       (10) IDISP = Is dispersion/diffusion to be included in the transport simulation? 0 = No, 1 = yes.
3.      NODAL POINT COORDINATE
       This data set is required for both NSTRH and NSTRT equal to zero. Two subsets are included.
       Subset 1 contains one record telling the program total number of nodes in the region of interest.
       Subset 2 usually requires 2*NNP records, NNP record for the x-coordinate and NNP record for the
       z-coordinate. However, if a group of subsequent nodes appear in regular pattern, automatic generation
       can be made.

       A.  Subset 1 - FREE FORMAT: it contains one variable.
       (1) NNP = total number of global nodes in the region of interest.

       B.  Subset 2 - FREE FORMAT : Each record contains the following variables.
       (1) NI = Node number of the first node in the sequence.
       (2) NSEQ  = NSEQ subsequent nodes will be automatically generated.
       (3) NAD = Increment of node number for each of the NSEQ subsequent nodes.
       (4) XNI = x-coordinate of node NI, (L).
       (5) XAD = Increment of x-coordinate for each of the NSEQ subsequent nodes, (L).
       (6) XRD = Percentage of the increase of the increment over its preceding increment, (Decimal point);
               If XRD.EQ.O., all increments XAD's are the same. If XRD .GT. 0.0, the first increment is
               XAD*(1+XRD),  the second increment  is XAD*(1+XRD)**2, the third increment is
               XAD*(1+XRD)**3, and so on.

         *** NOTE: A line with 6 O's must be used to signal the end of this data set.

             A similar group of records is used to input z-coordinate *****
4.      ELEMENT INCIDENCE
       This data set is required for both NSTRH and NSTRT equal to zero. Two subsets are included.

                                             106

-------
Subset 1 contains three variables telling the program the total number of global elements and material
types in the region of interest, and the control index of material property.  Subset 2 usually requires
NEL records. However, if a group of elements appear in a regular pattern, automatic generation is
made.

A. Subset 1 - FREE FORMAT: it contains three variables.
(1) NEL = total number of elements in the region of interest.
(2) NMAT = No. of materials.
(3) KCP = Permeability input control:
       0 = Input hydraulic saturated conductivity,
       1 = Input intrinsic saturated permeability.

B. Subset 2 - FREE FORMAT: Each record contains the following variables.
(1) MI = Global element number.
(2) IE(MI,1) = Global node number of the first node of element MI.
(3) IE(MI,2) = Global node number of the second node of element MI.
(4) IE(MI,3) = Global node number of the third node of element MI.
(5) IE(MI,4) = Global node number of the forth node number MI.
(6) IE(MI,5) = Material type to be applied to element MI.
(7) MODL = Number of elements in the direction of most rapidly numbered nodes.
(8) NLAY = Number of elements in the direction of least rapidly numbered nodes.

IE(MI,1) IE(MI,4) are numbered beginning with the lower left corner and progressing around the
element in a counterclockwise direction. For a rectangular block of elements, it is only necessary to
specify the first element, the width MODL and the length NLAY, where MODL and NLAY are
measured in elements. The following figure provides an example.  The object is considered to be
rectangular since it has width MODL = 3 on two opposite sides and length NLAY = 5 on the other
two opposite sides. To generate definitions of elements 2 through 15 automatically, including both
the incidence and material type, only one card is necessary: Although all elements of this example will
be assumed to contain the same material type, MTYP =  1, this situation can easily be changed by
using material-correction facility.
                                 22
                                 21
                                      107

-------
        1   156235
5.      MATERIAL TYPE CORRECTION
       Two subsets of free-formatted data records are required for this data set.

       A. subset 1:
       (1) NCM = Number of elements with material corrections.

       B. subset 2: This set of data records is required only if NCM > 0.  Normally, NCM records are
              required. However, if a group of elements appear in a regular pattern, automatic generation
              may be made. Each record contains the following variables.
       (1) MI = Global element number of the first element in the sequence.
       (2) NSEQ = NSEQ subsequent elements will be generated automatically.
       (3) MAD = Increment of element number for each of the NSEQ subsequent elements.
       (4) MITYP = Type of material correction for element MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ subsequent elements.

       **** NOTE: A line with 5 O's must be used to signal the end of this data set.
6.      MATERIAL PROPERTIES
       Five subsets of free-formatted data records are required for this data set. There are 3*NMAT+2
       records in this data set.

       A. Subset 1 - FREE FORMAT: it contains three variables.
       (1) NMPPM = No. of parameters used for describing material properties. (It is set to be 10 originally).
       (2) NCC = No. of components in the system. Since the kinetic reaction model is built in the program
              according  to Eq. (2.9) through Eq. (2.15), NCC is assigned to 7 and IRXN is set to 1 if the
              microbial-chemical Monod type reactions are involved.  When NCC = 7 and the Monod
              kinetics are used, the first component must be Microbe # 1, the second component Microbe
              # 2, the third component Microbe #3, the fourth component the substrate, the fifth component
              Oxygen, the sixth component Nitrate, and the seventh component the nutrient.
              NCC can  be  1 for the single component simulation and NCC is equal to 2 by  using
              stochiometric model which results in IRXN = -1. NCC can be any value if users modify the
              kinetic model in the program (Subroutine ADVRX).
       (3)  IRXN  = the  index indicating the chemical-microbial kinetic reaction type.  -1 refers  to
              stochiometric reaction; 1 indicates Monod type reaction.

       **** Subsets 2, 3,  and 4 should be repeated NMAT times.
       B. Subset 2 - FREE FORMAT: it contains the following NMPPM variables for each material type.
       (1) PROP(I,1) = Saturated longitudinal hydraulic conductivity (intrinsic permeability) in X-direction
              of the I-th  material (L/T or L2).
       (2) PROP(I,2) = Saturated longitudinal hydraulic conductivity (intrinsic permeability) in Z-direction
              of the I-th  material (L/T or L2).
       (3)  PROP(I,3)  = Saturated transverse hydraulic conductivity (intrinsic permeability) of the I-th

                                             108

-------
              material (L/T or L2).
       (4) PROP(I,4) = Density of groundwater in the I-th material (M/L3).
       (5) PROP(I,5) = Gravity in the I-th material (L/T2).
       (6) PROP(I,6) = Longitudinal dispersivity in the I-th material (L).
       (7) PROP(I,7) = Transverse dispersivity in the I-th material (L).
       (8) PROP(I,8) = Diffusion coefficient in the I-th material (L2/T).
       (9) PROP(I,9) = Viscosity of groundwater in the I-th material (M/LT).
       (10) PROP(I,10) = Bulk density of the I-th material  (M/L3).
       (NMPPM) PROP(I,NMPPM) = The NMPPM-th material parameter of the I-th material.

       C. Subset 3 - FREE FORMAT: Each record contains NCC parameters (L3/M) as follows.
       (1) RKD(I,1) = Distribution coefficient of Microbe #1 in the I-th material.
       (NCC) RKD(I,NCC) = Distribution coefficient of the NCC-th component in the I-th material.

       D. Subset 4 - FREE FORMAT: each record contains NCC parameters (1/T) as follows.
       (1) TRANC(I,1) = Transformation rate constants of Microbe #1 in the I-th material.
       (NCC) TRANC(I,NCC) = Transformation rate constants of the NCC-th component in the  I-th
              material.

       E. Subset 5 - one record of FREE FORMAT: it contains NCC parameters (M/L3) as follows.
       (1) DINTS(l) = Intrinsic density of Microbe #1.
       (NCC) DINTS(NCC) = Intrinsic density of the NCC-th component.

       F. Subset 6 - one record of FREE FORMAT: it contains NCC parameters (L2/T) as follows.
       (1) RHOMU(l) = Viscosity effecting factor of Microbe #1.
       (NCC) RHOMU(NCC) = Viscosity effecting factor of the NCC-th component.


7.      SOIL PROPERTIES
       Three or five subsets of free-formatted data records are required for this data set depending on the
       given forms of the soil property functions.

       A. Subset 1: Soil property control parameters

                                            109

-------
(1) KSP = Soil property input control, 0 = analytical input, 1 = Tabular data input.
(2) NSPPM = Number of points in tabular soil property functions or number of parameters to specify
       analytical soil functions per material.

B. Subset 2a: Analytical soil parameters - This sub-data-set is needed if and only if KSP is 0. Two
       sets of records are required, one  for  moisture-content parameters and  the other for
       conductivity (permeability) parameters.
(1) SPP(J,I,1) = Analytical moisture-content parameter J of material I, J = 1..NSPPM. NMAT sets
       of these parameters are required for I = 1..NMAT. That is, if SPP(J,I,l)for J=  1..NSPPM can
       be put on a single line, we need NMAT consecutive line for the sets of parameters.
(2) SPP(J,I,2) = Analytical relative conductivity parameter J of material I.  Similar input data setting
       is required for these parameters as for SPP(J,I,1)

C. Subset 2b: Soil properties in tabular form - This sub-data-set is needed if and only if KSP not is 0.
       Four sets of records are needed — for pressure, water-content, relative conductivity (or relative
       permeability), and water capacity, respectively.
(1) SPP(J,I,4) = Tabular value of pressure head of the J-th point for material I. NMAT sets of these
       parameters are required for I = 1..NMAT. That is, if SPP(J,I,4) for J = 1..NSPPM can be put
       on a single line, we need NMAT consecutive line for the sets of parameters.
(2) SPP(J,I,1) = Tabular value of moisture-content of the J-th point in material I. Similar input data
       setting is required for these  parameters as for SPP(J,I,4)
(3) SPP(J,I,2) = Tabular value of relative conductivity of the J-th point in material I. Similar input
       data setting is required for these parameters as for SPP(J,I,4)
(4) SPP(J,I,3) = Tabular value of moisture-content capacity of the J-th point in material I.  Similar
       input data setting is required for these parameters as for SPP(J,I,4)
TIME CONTROL PARAMETERS
Five subsets of data records are required for this data set.

A. Subset 1: free format
(1) NTI = Number of time steps or time increments.
(2) NDTCHG = No. of times to reset time step size to initial time step size.

B. Subset 2: free format
(1) BELT = Initial time step size, (T).
(2) CHNG = Percentage of change in time step size in each of the subsequent time increment,
       (dimensionless in decimal point).
(3) DELMAX = Maximum value of DELT, (T).
(4) TMAX = Maximum simulation time, (T).

C. Subset 3: format = 8011
(1) KPRO = Printer control for steady state and initial conditions;
       0 = print nothing,
       1 = print FLOW, FRATE, and TFLOW,
       2 = print above (1) plus pressure head H,
       3 = print above (2) plus total head,
       4 = print above (3) plus moisture content,

                                       110

-------
              5 = print above (4) plus Darcy velocity.
       (2) KPR(I) = Printer control for the I-th (I = 1, 2, ..., NTI) time step similar to KPRO.

       D. Subset 4: format = 8011
       (1) KDSKO = Auxiliary storage control for steady state and initial condition;
              0 = no storage, 1 = store on Logical Unit 1.
       (2) KDSK(I) = Auxiliary storage control for the I-th (I = 1, 2, ..., NTI) time step similar to KDSKO.

       E. Subset 5: free format
       (1) TDTCH(I) = Time when the I-th (I = 1, 2,..., NDTCHG) step-size-resetting is needed.

       * * * * NOTE: NTI can be computed by NTI = 11 + 1+12 + 1, where 11= largest integer not exceeding
       Log(DELMAX/DELT)/Log(l+CHNG),       12   =    largest   integer   not   exceeding
       (RTIME-DELT*((1+CHNG)**(I1+1)-1)/CHNG)/DELMAX,  RTIME = Real  simulation time,
       DELMAX,DELT,and CHNG are defined in Data Set 3.
The following data sets, data set 9 to data set 16, must be skipped as IMOD equals 1.

9.      BASIC REAL AND INTEGER PARAMETERS FOR FLOW
       Two subsets of free-formatted data records are required for this data set.

       A. Subset 1 - one record of basic integer for flow part (FREE FORMAT).
       (1) NCYL = No. of cycles allowed for iterating variable boundary conditions.
       (2) NITERH = No. of iterations allowed for solving nonlinear equations.
       (3) NPITERH = No.  of iterations allowed for pointwise iteration.
       (4) KSTRH = Auxiliary storage output control: 0 = No storage, 1 = Output stored on logical unit  11.
       (5) KSSH = Steady  state control for flow: 0 = Steady state solution is desired, 1 = Only transient
              solutions are desired.
       (6) KGRAV = Gravity term control: 0 = No gravity term (Used for horizontal flow only),
              1 = With gravity term.

       B. Subset 2 - one record of basic real number for flow part (FREE FORMAT).
       (1) TOLAH = Steady-state convergence criteria,
       (2) TOLBH = Transient-state convergence criteria.
       (3) WF = Time derivative weighting factor: 0.5 = Crank-Nicolson central, 1.0 = Backward difference
              and/or mid-difference.
       (4) OMEH = Iteration parameter for solving the nonlinear matrix equation:
              0.0 ~ 1.0 = Under-relaxation,
              1.0 = Exact relaxation,
              1.0 ~ 2.0 = Over-relaxation.
       (5) OMIH = Relaxation parameter for solving the matrix equation pointwisely:
              0.0 ~ 1.0 = Under-relaxation,
              1.0 = Exact relaxation,
              1.0 ~ 2.0 = Over-relaxation.
       (6) CNSTKR = Constraint on relative hydraulic conductivity;
              0 = no constraint,

                                             111

-------
              0.0001, 0.001, or 0.01 should be tried when nonconvergency occurs in solving the nonlinear
              flow equation.
10.  CARD INPUT FOR INITIAL OR PRE-INITIAL CONDITIONS OF FLOW PART
       Free-formatted data records are required for this data set.  Generally, NNP records are needed.
       However, if a group of nodes appear in a regular pattern, auto-generation is made.

       Initial pressure head - each record contains the following variables.
       (1) NI = Global node number of the first node in the sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) NAD = Increment of node number for each of the NSEQ nodes.
       (4) HNI = Initial or pre-initial pressure head of node NI, (L).
       (5) HAD = Increment of initial or pre-initial head for each of the NSEQ nodes, (L).
       (6) HRD = 0.0

       **** NOTE: A line with 6 O's must be used to signal the end of this data set.

       NOTE ON INITIAL CONDITIONS AND RESTARTING: The initial condition for a transient
       calculation may be obtained in two different ways: from card  input, or steady-state calculation using
       a time-invariant boundary condition that is different from those for transient computation.  In the latter
       case a card input of the pre-initial conditions is required  as the zero-th oder iterate of the steady state
       solution.

       NOTE ON STEADY STATE INPUT:  Steady state option may be used to provide either the final
       state of a system under study or the initial conditions for a transient state calculation. In former case
       KSSh = 0, KSSt > 0, and NTI = 0, and in the latter case KSSh = 0 or KSSt > 0 and NTI > 0.  If KSSh
       > 0, and KSSt > 0, there will be no steady state calculation.
11.  ELEMENT (DISTRIBUTED) SOURCE/SINK FOR FLOW SIMULATIONS
       Four subsets of free-formatted data records are required in this data set.

       A. Subset 1: control parameters
       (1) NSELH= No. of source/sink elements.
       (2) NSPRH= No. of source/sink profiles.
       (3) NSDPH= No. of data points in each of the NSPR source/sink profiles.
       (4) KSAIH= Is element-source/sink profile to be input analytically, 0 = no, 1 = yes.

       B. Subset 2: source/sink profiles - This group of data is needed if and only if NSELH .GT. 0. For
              each sub-data-record, NSDPH of the data pair (TSOSFH(J,I),SOSFH(J,I)) are required. If this
              sub-data-record can be fitted in a line, we will need NSPRH lines.
       (1)  TSOSFH(J,I) = Time of the J-th data point in the I-th profile, (T).
       (2)  SOSFH(1,I)  = Source/sink value of the J-th data point in the I-th profile, (L**3/T/L**2/L).

       C. Subset 3: global source/sink element number - This group of data is needed if and only if NSELH
              .GT. 0. NSELH data points are required for this record.
       (1)  MI = The first compressed element in the sequence.

                                             112

-------
       (2) NSEQ = NSEQ elements will be generated automatically.
       (3) MAD = Increment of element number for each of the NSEQ elements.
       (4) MITYP = Global element number of element MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ elements.

       **** NOTE: A line with 5 O's is used to signal the end of this data set.

       D. Subset 4: Source type assigned to each element - Usually one record per element.  However,
              automatic generation can be made. For I-th (I = 1, 2, ...,) record, it contains the following.
       (1) MI = Global element number of the first element in the sequence.
       (2) NSEQ = NSEQ elements will be generated automatically.
       (3) MAD = Increment of element number for each of the NSEQ elements.
       (4) MITYP = Source type in element MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ elements.

       **** NOTE: A line with 5 O's is used to signal the end of this data set.
12.  POINT (WELL) SOURCE/SINK DATA FOR FLOW SIMULATION
       Four subsets of free-formatted data records are required for this data set.

       A. Subset 1: control parameters
       (1) NWNPH = No. of well or point source/sink nodal points.
       (2) NWPRH = No. of well or point source/sink strength profiles.
       (3) NWDPH = No. of data points in each of the NWPR profiles.
       (4) KWAIH = Is well-source/sink profile to be input analytically, 0 = no, 1 = yes.

       B. Subset 2: source/sink profiles - This group of data is needed if and only if NWNPH .GT. 0. For
              each sub-data-record, NWDPH of the data pair (TWSSFH(J,I),WSSFH(J,I)) are required.
              If this sub-data-record can be fitted in a line, we will need NWPRH lines.
       (1) TWSSFH(J,I) = Time of the J-th data point in the I-th profile, (T).
       (2).  WSSFH(J,I) = Source/sink  value of the J-th data point in the I-th profile, (L**3/T/L).

       C. Subset 3: global source/sink node number - This group of data is needed if and only if NWNPH
              .GT. 0. NWNPH data points are required for this record.
       (1) NI = Compressed well node number of the first node in the sequence.
       (2) NSEQ = NSEQ nodes will be generated automatically.
       (3) NAD = Increment of well node number for each of the NSEQ nodes.
       (4) NITYP = Global node number of node NI.
       (5) NTYPAD = Increment of NITYP for each of the NSEQ nodes.

       **** NOTE: A line with 5 O's is used to signal the end of this data set.

       D. Subset 4:  Source type assigned to each well - Usually one record per element.  However, automatic
              generation can be made. For I-th (I = 1, 2, ...,) record, it contains the following.
       (1) NI = Compressed well node number of the first node in the sequence.
       (2) NSEQ = NSEQ nodes will be generated automatically.

                                            113

-------
       (3) NAD = Increment of well node number for each of the NSEQ nodes.
       (4) NITYP = Source type in node NI.
       (5) NTYPAD = Increment of NITYP for each of the NSEQ nodes.

       **** NOTE:  A line with 5 O's is used to signal the end of this data set.
13. RAINFALL/EVAPORATION-SEEPAGE BOUNDARY CONDITIONS
       Seven subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NRES = No. of variable boundary element sides.
       (2) NRNP = No. of variable boundary nodal points.
       (3) NRPR = No. of rainfall profiles,
       (4) NRDP = No. of rainfall data points in each of the NRPR rainfall profiles.
       (5) KRAIH = Is rainfall profile to be input analytically, 0 = no, 1 = yes.

       B. Subset 2: boundary profiles - This subset is required only when NRES is not 0. NRPR profiles are
              needed. For each profile, NRDP of the data pair (TRF(J,I),RF(J,I)) are required.  If these data
              pairs can fit in a line, we will need NRPR of data lines.
       (1) TRF(J,I) = Time of the J-th data point in the I-th profile, (T).
       (2) RF(J,I) = Rainfall/evaporation rate of the J-th data point in the I-th profile, (L/T).

       C. Subset 3: Global Node Number of All Compressed Variable Boundary (VB) Nodes.  At most,
              NRNP records are needed for this subdata set, one each for NRNP variable boundary nodes.
              For I-th (I = 1, 2, ...,) Record, it contains the following 5 variables.
       (1) NI = Compressed VB node number of the first node in the sequence.
       (2) NSEQ = NSEQ nodes will be generated automatically.
       (3) NIAD = Increment of NI for each of the NSEQ nodes.
       (4) NODE = Global node number of node NI,
       (5) NODEAD = Increment of NODE for each of the NSEQ nodes.

       **** NOTE: A line with 5 O's is used to signal the end of this data set.

       D. Subset 4:  boundary profile types assigned to each VB point. At most NRNP records are needed.
              However, automatic generation can be made. For I-th (1=1,2,...,) record, it  contains the
              following variables.
       (1) MI = Compressed VB node of the first node in the sequence.
       (2) NSEQ = NSEQ nodes will be generated automatically.
       (3) MIAD = Increment of NI for each of the NSEQ nodes.
       (4) MITYP = Type of rainfall/evaporation profiles assigned to node MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ nodes.

       **** NOTE: A line with 5 O's is used to signal the end of this data set.
       E. Subset 5: Ponding Depth Allowed in Each of NRNP Variable Boundary Nodes. Normally, NRNP
              records are needed, one for each of the NRNP nodes.  However, if a group of nodes has a

                                            114

-------
              regular pattern of ponding depth, automatic generation is made.  For I-th (1=1,2,...,) record,
              it contains the following variables.

       (1) NI = Compressed VB node number of the first node in a sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) NIAD = Increment of NI for each of the NSEQ subsequent nodes.
       (4) HCONNI = Ponding depth of node NI, (L).
       (5) HCONAD = Increment of HCONNI for each of the NSEQ nodes, (L).
       (6) HCONRD = 0.0

       **** NOTE: A line with 6 O's must be used to signal the end of this subdata set.

       F. Subset 6: Minimum  Pressure  Head Allowed in Each NRNP Variable Boundary Nodes. This
              subdata set is read-in similar to the above subdata set. For I-th (I = 1, 2, ..., ) record, it
              contains the following variables.
       (1) NI = Compressed VB node number of the first node in a sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) NIAD = Increment of NI for each of the NSEQ subsequent nodes.
       (4) HMINNI = Minimum pressure head allow for node NI, (L).
       (5) HMINAD  = Increment of HMINNI for each of the NSEQ nodes, (L).
       (6) HMINRD  = 0.0

       **** NOTE: A line with 6 O's must be used to signal the end of this subdata set.
       G. Subset 7:  Specification of Rainfall/evaporation-seepage Sides.  Normally, NRES records are
              required, one each for a variable boundary (VB) element side.  However, if a group of
              rainfall/evaporation-seepage element sides appears in a regular pattern, automatic generation
              may be made. For I-th (I = 1, 2,...,) record, it contains the following variables.
       (1) MI = Compressed VB element side number of the first element side in a sequence.
       (2) NSEQ = NSEQ subsequent VB element sides will be generated automatically.
       (3) M = Global element no. to which the Ml-th RS  side belongs.
       (4) IS1 = Compressed RS node number of the first  node of element side MI.
       (5) IS2 = Compressed RS node no. of the 2-nd node of the Ml-th RS side.
       (6) MIAD = Increment of MI for each of the NSEQ subsequent VB element sides.
       (7) MAD = Increment of M for each of the NSEQ subsequent RS side.
       (8) IS1AD = Increment of IS1 for each of the NSEQ subsequent RS side.
       (9) IS2AD = Increment of IS2 for each of the NSEQ subsequent RS side.

       **** NOTE:  A line with 9 O's must be used to signal the end of this subdata set.
14. DIRICHLET BOUNDARY CONDITIONS FOR FLOW SIMULATION
       Four subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NDNPH = No. of Dirichlet nodal points, should be .GE. 1.
       (2) NDPRH = No. of total Dirichlet-head profiles, should be .GE. 1.
                                            115

-------
       (3) NDDPH = No. of data points in each total head profiles, should be .GE. 1.
       (4) KDAIH = Is Dirichlet boundary value profile to be input analytically, 0 = no, 1 = yes

       B. Subset 2: Dirichlet-head profiles - This sub-data-set is required only if NDNPH is not 0.  NDPRH
              of profiles are needed.  For each profile, NDDPH of the data pair (THDBF(J,I),HDBF(J,I))
              are needed. If these data pairs can fit in a line, we will need NDPRH lines.
       (1) THDBF(J,I) = Time of the J-th data point in the I-th profile, (T).
       (2) HDBF(J,I) = Total head of the J-th data point in the I-th profile, (L).

       C. Subset 3: Global Node Number of All Compressed Dirichlet Boundary Nodes.  At most, NDNPH
              records are needed for this subdata set,  one each for NDNPH Dirichlet boundary nodes. For
              I-th (I = 1,2,  ...,) record, it contains the following 5 variables.
       (1) NI = Compressed Dirichlet node number of the first node in the sequence.
       (2) NSEQ = NSEQ subsequent Dirichlet nodes will be generated automatically.
       (3) NAD = Increment of NI for each of the NSEQ nodes.
       (4) NITYP = Global node number for node NI and NSEQ subsequent nodes.
       (5) NTYPAD = Increment of NITYP for each  of the NSEQ subsequent nodes.

       **** NOTE: A line with 5 O's must be used to signal the end of this  sub-data set.

       D. Subset 4: boundary profile type assign to each Dirichlet node - Normally one record per Dirichlet
              node, i. e. a total of NDNPH records. However, if the Dirichlet nodes appear in regular
              pattern, automatic generation may be made. For I-th (1=1,2,...,) record, it contains the
              following variables.
       (1) NI = Compressed Dirichlet node number of the first node in the sequence.
       (2) NSEQ = NSEQ subsequent Dirichlet nodes will be generated automatically.
       (3) NAD = Increment of NI for each of the NSEQ nodes.
       (4) NITYP = Type of total head profile for node NI and NSEQ subsequent nodes.
       (5) NTYPAD = Increment of NITYP for each  of the NSEQ subsequent nodes.

       **** NOTE: A line with 5 O's must be used to signal the end of this  sub-data set.
15.  CAUCHY BOUNDARY CONDITIONS FOR FLOW SIMULATIONS
       Five subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NCESH = No. of Cauchy boundary element sides.
       (2) NCNPH = No.  of Cauchy nodal points.
       (3) NCPRH = No. of Cauchy-flux profiles.
       (4) NCDPH = No.  of data points in each of the NCPR Cauchy-flux profiles.
       (5) KCAIH = Is Cauchy flux profile to be input analytically, 0 = no, 1 = yes

       B. Subset 2: prescribed Cauchy-flux profiles - This sub-data-set is required only if NCESH is not 0.
              NCPRH of profiles  are  needed.   For  each profile,  NCDPH  of the data  pair
              (TQCBF(J,I),QCBF(J,I)) are needed.  If these data pairs can fit in a line, we will need
              NDPRH lines.
       (1) TQCBF(J,I) = Time of the J-th data point in the I-th profile, (T).

                                            116

-------
       (2) QCBF(J,I) = Normal Cauchy flux of the J-th data point in the I-th profile, (L**3/T/L**2); positive
              out from the region, negative into the region.

       C.  Subset 3: global node number of all compressed Cauchy nodes - At most, NCNPH records are
              needed for this subdata set, one each for NCNPH Cauchy nodes.
       (1) NI = Compressed Cauchy node number of the first node in the sequence.
       (2) NSEQ = NSEQ nodes will be generated automatically.
       (3) NIAD = Increment of NI for each of the NSEQ sides.
       (4) NITYP = Global node number of node NI.
       (5) NTYPAD = Increment of NITYP for each of the NSEQ sides.

       **** NOTE:  A line with 5 O's is used to signal the end of this data set.

       D.  Subset 4: type of Cauchy flux profiles assigned to each of all NCNPH sides. At most NCNPH
              records are needed.  However, automatic generation can be made.  For I-th (I = 1, 2, ...,)
              record, it contains the following variables.
       (1) MI = Compressed Cauchy node number of the first node in the sequence.
       (2) NSEQ = NSEQ nodes will be generated automatically.
       (3) MIAD = Increment of MI for each of the NSEQ nodes.
       (4) MITYP = Type of Cauchy flux profile assigned to node MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ nodes.

       **** NOTE:  A line with 5 O's is used to signal the end of this data set.

       E. Subset 4: Cauchy boundary element  sides - Normally, NCESH records are required, one each for
              a  Cauchy boundary element side.  However, if a group of Cauchy boundary element sides
              appears in a regular pattern, automatic generation may be made. For I-th (I = 1, 2, ..., )
              record, it contains the following variables.
       (1) MI = Compressed Cauchy element  side number of the first element-side in a sequence.
       (2) NSEQ = NSEQ subsequent Cauchy element-sides will be generated automatically.
       (3) M = Global element no. to which the Ml-th Cauchy side belongs.
       (4) IS1 = Compressed Cauchy node number of the first node on the Cauchy element-side MI.
       (5) IS2 = Compressed Cauchy node number of the second node on the Cauchy element-side MI.
       (6) MIAD = Increment of MI for each of the NSEQ subsequent sides.
       (7) MAD = Increment of M for each of the NSEQ subsequent Cauchy sides.
       (8) IS IAD = Increment of IS1 for each of the NSEQ subsequent element-sides.
       (9) IS2AD = Increment of IS2 for each of the NSEQ subsequent element-sides.

       **** NOTE:  A line with 9 O's is used to end this data set input.
16.  NEUMANN BOUNDARY CONDITIONS FOR FLOW SIMULATIONS
       Five subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NNESH = No. of Neumann boundary element sides.
       (2) NNNPH = No. of Neuamnn nodal points.
       (3) NNPRH = No. of Neumann flux profiles.
                                           117

-------
(4) NNDPH = No. of data points in each of the NNPR Neumann-flux profiles.
(5) KNAIH = Is Neumann flux profile to be input analytically, 0 = no, 1 = yes

**** NOTE: A line with 5 O's is used to signal the end of this data set.

B. Subset 4: prescribed Neumann-flux profiles - This sub-data-set is required only if NNESH is not
       0.   NNPRH  of profiles are needed.   For each  profile,  NNDPH of the data pair
       (TQNBF(J,I),QNBF(J,I)) are needed.  If these data pairs can fit in a line, we will need
       NDPRH lines.
(1) TQNBF(J,I) = Time of the J-th data point in the I-th profile, (T).
(2)  QNBF(J,I) =  Normal Neumann flux of the J-th data point in the I-th profile, (L**3/T/L**2);
       positive out from the region, negative into the region.

C. Subset 3: global node number of all compressed Neumann nodes -At most, NNNPH records are
       needed for this sub-data set, one each for NNNPH Neumann nodes.
(1) NI = Compressed Neumann node number of the first node in the sequence.
(2) NSEQ  = NSEQ nodes will be generated automatically.
(3) NIAD = Increment of NI for each of the NSEQ sides.
(4) NITYP = Global node number of node NI.
(5) NTYPAD = Increment of NITYP for each of the NSEQ sides.

**** NOTE: A line with 5 O's is used to signal the end of this data set.

D. Subset 4: type of Neumann flux profiles assigned to each of all NNNPH nodes.  At most NNNPH
       records are needed.  However, automatic generation can be made. For I-th (I = 1, 2, ...,)
       record, it contains the following variables.
(1) MI = Compressed Neumann node number of the first node in the sequence.
(2) NSEQ  = NSEQ nodes will be generated automatically.
(3) MIAD = Increment of MI for each of the NSEQ nodes.
(4) MITYP = Type of Neumann flux profile assigned to node MI.
(5) MTYPAD = Increment of MITYP for each of the NSEQ nodes.

**** NOTE: A line with 5 O's is used to signal the end of this data set.

E. Subset 5: Neumann boundary element sides - Normally, NNESH records are required, one each
       for a Neumann boundary element side.  However, if a group of Neumann boundary element
       sides appears in a regular pattern, automatic generation may be made. For I-th (I = 1,2, ...,
       ) record, it contains the following variables.
(1) MI = Compressed Neumann side number of the first side in sequence.
(2) NSEQ  = NSEQ subsequent Neumann sides will be generated automatically.
(3) M = Global element no. to which the Ml-th Neumann side belongs.
(4) IS1 = Compressed Neumann node number of the first node  on the Neumann element-side MI. (5)
IS2 = Compressed Neumann node number of the second node on Neumann element-side MI.
(6) MIAD = Increment of MI for each of the NSEQ subsequent sides.
(7) MAD = Increment of M for each of the NSEQ subsequent Neumann side.
(8) IS 1 AD = Increment of IS 1 for each of the NSEQ subsequent Neumann side.
(9) IS2AD = Increment of IS2 for each of the NSEQ subsequent Neumann side.
**** NOTE: A line with 9 O's is used to end this data set input.

                                     118

-------
The following data sets, data set 17 to data set 26, must be skipped as IMOD equals 10.

17. BASIC REAL AND INTEGER PARAMETERS FOR TRANSPORT
       Four subsets of free-formatted data records are required for this data set.

       A. Subset 1: This group of data reads the integer control parameter for transport. Five variables are
               included.
       (1) NITERT = No. of iterations allowed for solving nonlinear equations.
       (2) NPITERT = No. of iterations allowed for pointwise iteration.
       (3) KSTRT = Auxiliary storage output control: 0 = No storage, 1 = Output stored on logical unit 11.
       (4) KSST = Steady state control for flow: 0 = Steady state solution is desired,
               1 = Only transient solutions are desired.
       (5) KVIT = Velocity input control:
               -1 = Card input for velocity and moisture content.
               1 = Steady-state velocity and moisture content input via logical unit 12.
               2 = transient velocity and moisture content input via auxiliary storage logical unit 12.
       (6) MICONF = Index of the simulation of microbial configuation
               0 = mobile microbes
               1 = immobile microbes

       B. Subset 2: This group of data has seven variables.
       (1) TOLAT = Steady-state convergence criteria.
       (2) TOLBT = Transient-state convergence criteria.
       (3) WT = Time derivative weighting factor:
               0.5 = Crank-Nicolson central,
               1.0 = Backward difference and/or mid-difference.
       (4) WV = Time integration factor for the velocity terms:
               0.0 = Forward difference  (explicit),
               0.5 = Central difference,
               1.0 = Backward difference (implicit).
       (5) OMET = Iteration parameter for solving the nonlinear matrix equation:
               0.0 ~ 1.0 = Under-relaxation,
               1.0 = Exact relaxation,
               1.0 ~ 2.0 = Over-relaxation.
       (6) OMIT = Relaxation parameter for solving the matrix equation pointwisely:
               0.0 ~ 1.0 = Under-relaxation,
               1.0 = Exact relaxation,
               1.0 ~ 2.0 = Over-relaxation.
       (7) APHAG = Upstream weighting factor used if IOPTIM.EQ.O.
               NOTE: (1) The value is between 0.0 and 1.5,
                      (2) If APHAG.GE.1.34, the program will choose  appropriate values of weighting
                             factor.

       C. Subset 3: It contains the following 10 variables  (FREE FORMAT)
       (1) IZOOM = Is zooming needed  for advection computation? 0 = No, 1 = Yes.

                                              119

-------
       (2) IDZOOM = Is zooming needed for dispersion computation? 0 = No, 1 = Yes.
       (3) IEPC = Is EPCOF scheme included? 0 = No, 1 = Yes.
       (4) NXA = No. of regularly refined grids for the advection step in the X-drection in an element.
       (5) NZA = No. of regularly refined grids for the advection step in the Z-drection in an element.
       (6) NXD = No. of dispersion fine grids in the X-direction.
       (7) NZD = No. of dispersion fine grids in the Z-direction.
       (8) NXW = No. of sub-grid points for the sub-element tracking in the X-direction in an element.
       (9) NZW = No. of sub-grid points for the sub-element tracking in the Z-direction in an element.
       (10)  IDETQ = Index of particle tracking patten;
               1 = Average velocity is used (more accurate);
              2 = Single velocity of the starting point is used (less computation).

       D. Subset 4: It reads the following 2 variables (FREE FORMAT)
       (1) ADPEPS = Error tolerance of relative concentration and nonlinear convergence criteria.
       (2) ADPARM = Error tolerance of concentration relative to maximum concentration.
18.  INITIAL/PRE-INITIAL CONDITIONS FOR TRANSPORT
       There are seven subsets in this data set:

       A. Subset 1: This group of data reads the initial conditions (M/L3) of transport for Microbe #1. There
              are 6 variables in each line of this subset (FREE FORMAT)
       (1) NI = The global node number of the 1-st node of the line.
       (2) NSEQ = No. of subsequent nodes to be counted in the line.
       (3) NAD = The increment of global node no. for each subsequent node in the line.
       (4) FNI = The initial/pre-initial total concentration of the 1-st node of the line.
       (5) FAD = The increment of total concentration for each subsequent node in the line.
       (6) FRD = The increment ratio of total concentration for each subsequent node in the line, (i.e., the
              increment increases with rate 1.0+FRD). Usually, it is set to be 0 for transport.

       NOTE: All the values of variables in the last line of the subset should be 0.

       B. Subset 2 ~ Subset NCC: The 2-nd to NCC-th subsets reads the initial conditions of transport for
              Microbe  #2,  Microbe #3, substrate, Oxygen, Nitrate, nutrient, and up to the NCC-th
              component.
19.  MICROBE-CHEMICAL INTERACTION CONSTANTS
       Eleven subsets of FREE-FORMATTED data are needed.

       A. Subset 1: Four parameters describing the specific growth rate of microbes (1/T) are needed in this
              record. If there are no microbes in this system, the following four numbers have to be set to
              zeros.
       (1) GRATE(l) = Maximum specific oxygen-based growth rate for microbe #1. (|4(1) in Eqs. (2.9) ~
              (2.15)).
       (2) GRATE(2) = Maximum specific nitrate-based growth rate for microbe #2. (^ in Eqs. (2.9) ~
              (2.15)).
       (3) GRATE(3) = Maximum specific oxygen-based growth rate for microbe #3. (|4(3) in Eqs. (2.9) ~

                                             120

-------
        (2.15)).
(4) GRATE(4) = Maximum specific nitrate-based growth rate for microbe #3. (^ in Eqs. (2.9) ~
        (2.15)).

B. Subset 2 : Four yield coefficient (M/M) are needed in this record and these four values can not be
        zeros.
(1) YCOEFF(l) = Yield coefficient for microbe #1 utilizing Oxygen. (Y0(1) in Eqs. (2.9) & (2.12)).
(2) YCOEFF(2) = Yield coefficient for microbe #2 utilizing Nitrate. (Yn(2) in Eqs. (2.9) & (2.12)).
(3) YCOEFF(3) = Yield coefficient for microbe #3 utilizing Oxygen. (Y0(3) in Eqs. (2.9) & (2.12)).
(4) YCOEFF(4) = Yield coefficient for microbe #3 utilizing Nitrate. (Yn(3) in Eqs. (2.9) & (2.12)).

C. Subset 3: Four retarded substrate saturation constants (M/L3) are needed in this record.
(1) RTARDS(l) = Retarded substrate saturation constant under aerobic conditions with respect to
        microbe #1.  (Kj1' in Eqs. (2.9)  ~ (2.15)).
(2) RTARDS(2) = Retarded substrate saturation constant under anaerobic conditions with respect to
        microbe #2.  (Kj-2* in Eqs. (2.9)  ~ (2.15)).
(3) RTARDS(3) = Retarded substrate saturation constant under aerobic conditions with respect to
        microbe #3.  (Kj3' in Eqs. (2.9)  ~ (2.15)).
(4) RTARDS(4) = Retarded substrate saturation constant under anaerobic conditions with respect to
        microbe #3.  (Kj3' in Eqs. (2.9)  ~ (2.15)).

D. Subset 4: Four retarded Oxygen or Nitrate saturation constant (M/L3) are needed in this record.
(1) RTARDO(l) = Retarded Oxygen saturation constant under aerobic conditions with respect to
        microbe #1.  (K0(1) in Eqs. (2.9) ~ (2.15)).
(2) RTARDO(2) = Retarded Nitrate saturation constant under anaerobic conditions with respect to
        microbe #2.  (K™ in Eqs. (2.9) ~ (2.15)).
(3) RTARDO(3) = Retarded Oxygen saturation constant under aerobic conditions with respect to
        microbe #3.  (K0(3) in Eqs. (2.9) ~ (2.15)).
(4) RTARDO(4) = Retarded Nitrate saturation constant under anaerobic conditions with respect to
        microbe #3.  (K^ in Eqs. (2.9) ~ (2.15)).

E. Subset 5: Four retarded nutrient saturation constant (M/L3) are needed in this record.
(1) RTARDN(l) = Retarded nutrient saturation constant under aerobic conditions with respect to
        microbe #1.  (Kp0(1) in Eqs. (2.9)  ~ (2.15)).
(2) RTARDN(2) = Retarded nutrient saturation constant under anaerobic conditions with respect to
        microbe #2.  (Kpn(2) in Eqs. (2.9)  ~ (2.15)).
(3) RTARDN(3) = Retarded nutrient saturation constant under aerobic conditions with respect to
        microbe #3.  (Kp0(3) in Eqs. (2.9)  ~ (2.15)).
(4) RTARDN(4) = Retarded nutrient saturation constant under anaerobic conditions with respect to
        microbe #3.  (Kpn(3) in Eqs. (2.9)  ~ (2.15)).

F. Subset 6: Four Oxygen-use of Nitrate-use coefficient for synthesis are needed in this record.
(1) SCOEFF(l) = Oxygen-use coefficient for synthesis by microbe #1.  (y0(1) m Eq. (2.10)).
(2) SCOEFF(2) = Nitrate-use coefficient for synthesis by microbe #2. (yn(2) in Eq. (2.11)).
(3) SCOEFF(3) = Oxygen-use coefficient for synthesis by microbe #3.  (y0(3) in Eq. (2.10)).
(4) SCOEFF(4) = Nitrate-use coefficient for synthesis by microbe #3. (yn(3) in Eq. (2.11)).

G. Subset 7: Four Oxygen-use of Nitrate-use coefficient for energy are needed in this record.

                                        121

-------
       (1) ECOEFF(l) = Oxygen-use coefficient for energy by microbe #1.  (ce0(1) in Eq. (2.10)).
       (2) ECOEFF(2) = Nitrate-use coefficient for energy by microbe #2.  (cen(2) in Eq. (2.11)).
       (3) ECOEFF(3) = Oxygen-use coefficient for energy by microbe #3.  (ce0(3) in Eq. (2.10)).
       (4) ECOEFF(4) = Nitrate-use coefficient for energy by microbe #3.  (cen(3) in Eq. (2.11)).

       H. Subset 8: Four microbial decay coefficient (1/T) are needed in this record.
       (1) DCOEFF(l) = Microbial decay coefficient of aerobic respiration of microbe #1.  (A0(1) in Eqs.
               (2.13) & (2.15)).
       (2) DCOEFF(2) = Microbial decay coefficient of anaerobic respiration of microbe #2.  (An(2) in Eqs.
               (2.14) & (2.15)).
       (3) DCOEFF(3) = Microbial decay coefficient of aerobic respiration of microbe #3.  (A0(3) in Eqs.
               (2.13) & (2.15)).
       (4) DCOEFF(4) = Microbial decay coefficient of anaerobic respiration of microbe #3.  (An(3) in Eqs.
               (2.14) & (2.15)).

       I. Subset 9: Four Oxygen or Nitrate saturation constants (M/L3) for decay are needed in this record.
       (1) SATURC(l) = Oxygen-saturation constant for decay with respect to microbe #1. (F0(1) in Eq.
               (2.10)).
       (2) SATURC(2) = Nitrate-saturation constant for decay with respect to microbe #2.  (Fn(2) in Eq.
               (2.H)).
       (3) SATURC(3) = Oxygen-saturation constant for decay with respect to microbe #3. (F0(3) in Eq.
               (2.10)).
       (4) SATURC(4) = Nitrate-saturation constant for decay with respect to microbe #3.  (Fn(3) in Eq.
               (2.H)).

       J. Subset 10: Four nutrient-use coefficient for the production are needed in this record.
       (1) PCOEFF(l) = Nutrient-use coefficient for the production of microbe #1 with aerobic respiration.
               (£„<'> in Eq. (2.12)).
       (2) PCOEFF(2)  = Nutrient-use coefficient for the production of microbe  #2 with anaerobic
               respiration. (en(2) in Eq. (2.12)).
       (3) PCOEFF(3) = Nutrient-use coefficient for the production of microbe #3 with aerobic respiration.
               (e0<3>inEq.(2.12)).
       (4) PCOEFF(4)  = Nutrient-use coefficient for the production of microbe  #3 with anaerobic
               respiration. (en(3) in Eq. (2.12)).

       K. Subset 11: One variable (M/L3) is included in this record.
       (1) COFK = Inhibition coefficient.
20. ELEMENT (DISTRIBUTED) SOURCE/SINK FOR TRANSPORT SIMULATIONS
       Ten subsets of free-formatted data records are required in this data set.

       A. Subset 1: control parameters
       (1) NSEL = No. of source/sink elements.
       (2) NSPR = No. of source profiles, should be .GE. 1.
       (3) NSDP = No. of data points in each profile, should be .GE. 2.
       (4) KSAI = Is element-source/sink profile to be input analytically? 0 = no, 1 = yes.
                                              122

-------
       B. Subset 2: source/sink profile - This sub-data-set is needed if and only if NSEL .GT. 0. For each
               sub-data-record, NSDP of the data group (TSOSF(J,I),SOSF(J,I,1),SOSF(J,I,2)) are required.
               If this sub-data-record can be fitted in a line, we will need NSPR lines.
       (1) TSOSF(J,I) = Time of J-th data point in I-th profile, (T).
       (2). SOSF(J,I,1) = Source/sink flow rate of the J-th data point in the I-th profile, (L**3/T/L**3);
               positive for source and negative for sink.
       (3) SOSF(J,I,2) =  Source/sink concentration of the J-th data point in the I-th profile, (M/L**3).

       C. Subset 3: global source/sink element number. NSEL data points are required for this record.
       (1) LES(I) = Global element number of the I-th compressed distributed source/sink element.

       D. Subsets 4: Source type assigned to each element for microbe #1 - Usually one record per element.
               However, automatic generation can be  made.  For I-th (I = 1, 2, ...,) record, it contains the
               following.
       (1) MI = Global element number of the first element in the sequence.
       (2) NSEQ = NSEQ elements will be generated automatically.
       (3) MAD = Increment of element number for each of the NSEQ elements.
       (4) MITYP = Source type in element MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ elements.

       **** NOTE: A line with 5 O's is used to signal the end of this data set.

       E. Subset 5 ~ Subset NCC+3: Source type assigned to each element for microbe #2, microbe #3,
               substrate,  Oxygen, Nitrate, nutrient, and up to the NCC-th component. The input format is
               the same as subset 4.
21. POINT (WELL) SOURCE/SINK DATA FOR TRANSPORT SIMULATION
       Ten subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NWNP = No. of well or point source/sink nodes.
       (2) NWPR = No. of well or point source/sink strength profiles.
       (3) NWDP = No. of data points in each of the NWPR profiles.
       (4) KWAI = Is well-source/sink profile to be input analytically? 0 = no, 1 = yes.

       B. Subset 2: source/sink profiles - This group of data is needed if and only if NWNP .GT. 0. For each
               sub-data-record, NWDP of the data group (TWSSF(J,I), WSSF(J,I,1), WSSF(J,I,2)) are
               required.  If this  sub-data-record can be fitted in a line, we will need NWPR lines.
       (1) TWSSF(J,I) = Time of J-th data point in I-th profile, (T).
       (2) WSSF(J,I,1) = Source/sink flow rate of the J-th data point in the I-th profile, (L**3/T/L**3);
               positive for source and negative for sink.
       (3) WSSF(J,I,2) = Source/sink concentration of the J-th data point in the I-th profile, (M/L**3).

       C. Subset 3: global source/sink node number - This group of data is needed if and only if NWNP .GT.
               0. NWNP data points are required for this record.
       (1) NPW(I) = Global node number of the I-th compressed point source/sink node.
                                              123

-------
       D. Subset 4: Source type assigned to each well for microbe #1 - Usually one record per element.
              However,automatic generation can be made.
       (1) NI = Compressed point source/sink node number of the first node in a sequence.
       (2) NSEQ = NSEQ nodes will contain the source types that will be automatically generated.
       (3) NIAD = Increment of NI for each of the NSEQ nodes.
       (4) NITYP = Source type in node NI.
       (5) NTYPAD = Increment of NITYP for each of the NSEQ subsequent nodes.

       **** NOTE: A record with 5 O's must be used to signal the end of this data set.

       E. Subset 5 ~ Subset NCC+3: Source type assigned to each well for microbe #2, microbe #3,
              substrate, Oxygen, Nitrate, nutrient, and up to the NCC-th component.
22. RUN-IN/FLOW-OUT (VARIABLE) BOUNDARY CONDITIONS FOR TRANSPORT SIMULATIONS

       Ten subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NVES = No. of variable boundary element sides
       (2) NVNP = No. of variable boundary nodal points
       (3) NRPR = No. of incoming fluid concentration profiles to be applied to variable boundary element
              sides.
       (4) NRDP = No. of data points in each of the NRPR profiles.
       (5) KRAI = Is incoming concentration profile to be input analytically? 0 = no, 1 = yes.

       B. Subset 2: variable boundary flux profile - NRPR records are needed. Each record contains NRDP
              data points and is FREE-FORMATTED. Each data point has 2 numbers representing the
              time and run-in flow-out concentrations, respectively as follows:
       (1) TCRSF(J,I) = Time of the J-th data point on the I-th run-in concentration profile, (T),
       (2) CRSF(J,I) = Concentration of the J-th data point on the I-th profile, (M/L**3).

       C. Subset 3: Run-in concentration type assigned to each of all NVES sides for microbe #1. Usually
              one record per variable element side. However, automatic generation can be made. Each
              record contains 5 variable and is FREE-FORMATTED.
       (1) MI = Compressed VB element side of the first side in a sequence.
       (2) NSEQ = NSEQ subsequent sides will be generated automatically.
       (3) MIAD = Increment of MI for each of NSEQ subsequent sides.
       (4) MITYP = Type of concentration profile assigned to side MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ subsequent sides.

       **** NOTE: A record with 5 O's must be used to signal the end of this subdata set.

       D. Subset 4 ~ Subsets 4 through NCC+2: Run-in concentration type assigned to each of all NVES
              sides for microbe #2, microbe #3, substrate, Oxygen, Nitrate,  nutrient, and up to the NCC-th
              component. The input format is the same as subset 3.

       E. Subset NCC+3: global nodal number of all run-in flow-out boundary nodes Usually NVNP records

                                            124

-------
              are needed for this sub-data set. However, automatic generation can be made. Each record
              contains 5 variables and is FREE-FORMATTED.
       (1) NI = Compressed VB node number of the first node in a sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) NIAD = Increment for NI for each of the NSEQ nodes.
       (4) NODE = Global nodal number of the node NI.
       (5) NODEAD = Increment of NODE for each of the NSEQ subsequent nodes.

       **** NOTE: A record with 5 O's is used to signal end of this subdata set.

       F. Subset NCC+4: Specification of run-in boundary element sides - Normally, NVES records are
              required, one each for a VB element side. However, if a group of VB element sides appears
              in a regular pattern, automatic generation may be made. Each record contains 11 variables
              and is FREE-FORMATTED.
       (1) MI = Compressed VB element side number of the first side in a sequence.
       (2) NSEQ = NSEQ subsequent VB element sides will be generated automatically.
       (3) M = Global element no. to which the Ml-th VB side belongs.
       (4) IS1 = Global node number of the first node of element side MI.
       (5) IS2 = Global node number of the second node of element side MI.
       (6) MIAD = Increment of MI for each of the NSEQ subsequent VB element sides.
       (7) MAD = Increment of M for each of the NSEQ subsequent VB element sides.
       (8) IS IAD = Increment of IS1 for each of the  NSEQ subsequent element sides.
       (9) IS2AD = Increment of IS2 for each of the  NSEQ subsequent element sides.

       **** NOTE: A record with 9 O's is used to signal the end of this subdata set.
23. DIRICHLET BOUNDARY CONDITIONS FOR TRANSPORT SIMULATIONS

       Five subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1)  NDNP = No. of Dirichlet nodes, should be .GE. 1.
       (2)  NDPR = No. of Dirichlet profiles, should be .GE. 1.
       (3)  NDDP = No. of data points in each profile, should be .GE. 2.
       (4)  KDAI = Is Dirichlet boundary value profile to be input analytically? 0 = no, 1 = yes.

       B. Subset 2: Dirichlet-concentration profiles - NDPR  records are needed.  Each record contains
              NDDP data points and is FREE-FORMATTED. Each data point has 2 numbers representing
              the time and Dirichlet concentrations, respectively as follows:
       (1) TCDBF(J,I) = Time of J-th data point in I-th Dirichlet-concentration profile, (T).
       (2) CDBF(J,I) = Concentration of J-th data point in I-th Dirichlet-concentration profile, (M/L**3).

       C. Subset 3: Dirichlet concentration types assigned to Dirichlet nodes for Microbe # 1 - Normally one
              record per Dirichlet node, i.e., a total of NDNP records, is needed.  However, if the Dirichlet
              nodes appear in a regular pattern, automatic generation may be made. Each record contains
              5 variables and is FREE-FORMATTED.
       (1)  NI = Compressed Dirichlet node number of the first node in the sequence.

                                            125

-------
       (2)  NSEQ = NSEQ nodes will contain the Dirichlet concentration types that will be automatically
              generated.
       (3)  NIAD = Increment of NI for each of the NSEQ nodes.
       (4)  NITYP = Dirichlet concentration type in node NI.
       (5)  NTYPAD = Increment of NITYP for each of the NSEQ subsequent nodes.

       **** NOTE:  A record with 5 O's must be used to signal the end of this data set.

       D. Subset 4 ~ Subset NCC+2: Dirichlet concentration type assigned to each of all NDNP points for
              microbe #2,  microbe  #3, substrate, Oxygen, Nitrate,  nutrient,  and up to the NCC-th
              component. The input  format is the same  as subset.

       E. Subset NCC+3: global node number of compressed Dirichlet nodes - One record is needed for this
              sub-data set, which contains NDNP variables and is FREE-FORMATTED.
       (1)  NPDB(I) = Global node number of the I-th compressed Dirichlet node.
24. CAUCHY BOUNDARY CONDITIONS FOR TRANSPORT SIMULATION

       Six subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NCES = No. of Cauchy element sides.
       (2) NCNP = No. of Cauchy nodal points.
       (3) NCPR = No. of Cauchy-flux profiles.
       (4) NCDP = No. of data points on each Cauchy-flux profile.
       (5) KCAI = Is Cauchy flux profile to be input analytically? 0 = no, 1 = yes.


       B. Subset 2: Cauchy flux profiles - NCPR records are needed.  Each record contains NCDP data
              points and is FREE-FORMATTED. Each data point has 2 numbers representing the time and
              Cauchy flux, respectively, as follows:
       (1) TQCBF(J,I) = Time of the J-th data point in the I-th Cauchy flux profile, (T).
       (2) QCBF(J,I) = Value of Cauchy flux of the  J-th data point in the I-th Cauchy flux profile,
       (M/T/L**2).

       C. Subset 3: Cauchy flux type assigned to each of all NCNP nodes -Usually one record per Cauchy
              node. However, automatic generation can be made.  Each record contains 5 variables and is
              FREE-FORMATTED.
       (1) MI = Compressed Cauchy boundary node number of the first node in a sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) MIAD = Increment of MI for each of NSEQ subsequent nodes.
       (4) MITYP = Type of Cauchy flux profile assigned to node MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ subsequent nodes.

       **** NOTE: A record with 5 O's must be used to signal the end of this subdata set.

       D. Subset 4 ~ Subset NCC+2: Cauchy flux type assigned to each of all NCNP nodes for microbe #2,

                                           126

-------
              microbe #3, substrate, Oxygen, Nitrate, nutrient, and up to the NCC-th component. The input
              format is the same as subset 3.

       E. Subset NCC+3: global nodal number of all Cauchy boundary nodes -Usually NCNP records are
              needed for this sub-data set. However, automatic generation can be made. Each record
              contains 5 variables and is FREE-FORMATTED.
       (1) NI = Compressed Cauchy boundary node number of the first node in a sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) NIAD = Increment for NI for each of the NSEQ nodes.
       (4) NODE = Global nodal number of the node NI.
       (5) NODEAD = Increment of the global nodal number for each of the NSEQ subsequent nodes.

       **** NOTE: A record with 5 O's is used to signal end of this subdata set.

       F. Subset NCC+4: specification of Cauchy boundary element sides -Normally, NCES  records are
              required, one each for a Cauchy boundary element side.  However, if a group of Cauchy
              element sides appears in a regular pattern, automatic generation may be made. Each record
              contains 9 variables and is FREE-FORMATTED.
       (1) MI = Compressed Cauchy boundary element side number of the first element side in a sequence.
       (2) NSEQ = NSEQ subsequent Cauchy boundary element sides will be generated automatically.
       (3) M = Global element number to which the Ml-th Cauchy element side belongs.
       (4) IS1 = Compressed Cauchy node number of the first node of element side MI.
       (5) IS2 = Compressed Cauchy node number of the second node of element side MI.
       (6) MIAD = Increment of MI for each of the NSEQ subsequent Cauchy boundary element sides.
       (7) MAD = Increment of M for each of the NSEQ subsequent Cauchy boundary element sides.
       (8) IS IAD = Increment of IS1 for each of the NSEQ subsequent element sides.
       (9) IS2AD = Increment of IS2 for each of the NSEQ subsequent element sides.

       **** NOTE: A record with 9 O's is used to signal the end of this subdata set.
25. NEUMANN BOUNDARY CONDITIONS FOR TRANSPORT SIMULATIONS

       Six subsets of data records are required for this data set.

       A. Subset 1: control parameters
       (1) NNES = No. of Neumann element sides.
       (2) NNNP = No. of Neumann nodal points.
       (3) NNPR = No. of Neumann-flux profiles.
       (4) NNDP = No. of data points on each Neumann-flux profile.
       (5) KSAI = Is Neumann flux profile to be input analytically? 0 = no, 1 = yes.

       B. Subset 2: Neumann flux profiles - NNPR records are needed.  Each record contains NNDP data
              points and is FREE-FORMATTED. Each data point has 2 numbers representing the time and
              Neumann flux,  respectively, as follows:
       (1) TQNBF(J,I) = Time of the J-th data point in the I-th Neumann flux profile, (T).
       (2) QNBF(J,I) = Value of Neumann flux of the  J-th data point in the I-th Neumann flux profile,
              (M/T/L**2).

                                           127

-------
       C. Subset 3: Neumann flux type assigned to each of all NNNP nodes -Usually one record per
              Neumann boundary node.  However, automatic generation can be made.  Each record
              contains 5 variables and is FREE-FORMATTED.
       (1) MI = Compressed Neumann boundary node nmber of the first node in a sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) MIAD = Increment of MI for each of NSEQ subsequent nodes.
       (4) MITYP  = Type of Neumann flux profile assigned to node MI.
       (5) MTYPAD = Increment of MITYP for each of the NSEQ subsequent nodes.

       **** NOTE: A record with 5 O's must be used to signal the end of this sub-data set.

       D. Subset 4 ~ Subset NCC+2: Neumann flux type assigned to each of all NNNP nodes for microbe
              #2, microbe #3, substrate, Oxygen, Nitrate, nutrient, and up to the NCC-th component. The
              input format is the same as subset 3.

       E. Subset NCC+3: global nodal number of all Neumann boundary nodes -Usually NNNP records are
              needed for this sub-data set.  However, automatic  generation can be made. Each record
              contains 5 variables and is FREE-FORMATTED.
       (1) NI = Compressed Neumann boundary node number of the first node in a sequence.
       (2) NSEQ = NSEQ subsequent nodes will be generated automatically.
       (3) NIAD = Increment for NI for each of the NSEQ nodes.
       (4) NODE = Global nodal number of the node NI.
       (5) NODEAD = Increment of the global nodal number for each of the NSEQ subsequent nodes.

       **** NOTE: A record with 5 O's is used to signal end of this subdata set.

       F. Subset NCC+4: specification of Neumann boundary element sides -Normally, NNES records are
              required, one each for a Neumann boundary element side. However, if a group of Neumann
              element sides appears in a regular pattern, automatic generation may be made. Each record
              contains 9 variable and is FREE-FORMATTED.
       (1)  MI = Compressed Neumann boundary element side number of the first element side in a
       sequence.
       (2) NSEQ = NSEQ subsequent Neumann boundary element sides will be generated automatically.
       (3) M = Global element no. to which the Ml-th Neumann boundary side.
       (4) IS1 = Compressed Neumann node number of the first node of element side MI.
       (5) IS2 = Compressed Neumann node number of the second node of element side MI.
       (6) MIAD = Increment of MI for each of the NSEQ subsequent ides.
       (7) MAD =  Increment of M for each of the NSEQ subsequent Neumann boundary sides.
       (8) IS IAD = Increment of IS1 for each of the NSEQ subsequent element sides.
       (9) IS2AD = Increment of IS2 for each of the NSEQ subsequent element sides.

       **** NOTE: A record with 9 O's is used to signal the end of this subdata set.
26. HYDROLOGICAL VARIABLES

       This data set is needed if and only if KVIT .LE. 0.  When KVIT .LE. 0, three sub-data sets are
       needed, two groups for the velocity field and the other group for the moisture content.

                                           128

-------
       A. Subset 1: velocity field (x-direction) - Usually NNP records are needed.  However, if velocity
              appears in a regular pattern, automatic generation can be made. Each record contains 6
              variables and is FREE-FORMATTED.
       (1) NI = Node number of the first node in a sequence
       (2) NSEQ = NSEQ subsequent nodes will be automatically generated,
       (3) NIAD = Increment of node number in each of the NSEQ subsequent nodes,
       (4) VXNI = x-velocity component at node NI, (L/T),
       (5) VXAD = Increment of VXNI for each of the NSEQ subsequent nodes, (L/T),
       (6)  0.0

       B. Subset 2: velocity field  (z-direction) - Usually NNP records are needed.  However, if velocity
              appears in a regular pattern, automatic generation can be made. Each record contains 6
              variables and is FREE-FORMATTED.
       (1) NI = Node number of the first node in a sequence
       (2) NSEQ = NSEQ subsequent nodes will be automatically generated,
       (3) NIAD = Increment of node number in each of the NSEQ subsequent nodes,
       (4) VZNI = z-velocity component at node NI, (L/T),
       (5) VZAD = Increment of VZNI for each of the NSEQ subsequent nodes, (L/T),
       (6)  0.0

       **** NOTE: A record with 6 O's is used to signal the end of this sub-data set.

       C. Subset 3: moisture content field - Usually, NEL records are needed.  However, if moisture content
              appears in a regular pattern,  automatic generation can be made. Each record contains 5
              variables and is FREE-FORMATTED.
       (1) MI = Element number of the first element in a sequence,
       (2) NSEQ = NSEQ subsequent elements will be automatically generated,
       (3) MIAD = Increment of MI for each of NSEQ subsequent elements,
       (4) THNI = Moisture content of element NI, (Decimal point),
       (5) THNIAD = Increment of THNI for NSEQ subsequent elements, (Decimal point).
       (6) 0.0

       **** NOTE: A record with 6 O's is used to signal the end of this subdata set.
27.  END OF JOB
       If another problem is to be run, then input begins again with input data set 1. If termination of the job
       is desired, a blank card must be inserted at the end of the data set.
                                            129

-------
130

-------
                                        REFERENCES
Bachelor, G. A., D. E. Cawlfield, F. T. Lindstrom, and L. Boersma, Denitrification in nonhomogeneous
       laboratory scale aquifers: 5: User's manual for the mathematical model LT3VSI, A draft report, 1990.

Benfield, L. D., and F. J. Molz, A model for the activated sludge process which considers wastewater
       characteristics, flux behavior, and microbial population, Biotechnol. Bioeng., 26, 352-361, 1984.

Cheng, J. R., H. P. Cheng, and G. T. Yeh. A Lagrangian-Eulerian Method with Adaptively Local Zooming
       and Peak/Valley Capturing Approach to Solve Two-Dimensional Advection-Diffusion Transport
       Equations. International J. of Numerical Methods in Engineering.  1995 (in press).

Freeze, R. A., Role of subsurface flow in generating surface runoff: 1. Base flow contribution to channel flow,
       Water Resour. Res., 8, 609-623, 1972a.

Freeze, R. A., Role of subsurface flow in generating surface runoff: 2. Upstream source areas, Water Resour.
       Res., 8, 1272-1283, 1972b.

Herbert, D., Some principles of continuous culture, in Recent Progress in Microbiology, edited by G.
       Tunevall, Blackwell Scientific Publishers, Oxford, England, 1958.

Huyakorn, P. S., E. P. Springer, V. Guvanasen, and T. D. Wadsworth, A three-dimensional finite-element
       model for simulating water flow in variably saturated porous media, Water Resources Research, Vol.
       22, No. 13, 1790-1808, 1986.

MacQuarrie, K. T. B. and E. A. Sudicky, Simulation of biodegradable organic contaminants in groundwater,
       2. plume behavior in uniform and random flow fields, Water Resour. Res., 26(2), 223-239, 1990.

Molz, F. J. M. A. Widdowson, and L. D. Benfield, Simulation of microbial growth dynamics coupled to
       nutrient and oxygen transport in porous media, Water Resour. Res., 22(8), 1207-1216, 1986.

van Genuchten, M. Th., A closed form equation for predicting the hydraulic conductivity of unsaturated soils,
       Soil Science Society of America Journal, 44, 892-898, 1980.

Widdowson, M. A., F. J. Molz, and L. D. Benfield, A numerical transport model for oxygen- and nitrate-based
       respiration linked to substrate and nutrient availability in porous media, Water Resour. Res., 24(9),
       1553-1565, 1988.

Yeh, G. T., and  D.  S. Ward,  FEMWATER: A finite-element model  of water flow through saturated-
       unsaturated porous media, Rep. ORNL-5567, Oak Ridge Nat. Lab., Oak Ridge, Term., 37831, 137
       pp., 1980.

Yeh, G. T., and D. S. Ward, FEMWASTE: A finite-element model of waste transport through saturated-
       unsaturated porous media, Rep. ORNL-5601, Oak Ridge Nat. Lab., Oak Ridge, Term., 37831, 137
       pp., 1981.
                                              131

-------
Yeh, G. T., FEMWATER: A finite element model of water flow through saturated-unsaturated porous media,
       First Revision, Rep. ORNL-5567/R1, Oak Ridge Nat. Lab., Oak Ridge, Tenn., 37831, 258 pp., 1987.

Yeh, G. T., A Lagrangian-Eulerian method with zoomable hidden fine mesh approach to solving advection-
       dispersion equations, Water Resources Research Vol. 26, No. 6, 1133-1144, 1990.

Yeh, G. T., J. R. Chang, and T. E. Short, An exact peak capturing and oscillation-free scheme to solve
       advection-dispersion transport equations, Water Resour. Res., 28(11), 2937-2951, 1992.

Yeh, G. T., Class notes: CE597C: Computational Subsurface Hydrology Part II, The Pennsylvania State
       University, University Park, PA., 16802, Spring semester 1992.

Yeh, G. T., 3DFEMWATER: A three-dimensional finite element model of WATER flow through saturated-
       unsaturated media:  Version 2.0,  Short  course notes of  Simulation of Subsurface Flow  and
       Contaminant Transport by Finite Element and Analytical Methods, The Pennsylvania State University,
       University Park, PA., 16802, May 1993.

Yeh, G. T., 3DLEWASTE: A three-dimensional hybrid Lagrangian-Eulerian finite element model of WASTE
       transport through saturated-unsaturated media: Version 2.0, Short course notes of Simulation of
       Subsurface Flow and Contaminant Transport by Finite  Element and Analytical  Methods,  The
       Pennsylvania State University, University Park, PA., 16802, May 1993.

Yeh, G. T., J. R. Chang, J. P. Gwo, H. C. Lin, D. R. Richards, and W. D. Martin, 3DSALT: A three-
       dimensional finite  element model of density-dependent flow and transport  through saturate-
       unsaturated media, Instruction Report HL-94-1, US Army Corps of Engineers, 1994.
                                             132

-------