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 h ------- 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 |