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.