COUPLING TRANSPORT CODES WITH GEOCHEMICAL MODELS Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. ____________________________________________ Anjani Kumar Certificate of Approval: ________________________ ___________________________ Mark O. Barnett T. Prabhakar Clement, Chair Associate Professor Associate Professor Civil Engineering Civil Engineering ________________________ ________________________ Dongye Zhao Stephen L. McFarland Associate Professor Acting Dean Civil Engineering Graduate School COUPLING TRANSPORT CODES WITH GEOCHEMICAL MODELS Anjani Kumar A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Master of Science Auburn, Alabama December 15th, 2006 iii COUPLING TRANSPORT CODES WITH GEOCHEMICAL MODELS Anjani Kumar Permission is granted to Auburn University to make copies of this thesis at its discretion, upon request of individuals or institutions and at their expense. The author reserves all publication rights. ________________________ Signature of Author ________________________ Date of Graduation iv THESIS ABSTRACT COUPLING TRANSPORT CODES WITH GEOCHEMICAL MODELS Anjani Kumar Master of Science, December 15, 2006 (M.S., Indian Institute of Science, Bangalore (INDIA), 2003) (B.S., Indian Institute of Technology, Roorkee (INDIA), 1999) 93 Typed Pages Directed by T. Prabhakar Clement Multi-dimensional reactive transport models are needed for understanding the fate and transport of the metal contaminants in groundwater aquifers. The goal of this research is to develop a framework for coupling the MINTEQ family of geochemical equilibrium models with the public domain reactive transport code RT3D. In this research effort, a suite of one-dimensional reactive transport models were built by coupling MICROQL, MINEQL and MINTEQA2 equilibrium geochemistry modules with a fully implicit, finite-difference, advection-dispersion solver. The one-dimensional reactive transport models were tested by solving four benchmark problems, which were identified as a part of this study. These benchmark problems simulate four distinct equilibrium-controlled processes that involve ion-exchange, surface-complexation, precipitation-dissolution, and redox reactions. The modular framework developed for coupling the one-dimensional codes with the equilibrium codes was used to implement the USEPA code MINTEQA2 within RT3D to develop a prototype model. The v prototype model was also tested by solving several benchmark problems. Finally, a scalable numerical model was developed to describe the reactive transport of arsenic in MnO2(s) containing soils. The model performance was validated by reproducing a set of column-scale data using scaled parameters obtained from the batch experiments. vi ACKNOWLEDGEMENTS I would like to express my gratitude to my research advisor Dr. T. Prabhakar Clement for his encouragement, tolerance and guidance during my master?s studies. I also highly value the patience and support he demonstrated during my trying times here. I would like to extend my special thanks to Dr. Barnett and Dr. Zhao for agreeing to be in my thesis committee. I learnt a lot from Dr. Mark Barnett in the area of basic chemistry, which ultimately formed the foundation for my thesis work. I am thankful to him for being a pleasant and challenging instructor. It was sheer joy to go through the two courses taught by Dr. Jacob Dane. Thanks a lot. I am also grateful to other faculty as well as staff members of my civil engineering department, office of international education and Auburn University for making my stay comfortable here for close to three years. I would like to acknowledge the precious company of my research group Matt, Gautham, Venkat, Rohit, Linzy, Vijay, Tanja and other graduate students in my dept. I also thank researchers at other universities and labs from whom I received timely technical support, in particular Dr. Henning Prommer at CSIRO (Australia). This study was financially supported by a grant from the Sustainable Water Resources Research Center of 21st Century Frontier Research Program, Korean Ministry of Science and Technology. This thesis is dedicated to Sri Chaitanya Mahaprabhu and His Associates. I reserve my special gratitude and affection for the members of Gita-class and Srimad- Bhagavatam program. I am grateful to my family members for their love. vii Style manual or journal used Auburn University Graduate School Guide to Preparare Master?s Thesis Computer software used Microsoft Word, Microsoft Excel, End Note 9.0 (Water Resources Research) viii TABLE OF CONTENTS LIST OF TABLES ?????.??????????????..???????.. x LIST OF FIGURES ????...?????????????????...???? xi CHAPTER ONE INTRODUCTION AND LITERATURE SURVEY 1.1 Background ???????????????????? ?????. 1 1.2 Equilibrium Geochemistry Models ??????????.........?.??? 2 1.3 Reactive Transport Models ??????????????..????? 5 1.4 Research Objectives and Organization of the Thesis ???? ?????10 CHAPTER 2 MODELING VARIOUS TYPES OF EQUILIBRIUM REACTIONS 2.1 Introduction ????????????????????????? 13 2.2 General description of the problem of chemical equilibrium in aqueous systems ?????????????????????????.?. 13 2.3 Modeling Ion-exchange reaction ???...????????????? 16 2.4 Modeling Surface-complexation reaction ?????????????.. 17 2.5 Modeling Precipitation-dissolution ?????? ????????.? 20 2.6 Modeling Redox reaction??????????????????..?.. 23 2.6.1. Effective internal approac??????????????..?. 23 2.6.2. External electron approach ????????????........? 23 2.6.3. Oxygen fugacity approach ???????????.???.... 23 2.6.4. Redox couple approach ????????????.????. 24 2.6.5 Typical techniques for dealing with convergence problems 2.6.5.1 Basis-switching ???????????????? 25 2.6.5.2 Log-based formulation and its application???...?? 26 2.7 Summary ?????????????????????????.... 27 CHAPTER 3 MODELING REACTIVE TRANSPORT BY COUPLING MICROQL/MINEQL TO TRANSPORT CODES 3.1 Introduction ???????????????? ????? ??? 28 3.2 Reactive transport equations and their general solution strategy ???? 29 3.2.1. Equations representing the reactive transport . ?????.?.? 29 3.2.2 Solution of ADRE ????????????????? 30 ix 3.3 Reaction 1: Ion-Exchange ???????????????????. 32 3.4. Reaction 2: Surface-complexation ??????????????.? 35 3.5. Reaction 3: Precipitation-dissolution ?????????????.? 40 3.6. Reaction 4: Redox ??????????????????.?..?... 44 3.7 Summary ??????????????????????..?....... 48 CHAPTER 4 MODELING REACTIVE TRANSPORT BY COUPLING MINTEQA2 TO RT3D 4.1 Introduction ????????????????????????. . 49 4.2 Three Dimensional Reactive Transport Equations and the General Solution Strategy ????????????????????????..?. .. 50 4.2.1 Reactive Transport Equations ?????????????? . 50 4.2.2 Integration of Geochemistry Transport Routines ??.???..?. 51 4.3 Coupling of RT3D and MINTEQA2 Routines ????????.???. 52 4.4 Example Problem ????????????????????.?? . 55 4.5 Summary ??????.?????????????????.... ?.. 60 CHAPTER 5 DEVELOPMENT OF A SCALABLE MODEL FOR PREDICTING ARSENIC OXIDATION AND ADSORPTION AT PYROLUSITE SURFACES 5.1 Introduction ????????????????????????? 61 5.2 Background ???????????????????????..?. 61 5.3 Results from the batch experiments ??????????????..? 62 5.3.1 Results from the oxidation kinetics and adsorption experiments ... 62 5.4 Reactive transport modeling and designing the column experiments ?? 64 5.4.1 Conceptual model ??????????????..???.. 64 5.4.2 Reactive transport model ????????????..??... 65 5.4.3 Solution of the reactive transport equations ????????? 65 5.4.4 Moving from batch to column scale simulations ?????...?. 67 5.5 Results and discussion ?????????????..?????..?... 70 5.6 Summary ???????????????????????..??. 71 CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS 6.1 Conclusions ?????????.???????????????.73 6.2 Recommendations for future work ???????????????? 74 REFERENCES ???????????????????????????? 76 x LIST OF TABLES Table 2. 1 : An example of Basis-switching..................................................................... 25 Table 3 . 1 : Physical parameters used in the simulation.................................................. 33 Table 3 . 2 : Initial and boundary conditions used in the simulation................................ 33 Table 3 . 3 : Physical parameters used in the simulation.................................................. 36 Table 3 . 4 : Chemical sorption parameters used in the simulation .................................. 37 Table 3 . 5 : Concentrations of components for the test cases: Initial and Boundary conditions.................................................................................................................. 37 Table 3 . 6: Stoichiometry of Reactions............................................................................ 38 Table 3 . 7 : Physical parameters used in the simulation.................................................. 41 Table 3 . 8 : Initial and boundary conditions used in the simulation................................ 41 Table 3 . 9 : Physical parameters used in the simulation.................................................. 45 Table 3 . 10 : Initial and boundary conditions used in the simulation.............................. 45 Table 3 . 11 : Summary of the model development work in chapter 3............................. 48 Table 4 . 1 : Physical parameters used in the simulation.................................................. 56 Table 4 . 2 : Initial and boundary conditions used in the simulation................................ 56 Table 5 . 1 : Model data common to all simulations......................................................... 63 xi LIST OF FIGURES Figure 2 . 1 : MICROQL solution process for chemical equilibrium problem in aqueous systems........................................................................................................ 15 Figure 2 . 2 : MICROQL solution process for chemical equilibrium problem in aqueous systems, involving surface-complexation reactions. .................................. 19 Figure 2 . 3 : MINEQL solution process for chemical equilibrium problem in aqueous systems, involving precipitation-dissolution reactions............................................. 22 Figure 3. 1 : Flow chart for generalized solution process of a reactive transport problem. ............................................................................................................................ 31 Figure 3 . 2: Breakthrough profile of Na+......................................................................... 34 Figure 3 . 3: Breakthrough profile of Mg2+....................................................................... 34 Figure 3 . 4 : Breakthrough profile of Ca2+....................................................................... 35 Figure 3 . 5 : Profiles obtained at t=15hrs during simulation for case I . ......................... 39 Figure 3 . 6 : Profiles obtained at t=15hrs during simulation for cases II and V.............. 39 Figure 3 . 7 : Profiles obtained at t=15hrs during simulation for cases III and VI. .......... 40 Figure 3 . 8 : Calcite and Dolomite profiles at 21,000 sec................................................ 42 Figure 3 . 9 : pH profile at 21,000 sec. ............................................................................. 43 Figure 3 . 10 : Aqueous concentration profiles at 21,000 sec........................................... 43 Figure 3 . 11 : Aqueous concentration profiles at 21,000 sec........................................... 44 Figure 3 . 12 : pH profile at 2.0E6 s.................................................................................. 46 Figure 3 . 13 : pE profile at 2.0E6 s.................................................................................. 47 Figure 3 . 14 : Mineral Profiles at 2.0E6 seconds............................................................. 47 Figure 4 . 1 : GMS user interface for the new geochemistry module????????54 Figure 4 . 2 : Concentration profile of calcite mineral along the column......................... 57 Figure 4 . 3 : Concentration profile of gibbsite mineral along the column....................... 58 xii Figure 4 . 4 : Predicted pH profile along the column........................................................ 58 Figure 4 . 5 : Concentration profile of Ca2+ (aqueous) along the column. ....................... 59 Figure 4 . 6 : Concentration profile of SO42- (aqueous) along the column ....................... 59 Figure 5 . 1 : For the column fully packed with MnO2(s), comparing the results of As(III) (Exp A) and As(V) (Exp B) transport with the model predictions............... 69 Figure 5 . 2 : For a column with MnO2(s) 50% by wt. in a sand-MnO2(s) mixture, comparison of experimental and model results for As(III) transport (Exp C),......... 69 Figure 5 . 3 : For a column with very low MnO2(s) content, comparison of experimental and model results for As(III) transport (Exp D); partial oxidation of As(III) ....................................................................................................................... 70 1 CHAPTER ONE INTRODUCTION AND LITERATURE SURVEY 1.1 BACKGROUND About 95% percent of the world?s available fresh water is stored as groundwater [Freeze and Cherry, 1979] . Currently, groundwater accounts for approximately twenty one percent of annual water supply within United States [USGS, 2000]. Therefore, the possibility of contamination of groundwater aquifers by toxic chemicals poses considerable risk to an important resource. Groundwater systems can be contaminated by several possible sources. These sources include leachates from mine tailings and landfills that may contain a number of toxic chemical wastes. These toxic wastes may include heavy metals such as arsenic, cadmium, chromium, copper, iron, mercury, molybdenum, nickel, lead, selenium, uranium and zinc. The cost of groundwater remediation efforts can be enormous. According to an USDOE estimate, cleanup of hazardous waste sites, using currently available technologies may take over 70 years at an estimated total cost of 300 billion dollars [Frazier and Johnson, 2005]. Understanding the full range of the physical and chemical processes occurring during the contamination of the groundwater environment will help implement efficient remediation techniques and minimize the overall costs. Computer models that can simulate physical transport processes coupled to the geochemical 2 processes can play a key role in developing this understanding. Since these models simulate simultaneous transport and reaction of several chemical components, they are commonly referred to as multi-component reactive transport models. The focus of this work is to develop a set of new computer codes that can simulate multi-component reactive transport in porous media. 1.2 EQUILIBRIUM GEOCHEMISTRY MODELS Geochemical reactions can be mediated by either kinetically-controlled or equilibrium- controlled processes. The first part of the thesis (up to chapter 4) deals exclusively with equilibrium reactions. Therefore, before describing reactive transport models, the details of some widely-used equilibrium geochemistry models are reviewed. The modeling of equilibrium geochemistry of natural waters is a well established field. Computer models for solving equilibrium geochemistry problems that have been widely used include WATEQ [Truesdell and Jones, 1974], MICROQL [Westall, 1979a; 1979b], MINEQL [Westall J. C., et al., 1976], MINTEQA2 [Allison, et al., 1990], PHREEQE [Parkhurst, et al., 1980], PHREEQC [Parkhurst, 1995; , 1999], and EQ3/6 and its derivatives [Wolery, 1992]. WATEQ [Truesdell and Jones, 1974] was one of the first widely used geochemical equilibrium models. WATEQ uses a simple iteration scheme to solve aqueous speciation problems. While WATEQ cannot explicitly handle heterogeneous reactions such as precipitation-dissolution, it can compute mineral saturation indices. The model MINEQL was developed by Westall et al. [1976] by modifying REDEQL [McDuff and Morel, 1972]. MINEQL employs a tableau approach that reduces 3 the chemical equilibrium problem into a set of nonlinear equations, that can be solved by an iterative process. MINEQL uses the Newton-Raphson iteration scheme that systematically improves the rate of convergence and is superior to the simpler substitution method employed by WATEQ. MINEQL is also more advanced than WATEQ, since it can handle reactions such as mineral precipitation and dissolution reaction. MICROQL ? I [Westall, 1979a] is another chemical equilibrium program, but was written in the Basic language. Also MICROQL -1 computes chemical equilibrium in aqueous systems without sorbed or solid phase species. MICROQL ? II [Westall, 1979b] includes additional computation routines for solving adsorption equilibria in aqueous systems, using constant capacitance, diffuse layer, stern layer and triple layer adsorption model. The chemical speciation model MINTEQA2 is a more versatile, state-of-the-art equilibrium chemistry program. MINTEQA2 is a geochemical equilibrium speciation model developed to analyze dilute aqueous systems. The original version of the MINTEQ [Felmy, et al., 1983] was developed at the Battelle Pacific Northwest Laboratory (PNL) by combining the fundamental mathematical structure of MINEQL [Westall J. C., et al., 1976] with a thermodynamic database of the USGS model WATEQ3 [Ball, et al., 1981]. MINTEQA2 is substantially different from the original MINTEQ code and includes comprehensive thermodynamic database. MINTEQA2 is complemented by PRODEFA2, an interactive program used for creating the input files. The original PRODEF was also a product of PNL but it was substantially modified to develop PRODEFA2. MINTEQA2 can be used to calculate the equilibrium composition of dilute aqueous solutions in the laboratory or natural systems. It can be used to calculate the mass distribution between 4 the dissolved, adsorbed, and multiple solid phases under a variety of conditions that can include a gas phase maintained at a constant partial pressure. Compared to MINEQL, MINTEQA2 has a larger thermodynamic database (based on the WATEQ database) and has many aqueous complexes as well as mineral phases. MINTEQA2 also offers a number of surface complexation models. Another useful capability is the flexibility of running with different redox couplesin on or off modes. The PHREEQE computer program and its derivatives PHREEQC [Parkhurst, 1995; , 1999] were developed by the U.S. Geological Survey. PHREEQE is written in the FORTRAN programming language and is designed to perform a variety of aqueous geochemical calculations. The current version of PHREEQE, known as PHREEQC is written in the programming language C. PHREEQC version-1 is based on an ion- association aqueous model and has capabilities for (1) speciation and saturation-index calculations; and (2) batch-reaction and one-dimensional (1D) transport calculations involving reversible reactions, which include aqueous, mineral, gas, solid-solution, surface-complexation, ion-exchange equilibria, and irreversible reactions. PHREEQC can model specified mole transfers of reactants, kinetically controlled reactions, mixing of solutions, and temperature changes. Furthermore, inverse modeling routines which can evaluate sets of mineral and gas mole transfers that account for differences in composition of different water samples are also available within PHREEQC. Additional features in PHREEQC version-2 include capabilities to simulate dispersion (or diffusion) and stagnant zones in 1-D transport calculations, model kinetic reactions with user- defined rate expressions, model the formation or dissolution of ideal, multi-component, or non-ideal binary solid solutions, model fixed-volume gas phases in addition to 5 fixed-pressure gas phases, allow the number of surface or exchange sites to vary with the dissolution or precipitation of minerals or kinetic reactants, and methods to automatically use multiple sets of convergence parameters. PHREEQC can be used as an all-purpose geochemical model for modeling groundwater systems. The current version, PHREEQC- 2, is versatile in that it can solve a wide range of problems including those involving surface chemistry and reaction kinetics. It can access thermodynamic data from several large, well-established databases or can use specified data. Other known geochemistry models are CHESS [van der Lee, 1998] and The Geochemist's Workbench [Bethke, 1994]. CHESS stands for ?CHemical Equilibrium of Species and Surfaces? and is used mainly for safety and performance assessment of waste repositories. Model was specifically developed to simulate the equilibrium state of complex aquatic systems including oxides or minerals, organics, colloids and gases. About seven different thermodynamic databases are available within CHESS and new species can be easily added. The Geochemist's Workbench [Bethke, 1994] includes a set of tools that can be used for manipulating chemical reactions, calculating stability diagrams, equilibrium states of natural waters, tracing reaction processes, modeling reactive transport, and for plotting the results of these calculations. Among all these equilibrium geochemistry models, MINTEQA2 and PHREEQC are most popular and widely used codes. 6 1.3 REACTIVE TRANSPORT MODELS WITH EQUILIBRIUM/KINETICS GEOCHEMISTRY There are several reactive transport models available in the literature that couple a transport model with an equilibrium geochemistry routine. The transport of solutes in groundwater systems is described by a set of linear partial differential equations where as the chemical equilibria is described by a set of nonlinear algebraic equations. Reactive transport models use different types of numerical approaches for coupling these two processes. Fundamentally, there are three distinct approaches available for coupling a transport model with geochemical equilibrium model: (1) mixed differential and algebraic equations approach (DAE), (2) direct substitution approach (DSA), and (3) sequential iterative/noniterative approach (SIA/SNIA). An extremely important consideration in any of these coupling approaches is the choice of primary dependent variables (PDVs). There have been six types of PDVs employed in the existing models: (1) concentration of all species, (2) concentration of all component species and precipitated species, (3) total analytical concentrations of aqueous components, (4) total dissolved concentrations of aqueous components, (5) concentration of aqueous component species, and (6) hybrid concentrations. In the DAE approach, the transport and chemical equilibrium equations are solved simultaneously. In this approach, PDVs are either (1) the concentration of all species [Miller and Benson, 1983] or (2) the concentrations of all components and precipitated species [Lichtner, 1985]. In DSA, the chemical equilibrium equations are substituted into the transport equations to result in a set of nonlinear PDEs, which are subsequently solved simultaneously for a set of PDVs. Several models have been proposed, each using 7 a different set of variables as the PDVs. Jennings et al. [1982] used the concentrations of aqueous component species as PDVs. Rubin and James [1973] and Valocchi et al. [1981] used total dissolved concentrations of aqueous components as PDVs. Yeh and Tripathi [1988] used the total analytical concentrations of aqueous components as the PDVs. Among others, DSA has been employed by Miller and Benson [1983], Carnahan [1988] and Steefel and Lasaga [1990] . Some of the investigators (Rubin [1983]; Lewis [1987]) have used Hybrid PDVs. Unlike other natural PDVs such as concentrations of all species or concentrations/activities of aqueous components, these hybrid PDVs are linear combinations of species concentrations. In the SIA/SNIA approach, the solution process is divided into two steps: transport step and reaction step. The aqueous components are transported individually by advection and dispersion and then the aqueous, solid and sorbed components are allowed to react with each other. These two steps are coupled either iteratively (SIA) or sequentially (SNIA). Various operator-splitting approaches for coupling theses two steps are discussed by Herzer and Kinzelbach [1989] and Yeh and Tripathi [1991]. Theis, et al [1982], Kirkner and Reeves [1988] and Kirkner, et al.[1985] used SIA approach with total dissolved concentrations of aqueous components as PDVs. Walsh [1984], Cederberg, et al [1985], Bryant, et al. [1986] and Walter, et al. [1994] used total analytical concentrations of aqueous components as PDVs. The SNI approach was used by Walter, et al.[1994], Hundsdorfer and Verwer [1995] and Appelo, et al. [1997] . Yeh and Tripathi [1989] critically evaluated all the approaches available for coupling transport to the chemistry and the use of several types of PDVs employed in the reactive transport models. Their first goal was to assess the computational burden posed 8 by the models. They found that both DAE and DSA models require excessive CPU memory and time for realistic two- and three-dimensional problems, and SIA models require the least amount of computer resources. Their second goal was to assess the flexibility of the model to represent various types of processes. They concluded that only those models that use the first three types of PDVs can treat the complete set of geochemical reactions simultaneously. Their third goal was how easily a model can be modified to deal with mixed kinetic and equilibrium reactions. They observed that DAE and SIA models can be modified with reasonable ease to handle mixed chemical kinetics and equilibria. DSA models require strenuous efforts during modification for treating mixed chemical kinetics and equilibria. They concluded that DSA and DAE models should remain as research tools for solving one-dimensional problems and recommended SIA models with the third type of PDVs for solving practical problems. In literature, a number of reactive transport models are available that couple the equilibrium speciation models PHREEQE and PHREEQC with a transport code. The model CHMTRNS developed by Noorishad, et al.[1987] used PHREEQE as the equilibrium module. Engesgaard and Kipp [1992] used this model to validate their work. PHREEQE was further used by Narasimhan, et al.[1986] in their reactive transport model DYNAMIX, which couples this geochemistry model with the transport code TRUMP. An improved version of DYNAMIX model was used by Liu and Narasimhan [1989] to solve several redox-controlled reactive transport problems. Another popular reactive transport model PHT3D [Prommer, et al., 2003] uses MT3DMS [Zheng and Wang, 1999] for simulating the three-dimensional, advective-dispersive, multi-species transport and the geochemical model PHREEQC-2 [Parkhurst, 1999] for simulating 9 reaction processes. PHT3D uses PHREEQC-2 database files to define equilibrium and kinetic reactions. Another similar model has been developed by Phanikumar and McGuire [2004], integrating PHREEQC-2 into RT3D [Clement, 1997; Clement, et al., 1998]. The reactive transport models that use MINTEQA2 as equilibrium speciation module are MINTRAN [Walter, et al., 1994], Smith and Jaffe [1998] and OTEQ (One- dimensional Transport with EQuilibrium Chemistry). Walter, et al.[1994] developed MINTRAN by coupling the finite-element transport module PLUME2D with MINTEQA2. They tested their model with respect to ion-exchange chemistry, and precipitation-dissolution chemistry involving multiple sharp fronts. In their companion paper, they used the model to complete two-dimensional simulations of heavy metal transport in an acidic mine tailings environment. Smith and Jaffe [1998] modeled the transport and reaction of trace-metals in saturated soil and sediments using a one? dimensional reactive transport model. The transport processes included advection, and diffusive/dispersive mixing. They coupled the one-dimensional transport module to a modified version of MINTEQA2. The OTEQ model was developed by coupling a stream transport model OTIS to the equilibrium model MINTEQA2. This model has been used in surface water modeling [Runkel, et al., 1996; Runkel, et al., 1999]. Several models coupling transport with mixed equilibrium/kinetic reactions have appeared since the mid-1990s (Steefel and Lasaga [1994], McNab and Narasimhan [1994], Chilakapati [1995], Salvage, et al. [1996], Steefel and Yabusaki [1996], Abrams, et al.[1998], Salvage and Yeh [1998], Tebes-Stevens, et al. [1998], MIN3P [Mayer, 1999], Chilakapati, et al.[2000], Yeh [2001], Barry, et al.[2002], Mayer, et al [2002]). 10 Most of these models have limited capabilities for defining user-specified kinetic and equilibrium reactions. The reactive transport model RT3D [Clement, 1997; Clement, et al., 1998] is a user-friendly, MODFLOW-family reactive transport code. The code was originally developed for solving bioremediation problems and was later published as a general purpose, public-domain reactive transport model. It provides a reaction framework with a set of predefined modules to describe common bioremediation kinetics and it also offers the flexibility to add other complex reaction kinetics. Although RT3D is designed primarily for handling kinetic reactions, it can be designed to indirectly solve reactive transport problems involving equilibrium geochemical reactions. Li, et al. [2006] used RT3D to solve a mixed equilibrium-kinetics problem and simulated porosity reduction caused by mineral fouling in a permeable reactive barrier. Since RT3D is a widely used public domain code, enhancing the capability of RT3D to directly solve equilibrium problems will be of benefit to both academic and industrial research communities. 1.4 RESEARCH OBJECTIVES AND ORGANIZATION OF THE THESIS Multi-dimensional reactive transport models are needed for understanding the fate and transport of the metallic contaminants in groundwater aquifers. The goal of this research is to develop a framework for coupling MINTEQ family geochemical equilibrium models with the RT3D code. RT3D is a Fortran 90-based software package used for simulating three-dimensional, multi-species, reactive transport of chemical compounds in groundwater. It provides an easy-to-use framework with predefined modules to simulate common bioremediation kinetics; in addition it provides a flexibile option to add other 11 complex reaction kinetics involving multiple aqueous and adsorbed phase species. This flexible option can be used to integrate an existing geochemistry equilibrium module within RT3D. The objectives of this project are: (1) to couple MICROQL and MINEQL as reaction modules within a modular one-dimensional transport codes; (2) identify a set of benchmark problems involving processes including ion-exchange, precipitation- dissolution, surface-complexation and redox reactions, and solve them using the newly developed codes; (3) develop and validate a prototype model for incorporating MINTEQA2 into RT3D and test the model performance; and (4) develop a scalable reactive transport model to predict arsenic transport and speciation in a manganese dioxide containing porous medium. Chapter 1 provides the background literature and objectives of the current research effort. Chapter 2 provides a summary of four major geochemical processes that include (1) ion-exchange, (2) precipitation-dissolution, (3) surface-complexation (4) redox reactions. This chapter also provides the details of the mathematical techniques used for modeling these four processes. Chapter 3 presents the details of coupling one-dimensional multi-species transport code with MICROQL and MINEQL geochemical models. The MICROQL coupled transport code was then used to solve one-dimensional problems with ion-exchange and surface-complexation problems, and the MINEQL-coupled transport code was used for solving precipitation-dissolution and redox problems. 12 Chapter 4 presents the details of a prototype model that integrated the USEPA code MINTEQA2 within the multi-dimensional reactive transport code RT3D. The coupling framework developed in Chapter 3 was employed to develop this model. Efforts were also made to implement the new prototype model within the groundwater user interface GMS (Groundwater Modeling System). Chapter 5 of the thesis deals with the development of a scalable reactive transport model to predict arsenic transport in manganese dioxide containing soil columns. The ability of the numerical model to accurately predict Arsenic (As) transport was verified by comparing the model predictions against experimental results. This work was completed in collaboration with a PhD candidate , who performed the laboratory experiments. The author developed the scalable reactive-transport model, which was used to design the laboratory experiments. He also completed the model validation exercises. The focus of this chapter is to provide a summary of the modeling efforts. Chapter 6 provides the conclusions of this research effort and several recommendations for future work. 13 CHAPTER 2 MODELING VARIOUS TYPES OF EQUILIBRIUM REACTIONS 2.1 INTRODUCTION Major equilibrium reactions commonly encountered in groundwater systems are ion- exchange, surface-complexation, precipitation-dissolution and redox reactions. Numerical techniques for modeling these reactions are different. This chapter provides a mathematical description of these four primary reactions and reviews the numerical techniques for solving them. However, before describing the individual reactions, a general description of the problem of chemical equilibrium in aqueous systems is presented along with its solution strategy. 2.2 GENERAL DESCRIPTION OF THE PROBLEM OF CHEMICAL EQUILIBRIUM IN AQUEOUS SYSTEMS AND ITS SOLUTION The problem of chemical equilibrium in an aqueous system can be described in mathematical terms by a set of mass action and mass balance equations [Westall, 1979a]. Mass action equation can be written as: ( , ) 1 n a i j i i j j C K X = = ? for i=1,m (2.1) Mass balance equation can be written as: 1 ( , ) 0 m j i j i Y a i j C T = = ? =? for j=1,n (2.2) 14 Where iC and iK are the concentration and formation constant of ith species respectively, jX is the activity guess value of jth component and ( , )a i j is stoichiometric coefficient of jth component in ith species. The solution to the chemical equilibrium problem can be obtained by making an initial guess of Xj and then keep on updating them over several iterations till convergence is achieved. The component activities Xj are updated using Newton-Raphson method in the following manner [Westall, 1979a]: 1 1 ( ( , ) ( , ) / ) updated old m j jk i k ik X X X X Z Y YZ a i j a i k C X X X Change in X Z Jacobian of Y with respect to X ? ? ? = = ? ? ? = = = ? = = ? (2.3) 0.0, /10.0 updated updated old When X then X X < = (2.4) The solution process is summarized in Figure 2.1. 15 Figure 2.1 : MICROQL solution process for chemical equilibrium problem in aqueous systems. Initialize Input Complexes Convergence Check: Set Flag Mass Balance Invert Jacobians Too Many Iterations? Improved Values of X Converged? Output Yes No Yes Stop Stop No 16 2.3 MODELING ION EXCHANGE REACTION Ion exchange is a type of sorption reaction where a chemical species is replaced by another at the solid surface. Ion exchange equations explicitly account for all the ions which compete for the exchange sites. There are several conventions to relate activities to the concentration in the case of adsorbed cations. This will be shown here by the exchange of Na2+ for Ca2+. The exchange reaction can be written Gaines-Thomas convention as: 21 1 22 2Na CaX NaX Ca + ++ = + (2.5) The above equation can be written using Gapon convention as: 21 0.5 2Na Ca X NaX Ca + ++ = + (2.6) In MINTEQA2, the ion-exchange reactions are simulated using Gaines-Thomas convention. This Gains and Thomas model assumes that the surface sites are initially occupied by an exchangeable ion that is released during the exchange process. The other assumptions of the model are that the charge on the surface remains constant and the number of surface sites available for sorption, expressed as the cation exchange capacity (CEC) is fixed. The ion exchange reactions can be modeled as any other chemical equilibrium problem in aqueous systems using the method shown in Figure 2.1. This can be accomplished by considering surface site as any other aqueous component with fixed total concentration of CEC and other ion-exchanged species as aqueous species. The selectivity coefficients of the different ions are taken to be the equivalent to the equilibrium constants for the respective ion exchange reactions. 17 2.4 MODELING SURFACE-COMPLEXATION REACTIONS The adsorption on a surface is calculated as a function of the surface charge ? during the calculation of chemical equilibria using electrostatic models. The adsorption reactions are referred to as surface-complexation reactions. The computation of adsorption equilibria is described very well in Westall [1979b]. The charge on the surface is defined as the excess of positive groups over negative groups. In molar quantity, charge T ? is calculated as s aT F? ?= (2.7) where s is specific surface area, a is concentration of solid and F is Faraday?s constant. The four commonly used electrostatic models are following: (a) Constant Capacitance Model (b) Diffuse Layer Model (c) Stern Model (d) Triple Layer Model Surface is assumed to be at a potential ? with respect to the bulk of the solution. Several expressions have been proposed to relate ? to ?. These expressions differ between model to model. In these problems, the potential ? has a representative component exp(-F?/RT) and the equivalent total mass of this component is directly related to the surface charge ?. The solution process for chemical equilibrium problems, involving surface- complexation reactions are slightly different from those with only aqueous phase reactions. Unlike any other component, charge on the surface is not known experimentally. Starting with a guess value of potential ?, iterations are carried out till the calculated charge on the surface doesn?t change beyond the accepted limit. When the 18 equilibrium problem is solved, the sum of the charge on all the surface species must be equal to the electrostatically calculated charge T ? . The process of chemical equilibrium computation involves updating the guess activity over several iterations till the convergence is achieved. For adsorption equilibria, the certain elements of the Jacobian in equation (2.3) need to be modified in the following manner [Westall, 1979b]: 1 1 ( ( , ) ( , ) / ) updated old m i i X X X X Z Y Y TZ a i a i C X X X ? ? ?? ? ? ? ? ? ? ? ? = = ?? ? = ?= = ? ?? (2.8) Note that the total concentration of the electrostatic component T ? is not experimentally determined and is a function of the potential ?. The derivative of T ? with respect to X? is not zero and must be included in the normal expression for Jacobian given in equation (2.3). Figure 2.2 describes the MICROQL solution process for chemical equilibrium problem in aqueous systems, involving surface-complexation reactions. 19 Figure 2 . 2 : MICROQL solution process for chemical equilibrium problem in aqueous systems, involving surface-complexation reactions. Fundamental Constants, Adjustable Parameters, Problem Size, Identification of Surface Components Input Complexes Convergence Check: Set Flag Mass Balance Invert Jacobians Too Many Iterations? Improved Values of X Converged? Output Yes No Yes Stop Stop No Compute ?, ?, ?SOH Jacobian Modification of Jacobian 20 2.5 MODELING PRECIPITATION-DISSOLUTION REACTIONS The mineral precipitation-dissolution reactions are described by the mass action equation for the solid and the reacting ions. The general form of the dissolution reaction of a solid is given by ( ) ( ) ( )a bA B s a A aq b B aq? + (2.9) Where s refers to the solid phase and aq refers to the aqueous phase. The equilibrium constant for the reaction is called solubility product ksp and is defined as: { } { } { } a b a b A B spA B K= (2.10) The activity of solid is assumed to be unity, so the above expression can be written as: { } { }a b spA B K= (2.11) The term { } { }a bA B is called ion activity product (IAP). The tendency of a solid to precipitate or dissolve is judged by an index known as saturation index (SI), which is defined as: ( )lo g s pI A PKS I = (2.12) When SI > 0, the mineral tends to precipitate, when SI = 0, the mineral is in equilibrium with the solution, and when SI < 0, the mineral tends to dissolve or is not at all present. In MINTEQA2, each of the relevant solids is ranked for its tendency to precipitate by dividing the SI by the number of ions in the precipitation reactions. The solid with the highest SI value is allowed to precipitate and then all of the remaining relevant solids are ranked again and sequentially precipitated or dissolved, until all the solids except those precipitated are undersaturated. 21 During the chemical equilibrium calculation, MINTEQA2 classifies solids into four categories: (1) the solids which are in infinite supply, (2) the solids which are present in finite amount, (3) the dissolved solids and (4) the excluded solids that are not included in the simulation depending upon the individual problems. Solution process for chemical equilibrium problem in aqueous systems involving precipitation-dissolution reactions is presented in Figure 2.3. 22 Figure 2 . 3 : MINEQL solution process for chemical equilibrium problem in aqueous systems, involving precipitation-dissolution reactions. Initialize Input Ionic Strength Modify for Solids Print Input Data Solve for Soluble Species Unmodify for Solids Precipitate or Dissolve Output Solved Problem Yes No Stop 23 2.6 MODELING REDOX REACTIONS There are four different approaches available for modeling equilibrium-controlled redox processes. The description of these four approaches and their users are summarized below: 2.6.1 Effective internal approach This method was developed by Thortenson, one of the authors of PHREEQE [Parkhurst, et al., 1980]. It is based on the principle of conservation of electrons. If one redox reaction involves losing an electron through oxidation, there should exist another reaction that must be simultaneously gaining an electron by reduction. In this approach, a parameter called redox state is defined for the system. This redox state is sum of valence state (or redox state) of all the redox sensitive species (note the redox states of all the non-redox species are defined as zero. Walsh [1984], Liu and Narasimhan [1989] Engesgaard and Kipp [1992] and PHT3D [Barry, et al., 2002] used this approach. 2.6.2 External electron approach This approach considers a hypothetical electron activity as an aqueous component. For multivalent components, a species in the highest oxidation state is used as an aqueous component. All other lower valent species are described as aqueous complexes formed from the half-cell reaction between the aqueous component and the hypothetical electron. [Xu, 1996] and Liu and Narasimhan [1989] used this approach to model redox reactions. 2.6.3 Oxygen fugacity approach This approach specifies oxygen fugacity as the redox parameter in the system. Thus oxygen is included as an aqueous component and an additional mass balance equation is 24 solved for oxygen in the system of chemical speciation involving redox reactions. The geochemical model EQ3/EQ6 [Wolery, 1992] used this approach. 2.6.4 Redox couple approach The redox couple approach treats the redox couples as two separate aqueous components. The pE is calculated from the activities of the individual components of the redox couple after the completion of the equilibrium speciation. At equilibrium, all the redox couples show the same pE value. This approach has been used in EQ3/6 [Wolery, 1992] as well as PHREEQE [Parkhurst, et al., 1980]. 2.6.5 Selection of the redox approach for modeling reactive transport using MINEQL/MINTEQA2 In MINTEQA2, redox reactions can be represented using either the external electron approach or the redox couple approach. MINTEQA2 neither defines the parameter, redox state of the system nor uses oxygen to represent redox reactions, therefore the effective internal approach and the oxygen fugacity approach cannot be used. Our analysis indicates that the redox couple approach fails to converge to a single pE even for a simple system with two redox couples. Therefore in this work, the external electron approach was selected for modeling redox reactions. 2.6.6 Typical techniques for dealing with convergence problems Severe convergence problems were encountered while solving redox problem. Therefore, redox problems need special consideration. This is because electron is the only component with possible activity >1 (even up to 1010) and has negligible free ion concentration in solution. A special basis-switching technique along with a log-based formulation was used to alleviate the convergence problem. 25 2.6.6.1 Basis-switching technique Basis is the set of components in terms of which stoichiometric coefficients of the species participating in equilibrium reactions are written. Table 2.1 presents an example to demonstrate basis switching. This example is taken form the MICROQL manual [Westall, 1979a]. This example describes how the basis is switched from the set of components [Ca2+, CO32-, H+] to [Ca2+, CaCO3, H+]. This approach should be used when the activity of one of the components is negligible (in this example CO32 is assumed to be negligible and is replaced with CaCO3) compared to its total concentration. Basis- switching technique is proven to be the most effective way for dealing with convergence problems under such conditions. Table 2. 1 : An example of Basis-switching Ca2+ CO32- H+ Ca2+ 1 0 0 CO32- 0 1 0 H+ 0 0 1 CaOH+ 1 0 -1 CaHCO3+ 1 1 1 CaCO3 (aq) 1 1 0 H2 CO3 0 1 2 HCO3- 0 1 1 OH- 0 0 -1 Ca2+ CaCO3 H+ Ca2+ 1 0 0 CO32- -1 1 0 H+ 0 0 1 CaOH+ 1 0 -1 CaHCO3+ 0 1 1 CaCO3 (aq) 0 1 0 H2 CO3 -1 1 2 HCO3- -1 1 1 OH- 0 0 -1 26 2.6.6.2 Log-based formulation and its application In this approach, the activities of the components are used after transforming them to log space. The most crucial step in a solution process is the activity update step. We implemented the log-based formulation in the activity update step. Application of log- based formulation allowed us to solve a complex problem described in Liu and Narasimhan [1989] problem. The log-based formulation can be summarized as: 1 log log log log (2.303 ( , ) ( , ) )log log log log updated old j jk i ik X X X X Z Y YZ a i j a i k C X X Change in X Value Z Jacobian of Y with respect to X ? ? ? = ? ? ? = = = ? = = ? (2.13) Implementation of Log-based formulation has its own challenges. In course of a particular iteration, if | log |X? is large, the solution will fail to converge. The following empirical approach (suggested by Dr. Hammond at PNNL, personal communications) was utilized to deal with this: ,| log | 5.0 ( log 0.0) log 5.0 , log 5.0 In case X If X X Otherwise X ? > ? < ? = ? ? = (2.14) The above approximation appears to work well for solving redox problems. 27 2.7 SUMMARY A general description of the problem of chemical equilibrium in aqueous systems is presented along with its solution strategy. Differences in numerical techniques for individual reactions are illustrated through the use of flow charts. The major equilibrium processes included here are: ion-exchange, surface-complexation, precipitation- dissolution and redox reactions. Due to its inherent complexity and due to the availability of multiple solution approaches, special consideration was given to redox reactions. Out of four available approaches the external electron approach was found to be the most suitable approach for solving redox-controlled transport problems using MINTEQ-family codes. Numerical techniques used for alleviating the convergence problems related to redox reactions are also discussed. 28 CHAPTER 3 MODELING REACTIVE TRANSPORT BY COUPLING MICROQL/MINEQL TO ONE-DIMENSIONAL TRANSPORT CODES 3.1 INTRODUCTION In this chapter, one-dimensional reactive transport problems involving equilibrium geochemical reactions are considered. Geochemical reactions considered are ion- exchange, surface-complexation, precipitation-dissolution and redox reactions. A set of new computer codes are developed and their performance are validated by reproducing the previously reported results in the available literature. Simultaneously, attempts were also made to develop a common framework for modeling such reactions. The approach here involves selection of a particular type of equilibrium reaction, finding a frequently used benchmark problem for the reaction, devising a solution strategy for the problem and comparing the solutions with those reported in the literature. Individual section is dedicated to each of these four types of reactions. All the models couple a one-dimensional transport code to the equilibrium module MICROQL or MINEQL. A copy of each of these two codes developed in this study is provided in Appendices A and B of the supplementary material of this thesis. 29 3.2 REACTIVE TRANSPORT EQUATIONS AND THEIR GENERAL SOLUTION STRATEGY Reactive transport includes transport due to advection and dispersion and the geochemical reactions. 3.2.1. Equations representing the reactive transport Under one-dimensional steady state flow conditions, the transport of multiple components dissolved in ground water is governed by the following equations: mk R + xCv - xCD = tC aqkkkk ...,2,12 2 =?????? (3.1) where m is the total number of aqueous-phase (mobile) components, Ck is the aqueous- phase concentration of the kth component [ML-3], D is the hydrodynamic dispersion coefficient [L2T-1], v is the pore water velocity [LT-1], aqkR represents the rate of accumulation in aqueous phase concentration per bulk volume of the kth components [ML-3T-1] due to chemical reaction and x and t refer to distance [L] and time [T] respectively. Immobile phase concentration of the kth component due to precipitation pkC~ is given by: mk R = tC pk p k ...,2,1, ~ =???? (3.2) Where, pkR is the rate of accumulation in the form of precipitate aqueous phase concentration [ML-3T-1] due to chemical reaction and ?? and are soil bulk density 30 [ML-3] and porosity respectively. Similarly, immobile phase concentration of the kth component due to sorption or surface complexation, skC~ is given by: mk R = tC sk s k ...,2,1, ~ =???? (3.3) Where skR is the rate of accumulation in the form of sorbed aqueous phase concentration [ML-3T-1] due to chemical reaction. The reaction terms aqkR , pkR and skR are calculated by equilibrium calculation of the species present in the aqueous phase. 3.2.2 Solution of ADRE The transport equations are solved using fully implicit finite-difference method. Depending on the type of reaction, equilibrium model MICROQL or MINEQL is used to solve for the equilibrium concentrations. For ion-exchange and surface-complexation reactions, MICROQL is used and for precipitation-dissolution and redox reactions, MINEQL is used. The transport and the equilibrium modules are coupled by the sequential non-iterative procedure. Figure 3.1 presents the flow chart for generalized solution process for a reactive transport problem. 31 Figure 3.1 : Flow chart for generalized solution process of a reactive transport problem. Initial conditions and boundary conditions for all the components Set time step counter Solve transport equation for all the components Add all the sorbed/solids species concentrations to their respective components to calculate the total concentrations of the components Find the new aqueous conc. of all the components and update sorbed/solids conc concentrations Solve chemical equilibrium equations End of simulation Output Yes No Stop 32 3.3 REACTION 1: ION-EXCHANGE This test problem was taken from Valocchi et al. [1981]. In this paper, the investigators present a field experiment involving nonlinear ion exchange reactions in Palo Alto Baylands, California. In the field experiment, a treated municipal effluent was injected into a shallow aquifer and the breakthrough of the major cations and anions were monitored at an observation well 16m from the injection well. The principal chemical mechanism involved is the heterovalent ion exchange of Na+, Mg2+ and Ca2+. Other chemical reactions were ignored. The physical parameters used in this test problem are given in Table 3.1 and the initial and boundary conditions are provided in Table 3.2. Figures 3.2, 3.3 and 3.4 show the breakthrough profiles of Na+, Mg2+ and Ca2+ ions respectively. The results are in good agreement with those reported in Valocchi et al. (1981). 33 Table 3 . 1 : Physical parameters used in the simulation Parameter Value Pore water velocity 1.01 m/day Length of the domain 16 m Dispersivity 1.0 m Porosity 0.25 Table 3 . 2 : Initial and boundary conditions used in the simulation Aqueous Component Background Concentration (M) Injection Water Concentration (M) Na 8.70E-02 9.40E-3 Mg 1.80E-02 4.94E-4 Ca 1.10E-02 2.12E-3 Na-S 1.61E-01 0.0 Mg-S 1.41E-01 0.0 Ca-S 1.54E-01 0.0 34 Na+ Profile 1.00E+02 1.00E+03 1.00E+04 100 1000 10000 100000 Cubic Meter Injected Co nc . (m g/l ) Na+ (Present Simulation) Valocchi Figure 3 . 2: Breakthrough profile of Na+. Mg+2 Profile 1.00E+01 1.00E+02 1.00E+03 100 1000 10000 100000 Cubic Meter Injected Co nc (m g/l ) Figure 3 . 3: Breakthrough profile of Mg2+. 35 Ca+2 Profile 1.00E+01 1.00E+02 1.00E+03 100 1000 10000 100000 Cubic Meter Injected Co nc (m g/l ) Present Simulation Valocchi Figure 3 . 4 : Breakthrough profile of Ca2+. 3.4. REACTION 2: SURFACE-COMPLEXATION This problem was taken from Cederberg et al. [1985]. This problem involves reactive transport of cadmium, chloride and bromide in a one-dimensional soil column. The reaction involves formation of aqueous and surface complexes and sorption of free cadmium to the soil matrix. The sorption process was described by constant capacitance model. The physical and chemical parameters used in the simulation are given in Tables 3.3 and 3.4. Table 3.5 provides the initial and the boundary conditions used in the 36 simulations. Table 3.6 provides the stoichiometry coefficients of the reactions occurring in the aqueous and the sorbed phases. Figure 3.5 shows the profiles of Cd(Aq), Cd (sorb) and Cl obtained at t=15hrs during simulation of CdAq transport for case1. Results are in agreement with those reported in the above paper. Figure 3.6 presents the results for the cases 2 and 5 and Figure 3.7 presents the result for cases 3 and 6. Table 3 . 3 : Physical parameters used in the simulation Parameter Value Pore water velocity 8 cm/day Length of the domain 10 cm Longitudinal dispersion coeff. 0.16 cm2/day Porosity 0.30 Bulk density of medium 2500 g/l 37 Table 3 . 4 : Chemical sorption parameters used in the simulation Parameter Value Total no. of sites 0.046 mol/l Ionic strength 0.1 mol/l pH (held constant) 7.0 Capacitance 1.06 F/m2 Solid surface area 1.0 m2/g Table 3 . 5 : Concentrations of components for the test cases: Initial and Boundary conditions Background Concentration (M) Injection Water Concentration (M) Cases CdT ClT BrT CdT ClT BrT Case 1 0.00001 0.003 0.0001 0.0001 0.003 0.0001 Case 2 0.0001 0.0003 0.0001 0.0001 0.003 0.001 Case 3 0.00001 0.0003 0.0001 0.0001 0.003 0.001 Case 4 0.00001 0.003 0.001 0.0001 0.003 0.001 Case 5 0.0001 0.003 0.001 0.0001 0.03 0.01 Case 6 0.00001 0.003 0.001 0.0001 0.03 0.01 38 *This should be +1; however Cederberg?s work used this incorrect value. To be consistent with their work and to reproduce the results, we also used -1 in our simulations. Table 3 . 6: Stoichiometry of Reactions Species Cd2+ Cl- Br- SOH Psi H+ LOG K 1 H+ 0.0 0.0 0.0 0.0 0.0 1.0 0.000 2 Cd2+ 1.0 0.0 0.0 0.0 0.0 0.0 0.000 3 Cl- 0.0 1.0 0.0 0.0 0.0 0.0 0.000 4 Br- 0.0 0.0 1.0 0.0 0.0 0.0 0.000 5 CdCl+ 1.0 1.0 0.0 0.0 0.0 0.0 1.800 6 CdCl2 1.0 2.0 0.0 0.0 0.0 0.0 2.600 7 CdBr+ 1.0 0.0 1.0 0.0 0.0 0.0 2.200 8 CdBr2 1.0 0.0 2.0 0.0 0.0 0.0 3.000 9 CdOH+ 1.0 0.0 0.0 0.0 0.0 -1.0 -12.690 10 OH- 0.0 0.0 0.0 0.0 0.0 -1.0 -13.910 11 SOH 0.0 0.0 0.0 1.0 0.0 0.0 0.000 12 SOH2+ 0.0 0.0 0.0 1.0 1.0 1.0 7.400 13 SO- 0.0 0.0 0.0 1.0 -1.0 -1.0 -9.240 14 SOCd+ 1.0 0.0 0.0 1.0 -1.0* -1.0 -7.000 39 Aqueous Cd Transport (Case I and IV) 1.00E-06 1.00E-05 1.00E-04 1.00E-03 1.00E-02 1.00E-01 0 2 4 6 8 10 X (cm) Co nc en tra tio n ( M) Cl_T (Case I) Cl_T (TRANQL) Cd_Aq (Case I) Cd_Aq (TRANQL) Cd_Sorbed (Case I) Cd_Sorbed (TRANQL) Figure 3 . 5 : Profiles obtained at t=15hrs during simulation for case I . Aqueous Cd Transport (Case II and V) 1.00E-06 1.00E-05 1.00E-04 1.00E-03 1.00E-02 1.00E-01 0 2 4 6 8 10 X (cm) Co nc en tra tio n ( M) Cl_T (Case II) Cd_Aq (Case II) Cd_Sorbed (Case II) Cl_T (Case V) Cd_Aq (Case V) Cd_Sorbed (case V) Figure 3. 6 : Profiles obtained at t=15hrs during simulation for cases 2 and 5. 40 Aqueous Cd Transport (Case III and VI) 1.00E-06 1.00E-05 1.00E-04 1.00E-03 1.00E-02 1.00E-01 0 2 4 6 8 10 X (cm) Co nc en tra tio n ( M) Cl_T (Case III) Cd_Aq (Case III) Cd_Sorbed (Case III) Cl_T (Case VI) Cd_Aq (Case VI) Cd_Sorbed (Case VI) Figure 3. 7 : Profiles obtained at t=15hrs during simulation for cases 3 and 6. 3.5. REACTION 3: PRECIPITATION-DISSOLUTION This problem has been taken from Engesgaard and Kipp [1992]. This system involves two minerals and five complexation reactions along with transport of a conservative tracer (Chloride). The parameters used in the simulation are given in Table 3.7 and the initial and boundary conditions are provided in Table 3.8. Initially, the calcite mineral has uniform presence in the soil column. As the incoming fluid does not contain Ca2+ and CO32-, the mineral becomes undersaturated and gets dissolved. The aqueous phase Ca2+ and CO32- obtained from the dissolution of calcite, reacts with Mg2+ bearing incoming fluid and precipitation of dolomite (CaMg(CO3)2) is initiated. Figures 3.8, 3.9, 3.10 and 3.11 show the concentration profiles 41 of aqueous and solid species at t=21000 sec. The results are completely in agreement with those reported in Engesgaard and Kipp [1992]. Table 3 . 7 : Physical parameters used in the simulation Parameter Value Pore water velocity 9.37E-6 m/s Length of the domain 0.5 m Dispersivity 0.0067 m Porosity 0.32 Bulk density of medium 1800 Kg/m3 Table 3 . 8 : Initial and boundary conditions used in the simulation Parameter Background Concentration (M) Injection Water Concentration (M) pH 9.91 7.06 Ca2+ 1.239E-4 0.0 CO32- 1.239E-4 0.0 Mg2+ 0.0 1.0E-3 Cl- 0.0 2.0E-3 CaCO3 (s) 1.2206E-4 0.0 CaMg(CO3)2 (s) 0.0 0.0 42 0.00E+00 5.00E-06 1.00E-05 1.50E-05 2.00E-05 2.50E-05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Distance(m) Mo le/ Kg of S oil Calcite(Present Simulation) Calcite(Kipps and Engesgaard) Dolomite (Present Simulation) Dolomite(Kipps and Engesgaard) Figure 3 . 8 : Calcite and Dolomite profiles at 21,000 sec. 43 7 7.5 8 8.5 9 9.5 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Distance (m) St an da rd U nit s pH (Present simulation) pH (Kipp and Engesgaard) Figure 3 . 9 : pH profile at 21,000 sec. 0.00E+00 1.00E-05 2.00E-05 3.00E-05 4.00E-05 5.00E-05 6.00E-05 7.00E-05 8.00E-05 9.00E-05 1.00E-04 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Distance(m) Co nc (M ol/ l) Carbonate (Present Simulation) Carbonate (Kipp and Engesgaard) Bicarbonate (Present Simulation) Bicarbonate (Kipp and Engesgaard) MgCo3(Present Simulation) MgCO3 (Kipp and Engesgaard) CaCo3(Present Simulation) CaCO3 (Kipp and Engesgaard) Figure 3 . 10 : Aqueous concentration profiles at 21,000 sec. 44 0.00E+00 2.00E-04 4.00E-04 6.00E-04 8.00E-04 1.00E-03 1.20E-03 1.40E-03 1.60E-03 1.80E-03 2.00E-03 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Distance (m) Co n ( M) Ca2+ (Present Simulation) Ca2+ (Kipp and Engesgaard) Mg2+(Present Simulation) Mg2+ (Kipp and Engesgaard) Cl-(Present Simulation) Cl- (Kipp and Engesgaard) Figure 3 . 11 : Aqueous concentration profiles at 21,000 sec. 3.6. REACTION 4: REDOX This problem has been taken from Liu and Narasimhan [1989] and it was also solved by Xu[1996]. This system involves redox reactions accompanied by precipitation and dissolution. A hypothetical one-dimensional column of length 3 m is used. The parameters used in the simulation are given in Table 3.9 and the initial and boundary conditions are provided in Table 3.10. 45 Table 3 . 9 : Physical parameters used in the simulation Parameter Value Pore water velocity 1.0E-6 m/s Length of the domain 3.0 m Longitudinal dispersivity 0.02 m Molecular diffusion coeff. 1.0E-10 m2/s Table 3 . 10 : Initial and boundary conditions used in the simulation Parameter Background Concentration (Mol/l) Injection Water Concentration (Mol/l) pH 9.0 6.9 pE -5.912 5.07 CO32- 2.4E-4 3.0E-3 Si(OH)4 1.1E-4 1.1E-4 Na+ 1.3E-2 1.0E-3 Ca2+ 1.0E-3 6.7E-4 UO22+ 0.0 4.315E-5 USiO4 (Coffinite) 0.0 0.0 CaUO4 (Calcium Uranate) 0.0 0.0 In this problem, initially the soil column was assumed to be free of reactive solid phases and uranium species, and the aqueous solution was under alkaline reducing condition. At time greater than zero, the uranium bearing incoming fluid enters the system under acidic 46 and oxidizing conditions, coffinite precipitates from the aqueous solution. At the later stage, as the aqueous solution attains moderately oxidizing condition, precipitation of calcium uranate was initiated along with the concomitant dissolution of coffinite. Figures 3.12, 3.13 and 3.14 present the pH, pE and minerals? profiles at 2.0E6 seconds. The results predicted by our model shows excellent agreement with those predicted using TRANQUI [Xu, 1996]. 6 6.5 7 7.5 8 8.5 9 9.5 10 0 0.5 1 1.5 2 2.5 3 Distance(m) pH pH Profile (Present Simulation) pH Profile(Xu) Figure 3 . 12 : pH profile at 2.0E6 sec. 47 -7 -5 -3 -1 1 3 5 7 0 0.5 1 1.5 2 2.5 3 Distance(m) pE pE Profile (Present Simulation)pE (Xu) Figure 3 . 13 : pE profile at 2.0E6 sec. Mineral Profiles 0.00E+00 1.00E-04 2.00E-04 3.00E-04 4.00E-04 5.00E-04 6.00E-04 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Distance (m) So lid E qu iva len t ( M) Calcium Uranate Coffinite Calcium Uranate(Xu) Coffinate (Xu) Figure 3 . 14 : Mineral Profiles at 2.0E6 sec. 48 3.7 SUMMARY In this chapter the performance of the numerical model was validated by solving four standard benchmark problems, identified from the published literature. The model was able to successfully simulate four types of processes including ion-exchange, surface- complexation, precipitation-dissolution, and redox reactions. The table below provides a summary of these model development efforts. Table 3. 11 : Summary of the model development work in Chapter 3 Reaction Type Research article used Equilibrium Model Comments Ion-exchange Valocchi et al.[1981] MICROQL- I Results in agreement with those in literature Surface- complexation Cederberg, et al.[1985] MICROQL- II Results in agreement with those in literature Precipitation- dissolution Engesgaard and Kipp [1992] MINEQL Results in agreement with those in literature Redox Liu and Narasimhan [1989]& Xu [1996] MINEQL Results in agreement with those in Xu [1996] 49 CHAPTER 4 MODELING REACTIVE TRANSPORT BY COUPLING MINTEQA2 TO RT3D 4.1 INTRODUCTION Modeling of transport mediated by geochemical reactions requires coupling of a transport code with an equilibrium geochemistry code. There are two widely used reaction modules, PHREEQC (USGS software) and MINTEQA2 (USEPA software), currently available in the public domain for solving equilibrium geochemistry problems. This chapter summarizes the implementation details for coupling the geochemistry module MINTEQA2 with the general three-dimensional reactive transport model RT3D modeling framework, the model development and its validation. MINTEQA2 computer code is an advanced version of the basic code MICROQL. The computer codes MICROQL, MINEQL and MINTQA2 employ almost identical solution procedure, but can model geochemical reactions with various levels of complexity. MINTEQA2 is the most comprehensive code, but it is also the most computationally demanding source code. To exploit the computational advantages of small source codes, it was decided to couple and test a one-dimensional multi-species code with the MICROQL to solve problems involving equilibrium reactions and surface-complexation reactions, and with the MINEQL code for solving problems involving precipitation-dissolution reactions. The details of this task were provided in chapter 3. These exercises helped us probe into 50 the individual routines of MICROQL and MINEQL. The understanding gained from these research tasks was utilized to develop a prototype-model that coupled the full version of MINTEQA2 with RT3D. Further efforts were made to integrate this prototype model into the groundwater user interface GMS (Groundwater Modeling System). 4.2 THREE-DIMENSIONAL REACTIVE TRANSPORT EQUATIONS AND THE SOLUTION STRATEGY 4.2.1 Reactive Transport Equations The reactive transport system considered includes transport due to advection and dispersion and the geochemical reactions. The general macroscopic equations that describe the fate and transport of multiple aqueous- and solid-phase reactive species in three dimensions are written as [Clement, 1997; Clement, et al., 1998]: ( ) kk k k s cqC C s = - C + C r , where, k 1,2,...mvDij it x x x i j i ? ?? ?? ?? ? + =? ?? ? ? ? ? ? ? (4.1) ( )im cC = , where, im = 1,2,..., n mrt? ??tildenosp tildenosp (4.2) where n is the total number of species, m is the total number of aqueous-phase (mobile) species (thus, n minus m is the total number of solid-phase or immobile species), Ck is the aqueous-phase concentration of the kth species [M L-3], imC~ is the solid-phase (sorbed or precipitated) concentration of the imth species [either M M-1 (contaminant mass per unit mass of porous media) or M L-3 (contaminant mass per unit aqueous-phase volume) unit basis can be used], Dij is the hydrodynamic dispersion coefficient [L2 T-1], ? is the soil porosity, qs is the volumetric flux of water per unit volume of aquifer representing sources and sinks [T-1], Cs is the concentration of source/sink [M L-3], rc represents the 51 rate of all reactions that occur in the aqueous phase [M L-3 T-1], cr~ represents the rate reactions occurring in the soil-phase, and v is the saturated groundwater flow velocity in the pores [L T-1]. The flow velocities are calculated from the hydraulic-head values computed using a flow model. 4.2.2 Integration of Transport and Geochemistry Routines The first step in the solution process involves solving the advection and dispersion steps for all the mobile species. The sequential non-iterative operator-split approach is used to split transport and reaction steps. All the species controlled by equilibrium reactions are grouped into a set of components. The component concentrations were selected as the independent variables in the equilibrium problem. Note that any component can exist either in the aqueous phase and/or in the solid/sorbed phase. Therefore, total component concentration Tk is: 1, ...........,= + =k k k cT C S k N (4.3) Where Nc is the number of components, kC is the aqueous phase portion and Sk is the solid-phase portion of the components, which are defined using the following equations: 1 1, . .. .. .. . , = = =? aN k l k l c l C A c k N (4.4) 1 1,..........., = = =? sN k lk l c l S B s k N (4.5) Where Na is the number of aqueous phase species, Ns is the number of solid phase species, cl is the concentrations of the aqueous phase species, sl is the concentrations of solid phase (sorbed or precipitated) species. lkA is the stoichiometric coefficient of 52 component k in the aqueous species l, lkB is the stoichiometric coefficient of component k in the sorbed/solid species l. Na and Ns are the number of aqueous and solid phase species, respectively and Nc is the number of components. The transport and reaction modules are coupled using the sequential non-iterative approach. Conceptually, the coupling process can be visualized as a physical step where the aqueous components are advected and dispersed through space without reaction, and a chemical step where the components react without undergoing any transport. The transport equations are solved for all the components to obtain their respective aqueous phase concentrations transkC at the time step (t+?t). The incoming fluid are be equilibrated with ambient solid phase concentration using MINTEQA2 routines. This process is repeated until the required time level. This reaction coupling approach is similar to the method used by Phanikumar and McGuire [2004]. 4.3 COUPLING OF RT3D AND MINTEQA2 ROUTINES The current EPA version of MINTEQA2 required several modifications to allow coupling with the RT3D software. It is important to note that the original MINTEQA2 can only analyze batch systems; this project effort allowed to fully integrate this batch model within a user-friendly reactive transport code RT3D. The standalone version of MINTEQA2 required two I/O files. The first file is the input file read by MINTEQA2, and the second file is the output file generated by MINTEQA2 after solving the batch chemistry. The input file allows the user to describe the geochemical problem by providing the data related to the components and species for the system of interest. The output provides the final equilibrated chemistry. In order to couple the reaction module 53 of MINTEQA2 with transport modules of RT3D, the I/O information should be seamlessly transferred between the codes; this can be accomplished using two options. The first option is generate MINTEQA2 input files from RT3D, and then run RT3D; later from RT3D read the required information from the MINTEQA2 output files. This option is relatively easy to code but our test runs indicated that this is a computationally inefficient option. The more robust option is the second option where all the I/O operations are internally transferred between RT3D to MINTEQA2 codes; our test results show that this option runs much faster than the first option. Extensive modifications were made to MINTEQA2 I/O routines to facilitate this internal transfer. The code used for coupling MINTEQA2 to RT3D is listed in Appendix C of the supplementary material of this thesis. In addition to the above improvement, modifications were also made to MINTEQA2 routines that interacted with the reaction database. When run in the original batch mode, MINTEQA2 automatically selects all the possible types of species (aqueous, solid, sorbed) from the database corresponding to the input components. However, when MINTEQA2 is used as a reaction module within a reactive transport code, this step is not required except for the first time. As a part of the code development efforts, MINTEQA2 was modified in such a way that it generates the relevant species only at the initial time step. The species information is stored and reused when MINTEQA2 is called during subsequent time steps; this avoids the search of database at every time step. This further helped to reduce the computational time by almost 50%. The current version of RT3D deals with two types of species/components: mobile (aqueous) and immobile (sorbed). The equilibrium reaction module MINTEQA2 requires 54 the transport of the total concentrations of all the components. To integrate MINTEQA2 into existing RT3D framework, all the MINTEQA2 components are treated as mobile components. Similarly, all the solid and sorbed species in equilibrium chemistry are treated as immobile components. The I/O operations for entering the initial and boundary values of all the components are fully integrated within the GMS software (via the user- defined reaction module). The figure below shows how the component information can be directly input using the GMS user interface. Figure 4 . 1 : GMS user interface for the new geochemistry module. 55 4.4 EXAMPLE PROBLEM The surface complexation problem and the precipitation-dissolution problems described in the previous section were solved again using the new RT3D module. In addition, another new test problem that involves multiple precipitation-dissolution and redox reactions was also solved. This test problem was originally proposed by Walter et al. (1994). The problem considers geochemical transport processes occurring in a carbonate aquifer due to the presence of acid mine tailings effluent. The chemical conditions at the site are characterized by distinct zones where the water chemistry is in equilibrium with respect to a series of mineral pH buffers. The simulation domain is a one-dimensional column of aquifer material. Information related to simulation domain and other physical parameters are presented in Table 4.1. Fourteen aqueous components and six mineral phases are considered in the simulation. Table 4.2 provides the initial and boundary conditions for all the components and minerals. In this test problem, as the acidic tailings water entered into the column, calcite dissolved and buffers the pH. As long as there is sufficient calcite in the system, the pH remained relatively constant. After exhausting the calcite minerals, the second pH-buffering mineral siderite started to precipitate within the system. At this stage, the pH was controlled by equilibrium with respect to siderite. Siderite now dissolved and the pH remained relatively constant until siderite was exhausted. After the depletion of siderite, the pH decreased once again to a new level controlled by the equilibrium chemistry of gibbsite. The precipitation-dissolution of these three minerals and the resulting changes in pH, fully controlled the changes in the 56 aqueous phase concentrations of other chemical components. The results after 12days of transport are shown in Figures 4.2 to 4.6. Table 4 . 1 : Physical parameters used in the simulation Parameter Value Pore water velocity 2.0 cm/day Length of the domain 40 cm Dispersivity 0.5 cm Porosity 0.35 Grid spacing 0.5 cm Table 4 . 2 : Initial and boundary conditions used in the simulation Components/Minerals Initial Condition (M) Boundary Condition (M) Ca 6.92e-3 1.08e-2 Mg 1.96e-3 9.69e-4 Na 1.30e-3 1.39e-3 K 6.65e-5 7.93e-4 Cl 1.03e-3 1.19e-4 CO3 3.94e-3 4.92e-4 SO4 7.48e-3 5.00e-2 Mn 4.73e-5 9.83e-6 H4SiO4 1.94e-3 2.08e-3 Fe(II) 5.39e-5 3.06e-2 Fe(III) 2.32e-8 1.99e-7 Al 1.27e-7 4.30e-3 pH 6.96 3.99 pE 1.67 7.69 Calcite 1.95e-2 0.0 Siderite 4.22e-3 0.0 Amorphous Silica 4.07e-1 0.0 57 Gibbsite 2.51e-3 0.0 Ferrihydrite 1.86e-3 0.0 Gypsum 0.0 0.0 CaCO3 0.00E+00 5.00E-03 1.00E-02 1.50E-02 2.00E-02 2.50E-02 3.00E-02 0 5 10 15 20 25 30 35 40 Distance (cm) mo l/l Present Simulation Walter Figure 4 . 2 : Concentration profile of calcite mineral along the column. 58 Al(OH)3 0 0.005 0.01 0.015 0.02 0.025 0.03 0 5 10 15 20 25 30 35 40 Distance (cm) mo l/l Present Simulation Walter Figure 4 . 3 : Concentration profile of gibbsite mineral along the column. pH 4 4.5 5 5.5 6 6.5 7 0 5 10 15 20 25 30 35 40 Distance (cm) pH Present Simulation Walter Figure 4 . 4 : Predicted pH profile along the column. 59 Ca 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0 5 10 15 20 25 30 35 40 Distance (cm) mo l/l Present Simulation Walter Figure 4 . 5 : Concentration profile of Ca2+ (aqueous) along the column. SO4 0 0.01 0.02 0.03 0.04 0.05 0.06 0 5 10 15 20 25 30 35 40 Distance (cm) mo l/l Present Simulation Walter Figure 4 . 6 : Concentration profile of SO42- (aqueous) along the column. 60 4.5 Summary The current version of the public-domain reactive transport model RT3D does not have the capability to model equilibrium-controlled reactions. This chapter summarizes an attempt to integrate the USEPA equilibrium speciation model MINTEQA2 within RT3D. The goal was to enhance the capability of RT3D to solve equilibrium-controlled reactive transport problems. The developed model is tested by solving three one-dimensional benchmarks problems. The model results show a good match with the published results. However, further efforts are needed to test its capability for solving more complex problems in multiple dimensions. 61 CHAPTER FIVE DEVELOPMENT OF A SCALABLE NUMERICAL MODEL FOR PREDICTING COLUMN-SCALE ARSENIC TRANSPORT USING BATCH DATA 5.1 INTRODUCTION This chapter summarizes the development of a scalable reactive transport model for predicting arsenic (As) transport in manganese dioxide containing soil columns. A kinetic model coupled with a one-dimensional transport code was developed to analyze this problem. The model capability was verified by comparing model predictions against experimental results. This work was completed in collaboration with a PhD candidate [Radu, 2006] who performed the laboratory experiments. The author developed the scalable reactive-transport model which was used to design the laboratory experiments. He also completed the model validation efforts. The focus of this chapter is to provide a summary of these modeling efforts. 5.2 BACKGROUND Arsenic (As) is known for its acute toxicity to humans. Due to its chronic toxicity effects, As maximum contaminant level has recently been reduced from 50 ppb to 10 ppb. Arsenic contamination of groundwater has been of great interest in several places, most notably in Bangladesh and West Bengal, where two-thirds of the population is 62 at risk of serious health effects due to high concentrations of arsenic in the shallow alluvial and deltaic aquifers that are the main sources of drinking . As(III) and As(V) are most predominant forms of As in groundwater. As(V) is less toxic than As(III) and adsorbs more strongly to the soils compared to As(III). Therefore the oxidation of As(III) to As(V) and As(V) adsorption by various minerals present in subsurface is crucial in dealing with the risks posed by As contamination. Due to their excellent oxidative and adsorptive capabilities, naturally present manganese oxides can play a significant role in As oxidation and adsorption. Pyrolusite, a manganese oxide mineral present in abundance in natural systems was chosen in this study. The objective of this research is to develop a reactive transport model for predicting the effects of pyrolusite on As transport in one- dimensional soil columns. A scalable transport model was developed to predict column- scale transport by scaling the parameters obtained from the batch experiments. Several column experiments were designed to test the predictive capability of the scalable reactive transport model. 5.3 RESULTS FROM THE BATCH EXPERIMENTS Detailed discussion of the batch experiments are given in Radu [2006]. This section provides a summary of the batch data relevant to the modeling study. 5.3.1 Results from the oxidation kinetics and adsorption experiments The As(III) oxidation rate was found to be linearly dependent on both its concentration and the MnO2(s) concentration, represented by the following rate expression [Radu, 2006]: 63 )]()][([)]([ 2 sMnOIIIAsKdtIIIAsd ?= (5.1) Where K is As(III) oxidation rate constant [M-1L3T-1]. K value obtained from the batch experiments was 0.0012 lg-1min-1. Adsorption of As(V) to pyrolusite was found to be an almost instantaneous process since the entire uptake took less than a minute [Radu, 2006]. Langmuir adsorption isotherm model was found to best represent our experimental data [Radu, 2006]. In the Langmuir adsorption model, the adsorbed concentration is related to the aqueous concentration by the following equation: CK CqKQ l l += 1 max (5.2) Where Kl and qmax are Langmuir adsorption parameters (values given in Table 5.1) and C is aqueous concentration of As(V). It was determined that MnO2(s) has a maximum adsorptive capacity of 4 ?g As/g of MnO2(s) [Radu, 2006]. Table 5 . 1 : Model Parameters (Radu, 2006) Parameter Value Inlet As conc (?M) 13 MnO2 porosity 0.54 Silica sand porosity 0.425 Column height (cm) 7 Column diameter (cm) 1 Pore water velocity (cm/min) 0.79 Peclet No. 0.07 kl (cm3/mol) 2.0 qmax(?g As/g of MnO2(s)) 4 64 5.4 Reactive transport modeling and designing the column experiments The batch experiments provided the necessary building blocks to develop a numerical model for predicting the behavior of As species in soil systems. The modeling efforts include development of a conceptual model derived from scaling the batch results and development of a numerical model. The developed model was then used to design a set of column experiments and the experimental observations were then compared with the model predictions to check the efficacy of the reactive transport model. 5.4.1 Conceptual model. Several key observations made during the batch study were utilized for developing the conceptual model for MnO2-As system. Based on these observations, following assumptions were made in the conceptual model : (i) As(III) is the only species that undergoes oxidation; (ii) Only As(V) adsorption was considered in the model; (iii) As(V) adsorption to the MnO2(s) surface follows Langmuir adsorption isotherm; (iv) MnO2(s) is uniformly distributed throughout the column. The transport of As(III) can be described by following processes: 1. As(III) enters the column and gets partially or fully oxidized to As(V), depending on the oxidation rate; 2. Adsorption of As(V) onto the MnO2(s) grains is assumed to be strong and is controlled by the Langmuir adsorption isotherm; 3. The non-adsorbed As(V) fraction gets transported further by the pore water. In the case of As(V) transport, the transport process can be described by just the last two steps. 65 5.4.2 Reactive transport model. The equation governing the transport and the reactions of As (III) is: )]([2 2 )]([)]([)]([ IIIAscolumnOx IIIAs x IIIAsD t IIIAs ? ? ?? ? ?= ? ? ? (5.3) Where aqueous-phase concentration of As(III) is given in units of [ML-3], D is the hydrodynamic dispersion coefficient [L2T-1], v is the pore water velocity [LT-1], and Ocolumn[As(III)] represents the rate of As(III) oxidation [T-1] in the column. The equation governing the transport of As(V) is: )()]([)]([)()]([ 2 2 2 IIIAsO x VAs x VAsD t qf t VAs column MnO + ? ?? ? ?= ? ?+ ? ? ? ? ? (5.4) where aqueous-phase concentration of As(V) is given in units of [ML-3], ?MnO2 is the bulk density of MnO2(s) in the column, q is the adsorbed phase concentration of As(V) related to the aqueous phase concentration by the equation (5.2), f is the fraction of pyrolusite and ? is the porosity. Using the adsorption isotherm expression, equation (5.4) can also be written as: ( ) )( )]([)]([ )]([11 )]([ 2 2 2 max2 IIIAsO x VAs x VAsD VAsK qKf t VAs column l lMnO + ? ?? ? ?= ?? ??? ?? ??? ??? ? ??? ? ++? ? ? ? ? (5.5) 5.4.3 Solution of the reactive transport equations A fully-implicit finite-difference method is used to discretize the advection-diffusion and reaction equations (5.3) and (5.5). Note the equation for As(V) transport involves nonlinear terms in As(V). The Picard iterative method is used to deal with this nonlinear term. The resultant set of linear equations for both As(III) and As(V) are solved using 66 tridiagonal matrix solver routine. The model parameters used in the simulations have been presented in Table 5.1. )]([)]([)]([)]([ 2 2 IIIAsOxIIIAsx IIIAsDtIIIAs column??????=?? ? (5.6) Equation (5.6) is discretized as 1 1 1 1 2 1 1 1 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) n n n n n j j j j j n n nj j column j t x v x As III As III As III As III As III As III As III As IIIO ? ? + = ? ? ? + + + + + + + ? ? ? (5.7) Defining Courant No. = v tx?? and Diffusion No. = 2D tx?? , the above equation becomes 1 1 1 11( ) ( ) ( ) n n n j jjA B C DAs III As III As III + + + +? + + = (5.8) Where A, B and C are coefficients represented by A = -Courant No. /2.0- Diffusion No. B = 1.0 + 2 * Diffusion No. + Ocolumn * AS(III) * ?t (5.10) C = Courant No. /2.0 ? Diffusion No. D = ( )njAs III ( ) )]([ )]([)]([ )]([11 )]([ 2 2 2 max2 IIIAsO x VAs x VAsD VAsK qKf t VAs column l lMnO + ? ?? ? ?= ?? ??? ?? ??? ??? ? ??? ? ++? ? ? ? ? (5.11) Equation (5.11) is discretized as 2 1 1 1 1 1 1max 2 2 1 1 11 1 2 1 (1 [ ( )] ) 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) n n n n n M nOj j j j jl l iter n n nj j column j k q t k As V x v x fAs V As V As V As V As V As V As V As IIIO ? ?? ? + = ? ? ? + + + + + ? + + ++ ? ? ?? ?? ? + ? ?? ?? ?? + ? ? ?? ?? ? ? (5.12) 67 Defining Courant No. and Diffusion No. as earlier and Sorpterm as 2 m a x 21 (1 [ ( )] ) M n O l l ite r k q k A s V f ? ?? ?? ?? ? + ? ?? ?? ?+ ? ?? ?? ? , the above equation becomes 1 1 1 1 1( ) ( ) ( ) n n n j j jA B C DA s V A s V A s V + + + ? ++ + = (5.13) Where A, B and C are coefficients represented by A = - Courant No. /2.0- Diffusion No. B=1.0+2* Diffusion No.+ Sorpterm (5.14) C= Courant No. /2.0- Diffusion No. D= ( )njAs V *(1+Sorpterm) + Ocolumn * delt * ( )njAs III 5.4.4 Moving from batch to column scale simulations In order to simulate the column scale transport, we first considered the factors which affect adsorption and oxidation processes, while moving from the batch to the column scale. The key factor is the amount of MnO2(s) in the system. The batch experiments showed that the parameters oxidation rate constant and adsorption capacity are linearly proportional to MnO2 (s) concentration. Therefore, the oxidation rate constant can be scaled from the batch to column as follows: Ocolumn =Obatch * (Mass of MnO2 (s)/ Nx) / Pore-volume (5.14) Where Ocolumn is As(III) oxidation rate constant, which is going to be used in predicting the results of column experiments. Obatch is the oxidation-rate constant obtained from the batch experiments. Mass of MnO2(s) (in grams) is its amount used in the entire column and Nx is the number of nodes in the column. Pore-volume is calculated as follows: Pore-volume = Cross-sectional area *Column length * Porosity (5.15) 68 Since adsorptive capacity doesn?t depend on the solid to solution ratio, its value need not be adjusted for change in solid to solution ratio. Simulations are run for several hypothetical experimental conditions and the following experimental conditions were selected to be tested in the lab: (a) Entire column packed with MnO2(s): To achieve highest oxidation and adsorption of As, entire column was packed with MnO2(s). To study the difference in As(III) and As(V) transport, two separate experiments were conducted with As(III) and As(V). The experiment with input solution containing As(III) only was designated as Exp A and that with As(V) only was designated as Exp B. (b) Column packed with sand-MnO2(s) mixture with MnO2(s) 50% by wt.: The purpose of choosing this experimental condition was to study the effect of change in sand- MnO2(s) ratio on As(III) oxidation and adsorption. This experiment was designated as Exp C. ? Experiment with very low MnO2(s) content: For any MnO2(s) content above a certain threshold level, As(III) is completely oxidized to As(V). In that case, the model?s ability to accurately predict oxidation of As(III) cannot be confirmed. Therefore an experiment with sufficiently low MnO2(s) content was designed, leading to only partial oxidation of As(III). This experiment is designated as Exp D. Since sand and MnO2(s) have different porosity values and MnO2 content varies from experiment to experiment, the pore-volume needs to be calculated for each experiment individually. 69 0 0.2 0.4 0.6 0.8 1 1.2 0 5 10 15 20 25 30 35 V/V0 C/ C 0 prediction curve Exp B Exp A Figure 5 . 1 : For the column fully packed with MnO2(s), comparing the results of As(III) (Exp A) and As(V) (Exp B) transport with the model predictions; complete oxidation of As(III). 0 0.2 0.4 0.6 0.8 1 1.2 0 5 10 15 20 25 30 35 V/V0 C/ C 0 prediction curve Exp C Figure 5 . 2 : For a column with MnO2(s) 50% by wt., comparison of experimental and model results for As(III) transport (Exp C); complete oxidation of As(III). 70 0 0.2 0.4 0.6 0.8 1 1.2 0 5 10 15 20 25 30 35 V/Vo C/ Co prediction curve As total prediction curve As(V) Exp D As(V) Figure 5 . 3 : For a column with very low MnO2(s) content, comparison of experimental and model results for As(III) transport (Exp D); partial oxidation of As(III) 5.5 Results and discussion Based on the results of batch experiments, a conceptual model was proposed. This was translated into a reactive transport model to be able to predict the results of column experiments. Parameters obtained from the batch experiments were scaled to the column scale by taking into account the change in solid to solution ratio. The results of experiments A and B are compared with the model predictions in the Figure 5.1. As(V) breakthrough profiles of experiment A and B match very closely. Since As(III) was completely oxidized to As(V), only As(V) breakthrough profiles are presented. Even in the case of experiment C, MnO2(s) concentration was high enough to oxidize As(III) completely to As(V). In Figure 5.2, the results of experiment C are compared with the 71 model predictions. The model again provides a very good prediction of the experimental results. In case of experiment D, due to low MnO2 content As(III) was only partially oxidized to As(V). As discussed earlier, experiment D was designed in such a way that only partial oxidation of As(III) occurs. Astotal and As(V) profiles obtained from experiment D show good match with the model predictions in Figure 5.3. Results from the experiments A, B, C, and D demonstrate that the developed model properly accounts for both oxidation as well as adsorption of As by MnO2(s). 5.6 Summary In this study, a scalable reactive transport model was developed to describe the fate and transport of As species in manganese dioxide containing soil columns. The model involves two mobile reactive species that include As(III) and As(V). The governing transport equations along the sorption reaction were non-linearly coupled through the reaction terms. A fully implicit finite-difference approach was used to solve the coupled set of partial differential equations. A Picard-iterative scheme was used successfully to deal with the nonlinear sorption term. A synthetic form of commonly occurring manganese dioxide, pyrolusite, was used as the reactant in the soil columns. The oxidation rate constant and adsorption capacity were obtained from the batch experiments. Based on the information obtained from the batch experiments, a conceptual model and a matching numerical model were developed for predicting As transport in one-dimensional soil columns. The ability of the 72 numerical model to accurately predict As transport was verified by designing several column experiments. In the present study, the porous media was well characterized and the flow was uniform and one dimensional. Due to these idealistic conditions, a simple scaling approach to move from batch to column based on the solid-solution ratio appears to be effective. However in the real field scale scenarios, chemical heterogeneities cannot be easily characterized and in addition the flow field will be non-uniform influenced by physical heterogeneities. Therefore, it will be difficult to attain the necessary parameter distribution to directly employ the scaling approach developed in this work. Nevertheless, this effort is certainly a first major step to developing scalable model for accurately predicting the transport of As species in groundwater systems. 73 CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS 6.1 CONCLUSIONS Based on the outcome of this research, the following conclusions can be made: 1) To model reactive transport problems involving equilibrium geochemical reactions, an equilibrium geochemistry module should be coupled to a transport code. In this research effort, a suite of reactive transport models were built by coupling MICROQL, MINEQL and MINTEQA2 equilibrium geochemistry modules with a fully implicit, finite-difference, advection-dispersion solver. 2) The one-dimensional reactive transport models were tested by solving four benchmark problems which were identified as a part of this study. These benchmark problems simulate four distinct equilibrium-controlled processes that involve ion-exchange, surface-complexation, precipitation-dissolution, and redox reactions. 3) In terms of code length and complexity, the equilibrium geochemistry modules can be arranged as follows: MICROQL < MINEQL < MINTEQA2. The MICROQL coupled code is recommended for solving ion-exchange and surface- complexation reactions, MINEQL coupled code for precipitation-dissolution and redox reactions and MINTEQA2 coupled code for any combination of these reactions. This provides the user the necessary flexibility for selecting the model 74 depending on the type of equilibrium reactions. This flexibility reduces the code execution time and simplifies code modification tasks. 4) MINTEQA2 was implemented within RT3D to develop a prototype version of a coupled model. 5) The prototype model was validated by solving several one-dimensional problems, including a complex benchmark problem that includes precipitation-dissolution and redox reactions. This prototype model was also integrated into the groundwater user interface GMS (Groundwater Modeling System). 6) A scalable numerical model was developed to describe As speciation coupled with oxidation in MnO2(s) soils. 7) Finally, the model performance was validated by reproducing a set of column- scale experimental data using scaled parameters obtained from the batch experiments. 6.2 RECOMMENDATIONS FOR FUTURE WORK 1) The current version of MINTEQA2-coupled RT3D model is a prototype and it needs to be further tested by solving multi-dimensional reactive transport problems. All three codes including MICROQL, MINEQL and MINTEQA2 should be integrated within a single framework to allow user to be able to select and run any of these models. 2) A user interface is required to input the data necessary for modeling equilibrium reaction processes. 75 3) After fully integrating MINETQA2 within RT3D the code will have the capability to individually model both equilibrium and kinetics types of chemical reactions. However, it still cannot solve reactive transport problems involving partial equilibrium reactions where some of the reactions are kinetic and some are at equilibrium. Future efforts should focus on developing this capability. 4) The surface complexation module in MINTEQA2 requires further testing and debugging before it can be fully integrate with RT3D. 5) Two-dimensional soil box experiments should be completed to further study Arsenic reactive transport in manganese dioxide containing systems. This data will help test the validity of the scaling approach in multiple dimensions. 76 REFERENCES 1. Abrams, R. H., et al. (1998), Development and testing of a compartmentalized reaction network model for redox zones in contaminated aquifers, Water Resources Research, 34, 1531-1542. 2. Allison, J. D., et al. (1990), MINTEQA2/PRODEFA2, A geochemical assessment model for environmental systems: Version 3.0. 3. Appelo, C. A. J., et al. (1997), A hydrogeochemical transport model for an oxidation experiment with pyrite/calcite/exchangers/organic matter containing sand, Applied Geochemistry, 13, 257-268. 4. Ball, J. W., et al. (1981), WATEQ3 A geochemical model with uranium added. 5. Barry, D. A., et al. (2002), Modelling the fate of oxidisable organic contaminants in groundwater, Advances in Water Resources, 25, 945-983. 6. Bethke, C. (1994), The Geochemist's Workbench. A User's Guide to Rxn, Act2, Tact, React and Gtplot. University of Illinois, Urbana-Champaign. 7. Bryant, S. L., et al. (1986), Interactions of precipitation/dissolution waves and ion exchange in flow through porous media, AIChE J, 32, 751-764. 8. Carnahan, C. L. (1988), Some Effects of Data-Base Variations on Numerical Simulations of Uranium Migration, Radiochimica Acta, 44-5, 349-354. 9. Cederberg, G. A., et al. (1985), A Groundwater Mass-Transport and Equilibrium Chemistry Model for Multicomponent Systems, Water Resources Research, 21, 1095- 1104. 10. Chilakapati, A. (1995), RAFT: A simulator for reactive flow and transport of groundwater contaminants, PNL Rep. 10636, Pac, Northwest Lab., Richland, Wash. 11. Chilakapati, A., et al. (2000), Groundwater flow, muticomponent transport and biogeochemistry: Development and application of a coupled process model, Journal of Contaminant Hydrology, 43, 303-325. 77 12. Clement, T. P. (1997), RT3D - A modular computer code for simulating reactive multi-species transport in 3-Dimensional groundwater aquifers, Battelle Pacific Northwest National Laboratory Research Report, PNNL-SA-28967. Available at: http://bioprocess.pnl.gov/rt3d.htm. 13. Clement, T. P., et al. (1998), Modeling Multi-species Reactive Transport in Groundwater Aquifers, Groundwater Monitoring & Remediation Journal, 18, 79-92. 14. Engesgaard, P., and K. L. Kipp (1992), A Geochemical Transport Model for Redox- Controlled Movement of Mineral Fronts in Groundwater-Flow Systems - a Case of Nitrate Removal by Oxidation of Pyrite, Water Resources Research, 28, 2829-2843. 15. Felmy, A. R., et al. (1983), MINTEQ: A computer program for calculating aqueous geochemical equilibria, report, 62 pp, USEPA, Washington, D.C. 16. Frazier, M., and G. Johnson (2005), U.S. Departmrnt of Energy, Office of Science http://doegenomestolife.org/cleanup.pdf. 17. Freeze, R. A., and J. A. Cherry (1979), GROUNDWATER, Prentice Hall, Englewood Cliffs, NJ 07362. 18. Herzer, J., and W. Kinzelbach (1989), Coupling of Transport and Chemical Processes in Numerical Transport Models, Geoderma, 44, 115-127. 19. Hundsdorfer, W., and J. G. Verwer (1995), A note on splitting errors for advection- reaction equations, Applied Numer. Math., 18, 191-199. 20. Jennings, A. A., et al. (1982), Multicomponent equilibrium chemistry in groundwater quality models, Water Resources Research, 18, 1089-1096. 21. Kirkner, D. J., and H. Reeves (1988), Multicomponent Mass-Transport with Homogeneous and Heterogeneous Chemical-Reactions - Effect of the Chemistry on the Choice of Numerical Algorithm .1. Theory, Water Resources Research, 24, 1719-1729. 22. Kirkner, D. J., et al. (1985), Multicomponent mass transport with chemical interaction kinetics, Journal of Hydrology, 76, 107-117. 23. Lewis, F. M., et al. (1987), Solute transport with equilibrium aqueous complexation and either sorption or ion exchange: Simulation methodology and applications, Journal of Hydrology, 90, 81-115. 24. Li, L., et al. (2006), Modeling porosity reductions caused by mineral fouling in continuous-wall permeable reactive barriers, Journal of Contaminant Hydrology, 83, 89- 121. 78 25. Lichtner, P. C. (1985), Continuum Model for Simultaneous Chemical-Reactions and Mass-Transport in Hydrothermal Systems, Geochimica Et Cosmochimica Acta, 49, 779- 800. 26. Liu, C. W., and T. N. Narasimhan (1989), Redox-Controlled Multiple-Species Reactive Chemical-Transport .2. Verification and Application, Water Resources Research, 25, 883-910. 27. Mayer, K. U. (1999), A numerical model for multicomponent reactive transport in variably saturated porous media. PhD thesis, University of Waterloo, Waterloo, ON, Canada. 28. Mayer, K. U., et al. (2002), Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions, Water Resources Research, 38. 29. McDuff, R. E., and F. M. M. Morel (1972), REDEQL, a general program for the computation of chemical equilibrium in aqueous systems. Keck Laboratory of Environmental Engineering Science, California Institute of Technology, Pasadena, CA. 30. McNab, W. W., and T. N. Narasimhan (1994), Modeling Reactive Transport of Organic-Compounds in Groundwater Using a Partial Redox Disequilibrium Approach, Water Resources Research, 30, 2619-2635. 31. Miller, C. W., and L. V. Benson (1983), Simulation of solute transport in a chemically reactive heterogeneous systems: Model development and application, Water Resources Research, 19, 381-391. 32. Narasimhan, T. N., et al. (1986), Groundwater Contamination from an Inactive Uranium Mill Tailings Pile .2. Application of a Dynamic Mixing Model, Water Resources Research, 22, 1820-1834. 33. Noorishad, J., et al. (1987), Development of the nonequilibrium reactive chemical transport code CHMTRNS, Rep. LBL-22361, Lawrence Berkeley Lab., Univ. of Calif., Berkeley. 34. Parkhurst, D. L. (1995), User's guide to PHREEQC, A computer model for speciation, reaction-path, advective-dispersive transport and inverse geochemical calculations U.S. Geol. Surv. Water Resour. Invest. Rep. 95-4227. 35. Parkhurst, D. L., and C. A. J. Appelo (1999), User's guide to PHREEQC (version 2) - A computer program for speciation, batch-reaction, one-dimensional transport, reaction path, and inverse geochemical calculations, U.S. Geol. Surv. Water Resour. Invest. Rep., 99, 4259, 312 pp. 36. Parkhurst, D. L., et al. (1980), PHREEQE - A computer program for geochemical calculations, U.S. Geol. Surv. Water Resour. Invest. Rep., 80-96, 195 pp. 79 37. Phanikumar, M. S., and J. T. McGuire (2004), A 3D partial-equilibrium model to simulate coupled hydrogeological, microbiological, and geochemical processes in subsurface systems, Geophysical Research Letters, 31. 38. Prommer, H., et al. (2003), MODFLOW/MT3DMS-Based Reactive Multicomponent Transport Modeling, Ground Water, 41, 247-257. 39. Radu, T. (2006), Transport of arsenic in experimental subsurface systems, PhD Thesis Draft, Auburn University. 40. Rubin, J. (1983), Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions, Water Resources Research, 19, 1231-1252. 41. Rubin, J., and R. V. James (1973), Dispersion-affected transport of reacting solutes in saturated porous media: Galerkin method applied to equilibrium-controlled exchange in unidirectional steady water flow, Water Resources Research, 9, 1332-1356. 42. Runkel, R. L., et al. (1996), Reactive solute transport in streams .1. Development of an equilibrium-based model, Water Resources Research, 32, 409-418. 43. Runkel, R. L., et al. (1999), Reactive solute transport in streams: A surface complexation approach for trace metal sorption, Water Resources Research, 35, 3829- 3840. 44. Salvage, K. M., and G. T. Yeh (1998), Development and application of a numerical model of kinetic and equilibrium microbiological and geochemical reactions (BIOKEMOD), Journal of Hydrology, 209, 27-52. 45. Salvage, K. M., et al. (1996), Development of a model of subsurface hydrologic transport and biogeochemical reactions (HYDROBIOGEOCHEM), in Computational Methods in Water Resources, XI, vol. 2, Computation Methods in Surface Flow and Transport Problems, pp. 517-524, Comput. Mech., Boston, Mass. 46. Smith, S. L., and P. R. Jaffe (1998), Modeling the transport and reaction of trace metals in water-saturated soils and sediments, Water Resources Research, 34, 3135-3147. 47. Steefel, C. I., and A. C. Lasaga (1990), Evolution of dissolution patterns, permeability changes due to coupled flow and reaction, in Chemical Modeling of Aqueous Syetems II, edited by D. C. Melchior and R. L. Bassett, pp 213-225, American Chemical Society, Washington D. C. 48. Steefel, C. I., and A. C. Lasaga (1994), A Coupled Model for Transport of Multiple Chemical-Species and Kinetic Precipitation Dissolution Reactions with Application to Reactive Flow in Single-Phase Hydrothermal Systems, American Journal of Science, 294, 529-592. 80 49. Steefel, C. I., and S. B. Yabusaki (1996), OS3D/GIMRT, Software for modeling multi-component-multidimensional reactive transport, user's manual and programmer's guide, PNL-11166, Pac. Northwest Lab., Richland, Wash. 50. Tebes-Stevens, C., et al. (1998), Multicomponent transport with coupled geochemical and microbiological reactions: model description and example simulations, Journal of Hydrology, 209, 8-26. 51. Theis, T. L., et al. (1982), Multi-solute subsurface transport modeling for energy solid wastes, Tech. Progr. Rep., C00-10253-3, Dep. of Civ. Eng., Univ. of Notre Dame, Notre Dame, Ind. 52. Truesdell, A. H., and B. F. Jones (1974), WATEQ, A computer program for calculating chemical equilibria of natural water, U.S. Geol. Surv. J. Res., 2, 233-248. 53. USGS (2000), Web Page http://ga.water.usgs.gov/edu/wugw.html. 54. valocchi, A. J., et al. (1981), Transport of ion-exchanging solutes in groundwater: Chromatographic theory and fiels simulations, Water Resources Research, 17, 1517- 1527. 55. van der Lee, J. (1998), Thermodynamic and mathematical concepts of CHESS, Technical Report LHM/RD/98/39, CIG, Ecole de Mines de Paris, Fontainbleau, France. September. 56. Walsh, M. P., S. L. Bryant, and L. W. Lake (1984), Precipitation and dissolution of solids attending flow through porous media, AIChE J, 30, 317-328. 57. Walter, A. L., et al. (1994), Modeling of Multicomponent Reactive Transport in Groundwater .1. Model Development and Evaluation, Water Resources Research, 30, 3137-3148. 58. Westall J. C., et al. (1976), MINEQL: A computer program for the calculation of the chemical equilibrium composition of aqueous system, Tech Note 18, 91 pp, Dept. of Civ. Eng., Mass. Inst. of Technol., Cambridge. 59. Westall, J. C. (1979a), "MICROQL. I. A Chemical Equilibrium Program in BASIC," Technical Report, Swiss Federal Institute of Technology, EAWAG, Dubendorf, Switzerland. 60. Westall, J. C. (1979b), "MICROQL. II. Computation of Adsorption Equilibria in BASIC," Technical Report, Swiss Federal Institute of Technology, EAWAG, Dubendorf, Switzerland. 61. Wolery, T. J. (1992), EQ3/6, A software package for geochemical modeling of aqueous systems: Package overview and installation guide, UCRL-MA-110662 PT 1, Lawrence Livermore National Laboratory, California. 81 62. Xu, T. (1996), MODELING NONISOTHERMAL MULTICOMPONENT REACTIVE SOLUTE TRANSPORT THROUGH VARIABLY SATURATED POROUS MEDIA PhD Thesis, University of La Coruna. 63. Yeh, G. T., M. D. Siegel, and M. H. Li (2001), Numerical modeling of coupled fluid flows and reactive transport including fast and slow chemical reactions, Journal of Contaminant Hydrology, 47, 379-390. 64. Yeh, G. T., and V. S. Tripathi (1988), Hydrogeochem: A coupled model of HYDROgeological transport and GEOCHEMical Equilibrium of Multi-component Systems, Rep. ORNL-6371, Oak Ridge Nat. Lab., Oak Ridge, Tenn. 65. Yeh, G. T., and V. S. Tripathi (1989), A Critical-Evaluation of Recent Developments in Hydrogeochemical Transport Models of Reactive Multichemical Components, Water Resources Research, 25, 93-108. 66. Yeh, G. T., and V. S. Tripathi (1991), A Model for Simulating Transport of Reactive Multispecies Components - Model Development and Demonstration, Water Resources Research, 27, 3075-3094. 67. Zheng, C., and P. P. Wang (1999), MT3DMS, A modular three-dimensional multi- species transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and user?s guide, U.S. Army Engineer Research and Development Center Contract Report SERDP-99-1, Vicksburg, MS, 202 pp.