Published online 25 February 2008
Published in Vadose Zone J 7:305-315 (2008)
DOI: 10.2136/vzj2006.0130
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: TOUGH2
Incorporating Aqueous Reaction Kinetics and Biodegradation into TOUGHREACT: Applying a Multiregion Model to Hydrobiogeochemical Transport of Denitrification and Sulfate Reduction
Tianfu Xu*
Earth Sciences Division, Lawrence Berkeley National Lab., One Cyclotron Rd., Berkeley, CA 94720
* Corresponding author (tianfu_xu{at}lbl.gov).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 5 September 2006.
 |
ABSTRACT
|
|---|
The need to consider aqueous and sorption reaction kinetics and microbiological processes arises in many subsurface problems. By adding these processes to the TOUGH family code, more complex problems can be addressed that involve multiphase fluid and heat flow, and geochemical interaction. This paper presents a formulation for incorporating kinetic rates among primary species into mass-balance equations. The space discretization used is based on a flexible integral finite difference approach that uses irregular gridding to model biogeologic structures. A general multiregion model for hydrological transport interacted with microbiological and geochemical processes is proposed. A one-dimensional reactive transport problem with kinetic biodegradation and sorption was used to test the enhanced simulator, which involves the processes that occur when a pulse of water containing nitrylotriacetate (NTA) and cobalt is injected into a column. The current simulation results agree very well with those obtained with other simulators. The applicability of this general multiregion model was validated by results from a published column experiment of denitrification and sulfate reduction. The matches with measured nitrate and sulfate concentrations were adjusted with the interfacial area between mobile and immobile regions. Results suggest that TOUGHREACT not only is useful interpretative tool for biogeochemical experiments but also can produce insight into processes and parameters of microscopic diffusion and their interplay with biogeochemical reactions.
Abbreviations: DOC, dissolved organic carbon • DOM, dissolved organic matter • EPS, extracellular polymeric substances • IFD, integral finite differences • MINC, multiple interacting continua • NTA, nitrylotriacetate
|
INTRODUCTION
|
|---|
In natural systems, many redox reactions (such as SO42–/HS– and CO2/CH4) that take place homogeneously within the aqueous phase are slow to achieve equilibrium and therefore violate the local equilibrium assumption. Studies of many subsurface problems, including acidic mine drainage remediation, organic matter decomposition, nitrate cycle, oil and gas maturation, sulfate reduction in oil fields, and effective environmental remediation of groundwater contamination, require biodegradation simulation capabilities. Here, the formulation for intra-aqueous kinetic reactions and biodegradation has been incorporated into the general numerical simulation program for chemically reactive nonisothermal flows of multiphase fluids in porous and fractured media, TOUGHREACT (Xu and Pruess, 2001; Xu et al., 2006). Most previous models with reaction kinetics have simplified descriptions of groundwater flow. By adding intra-aqueous reaction kinetics and biodegradation to the TOUGH family code, more complex problems can be addressed that involve in particular multiphase fluid and heat flow, and geochemical interaction.
The TOUGHREACT program is written in Fortran 77 and was developed by introducing reactive geochemical transport processes into the multiphase fluid and heat flow simulator TOUGH2 (Pruess et al., 1999). A variety of subsurface thermo-physical-chemical processes have been already considered under a wide range of conditions of pressure, temperature, water saturation, ionic strength, pH, and redox potential (Eh). Interactions between mineral assemblages and fluids can occur under local equilibrium or kinetic rates. The gas phase can be chemically active. Precipitation and dissolution reactions can change formation porosity and permeability. The program can be applied to many geologic systems and environmental problems, including geothermal systems, diagenetic and weathering processes, subsurface waste disposal, acid mine drainage remediation, contaminant transport, and groundwater quality. The capabilities for reactive fluid flow and geochemical transport have been demonstrated in Spycher et al. (2003), Sonnenthal et al. (2005), and Xu et al. (2006).
The added process capabilities for intra-aqueous kinetic reactions and biodegradation were first tested with other simulators for a one-dimensional reactive transport problem. The flexible space discretization, based on an integral finite difference approach, can use irregular gridding to model biogeologic structures. TOUGHREACT can deal directly with general multiregion and multiple interacting continua models for saturated and unsaturated porous and fractured media. To test the applicability of the multiregion model embodied in TOUGHREACT to reactive transport of denitrification and sulfate reduction, the column experiments of von Gunten and Zobrist (1993) were modeled. The experimental data were first used by Zysset et al. (1994) for their macroscopic model involving the transport of dissolved substances in groundwater–biofilm systems. Zysset's model conceptualizes the diffusion-dominated microscopic transport processes within the biofilm by using a logistic approach based on a diffusional limitation. Biological processes were explicitly taken into account. Later, Doussan et al. (1997) used the same experimental data set for their model testing. They treated the diffusion dominated microscopic transport processes within the biofilm by using a mass-transfer coefficient in the macroscopic transport equations.
The objectives of this paper are (i) to present a general biogeochemical transport program that fully uses the comprehensive capabilities of TOUGHREACT on nonisothermal multiphase fluid flow and complex geochemical reactions, and (ii) to present a general multiregion model for producing insight into microscopic processes and parameters for denitrification and sulfate reduction.
|
Formulation
|
|---|
Flow and transport in geologic media are modeled based on space discretization by means of integral finite differences (IFD) (Narasimhan and Witherspoon, 1976). The IFD method provides a flexible discretization for geologic media that allows the use of irregular grids, which is well suited for simulation of flow, transport, and fluid-microbial-rock interactions in multi-region heterogeneous and fractured rock systems (see section "Multiregion Model" for details). For regular grids, IFD is equivalent to conventional finite differences. An implicit time-weighting scheme is used for the individual components of flow, transport, and kinetic biogeochemical reactions. TOUGHREACT uses a sequential iteration approach similar to Yeh and Tripathi (1991) and 
im
nek and Suarez (1994). The transport of chemical and microbial species is solved on a component-by-component basis. The system of mixed equilibrium–kinetic biogeochemical reaction equations is solved on a grid-block–by–grid-block basis by Newton–Raphson iteration. The governing equations for fluid and heat flow and geochemical transport are given in Xu and Pruess (2001). Here, (i) mass-balance equations for batch biogeochemical system and (ii) the rate law used for aqueous reaction kinetics and biodegradation are presented.
Mass Balance for Batch Biogeochemical System
The TOUGHREACT formulation for an aqueous equilibrium system was presented in Xu and Pruess (2001), which is based on mass balance in terms of primary (basis) species. In contrast to aqueous equilibrium, species involved in kinetic reactions, such as redox couples, are independent and must be considered as primary species (Steefel and MacQuarrie, 1996). For example, for the reaction
 | [1] |
under kinetic conditions, both HS– and SO42– must be placed in the primary species list (Steefel, 1997). Thus, all redox reactions making use of these species must be decoupled in the input of the thermodynamic database. Based on the previous formulation (Xu and Pruess, 2001), a general rate expression for intra-aqueous kinetic reaction and biodegradation has been incorporated into TOUGHREACT (see the next section).
At time zero (initial), the total concentrations of primary species j (Tj0) are assumed to be known and are given by
 | [2] |
where superscript 0 represents time zero; c is concentration; subscripts j, k, m, and n are the indices of primary species, aqueous complexes, minerals at equilibrium, and minerals under kinetic constraints, respectively; Nc, Nx, Np, and Nq are the number of corresponding species and minerals; and 
kj, mj, and nj are stoichiometric coefficients of the primary species in the aqueous complexes, equilibrium, and kinetic minerals, respectively.
After a time step 
t, the total concentration of primary species j (Tjt) is given by
 | [3] |
where rn is the kinetic rate of mineral dissolution and precipitation (positive for dissolution and negative for precipitation; units used here are moles of mineral per kilogram of water per time), for which a general multimechanism rate law was used (Xu et al., 2006). Tj0 and Tjt are related through generation of j among primary species as follows:
 | [4] |
where l is the aqueous kinetic reaction (including biodegradation) index, Na is total number of kinetic reactions among primary species, and rl is the kinetic rate which is in terms of one mole of product species, such as SO42– per unit time. Therefore, for product species, the stoichometric coefficients lj are positive, and for reactant species they are negative. For Eq. [1], lj for SO42– and H+ are 1; for HS– it is –1, and for O2(aq) it is –2.
By substituting Eq. [2] and Eq. [3] into Eq. [4], and denoting residuals as Fjc (which are zero in the limit of convergence), we have
 | [5] |
According to mass-action equations, concentrations of aqueous complexes ck can be expressed as functions of concentrations of the primary species cj. Kinetic rates rn and rl are functions of cj. The expression for rn is given in Xu et al. (2006); rl is presented in the next section. No explicit expressions relate equilibrium mineral concentrations cm to cj. Therefore, Np additional mass-action equations (one per mineral) are needed. Details on solution of the nonlinear system of equations by the Newton–Raphson iterative method are given in Xu et al. (1999). Notice that gas dissolution–exsolution, cation exchange, and surface complexation are included in TOUGHREACT, but for simplicity, the formulation is not shown here.
Rate Expressions
Following the expression of Curtis (2003) and adding multiple mechanisms (or pathways), a general-rate law for reactions among primary species (including intra-aqueous and sorption reaction kinetics and biodegradation) is used:
 | [6] |
where ri is the reaction rate of the ith kinetic reaction, Mech is the number of mechanisms or pathways and s is the mechanism counter, ki is a rate constant (often denoted vmax, maximum specific growth constant for biodegradation), 
j is the activity coefficient of species j, Cj is the concentration of species j, ni,j are power terms, Nl is the number of reacting species used in the forward rate term (called product terms), Nm is the number of Monod factors (Monod terms), Ci,k is the concentration of the kth Monod species, Ci,p is the concentration of the pth inhibiting species, KMi,k is the kth Monod half-saturation constant of the ith species, Np is the number of inhibition factors (inhibition terms), and Ii,p is the pth inhibition constant. Equation [6] accounts for multiple mechanisms and multiple products, Monod, and inhibition terms, which can cover many rate expressions (examples of rate expression will be given in sections below, "Kinetic Biodegradation and Sorption" and "Denitrification and Sulfate Reduction").
|
Multiregion Model
|
|---|
In vadose zones and groundwater systems, most redox reactions are catalyzed by bacteria. Most of these bacteria are fixed on the solid phase within aquifers (Harvey et al., 1989; Zysset et al., 1994) or in river sediments (van Beelen and Fleuren Kemila, 1989). Formation around bacteria of microzones of polysaccharidic gels and other bacteria limits the dispersion of metabolites and nutrients (Doussan et al., 1997). It may also protect microorganisms against desiccation or hagocytosis, toxins, and so on. This conglomerate, consisting of the fully hydrated polymeric gel and bacteria, is called biofilm. Most bacteria grow within an immobile biofilm region residing on the solid surface.
A multiregion concept was adopted here for modeling the present biogeochemical transport problem. Gwo et al. (1995) presented a multiple-pore-region conceptual model (in short or multiregion model) to account for pore structures as well as the resultant widely distributed pore water velocities in porous media. Pore regions can either be physically identified as discrete features, such as fractures and rock matrices, or be experimentally determined by separation of water retention curves according to pore classification schemes. A multiregion mechanism is proposed to account for the effect of local-scale and field-scale heterogeneities on mass transport under variably saturated conditions. The applicability of the concept was demonstrated using laboratory soil column tracer injection data. Later, the multiregion concept was used by other investigators such as imnek et al. (2003). Notice that the concept of region here is not a "regional model" that generally imply a scale measured in square kilometers. The concept being used is a sub-REV concept.
In this paper, we use a general three-region model (Fig. 1
) consisting of (i) a mobile region (a fraction of the porosity), (ii) an immobile region (another fraction of the porosity), and (iii) a solid particle region where mineral dissolution–precipitation and surface reactions may occur. Here, the immobile region contains stagnant water and biomass. According to Thullner et al. (2004), only minor amount of the total organic carbon was present as bacterial biomass, whereas the most remaining were attributed to extracellular polymeric substances (EPS). Due to EPS production, a relatively small amount of bacteria may change the volume of the immobile region and the hydraulic properties of the porous media. Currently, TOUGHREACT does not consider the dynamic changes in the volume of immobile region. The dynamic changes in porosity and volume will be implemented in the future. Several models (Clement et al., 1996; Thullner et al., 2004) have been published to simulate the experimentally observed interaction between biomass growth and changes in hydraulic properties of porous media. Note that TOUGHREACT can handle the transfer of bacteria between mobile and immobile regions if one wants to simulate this process.
In a principle similar to the multiregion model, Pruess and Narasimhan (1985) presented a method of "multiple interacting continua" (MINC) for resolving "global" fluid and heat flow and diffusion of chemicals in fractured rock and their interaction with the "local" exchange between fracture and matrix rock. This method was first developed for fluid and heat flow in fractured porous media. The extension of the MINC method to reactive geochemical transport is described by Xu and Pruess (2001). Similarly, a "shrinking core" model was proposed for sulfide mineral oxidation by Wunderly et al. (1996). As the reaction between oxygen and sulfide minerals within the particles progresses, the radius of the unreacted core will gradually decrease, while the thickness of the oxidized shell increases. Thus, the oxygen flux from the outside of the particle surface to the unreacted core decreases with time.
In general, it is not necessary to consider explicitly every pore and particle in a macrogrid block. All mobile portions of pores are lumped into Continuum no. 1. All immobile portions of pores are lumped into Continuum no. 2. All solid particle volumes are lumped into Continuum no. 3. If necessary, each region can be further divided. For example, geochemical reactions often occur close to the particle surface, such that the portion close to the surface can be considered as one solid region and the inner portion of the particle can be another solid region. The regions are specified by means of a set of volume fractions VOL(j), j = 1,..., J (where j is the region counter) into which the "primary" porous medium (macro-)grid blocks are partitioned. For microscopic diffusion processes, distances (d) and interfacial areas (A) between mobile and immobile regions are required to be specified or can be calibrated by matching observation data. The shapes, sizes, and distributions of pores and particles are not required for the model.
|
Kinetic Biodegradation and Sorption
|
|---|
Problem Setup
A biodegradation and sorption problem was used for verifying TOUGHREACT against other simulators. The main focus of this verification problem is on model results rather than experimental conditions. A column model of reactive biogeochemical transport originally developed by Tebes-Stevens et al. (1998) was employed. This problem has also been used by others for verification of PHREEQC (Parkhurst and Appelo, 1999) and Bio-CORE2D (Zhang, 2001). This problem has several interacting chemical processes common to many environmental problems: biologically mediated degradation of an organic substrate, bacterial cell growth and decay, metal sorption, and aqueous speciation including metal-ligand complexation. The problem involves the transport processes when a pulse of water containing NTA and cobalt is injected into a column. The problem includes advection and dispersion in the column, aqueous equilibrium reactions, and kinetic reactions for NTA degradation, growth of biomass, and cobalt sorption. The dimension and hydrological properties of the column given in Table 1
were from Parkhurst and Appelo (1999). Aqueous complexes considered and their stoichiometric coefficients and equilibrium constants are summarized in Table 2
.
Biodegradation
The initial concentrations in the column are listed in Table 3
. The column contains no NTA and cobalt initially but has a biomass of 1.36 x 10–4 g/L. Injection water with NTA and cobalt is applied at the inlet of the column for 20 h. Injection concentrations are also given in Table 3. After 20 h, the background (initial) water is introduced at the inlet until the experiment ends at 75 h. Nitrylotriacetate is assumed to degrade in the presence of biomass and oxygen according to the following reaction (Parkhurst and Appelo, 1999):
 | [7] |
TOUGHREACT requires kinetic reactants and products to be defined solely by the number of mole or milligrams or grams of each that enter or leave the solution. Equation [7] converts 1 mol HNTA2– and 1.62 mol O2(aq) to 0.576 mol C5H7O2N and 65.14 g biomass. For biomass, the generation amount is in grams, so it is not explicitly written in the right-hand side of Eq. [7]. The yield coefficient is not directly used and must be converted into moles or grams or milligrams or other units. The following double Monod rate expression is used to describe the rate of NTA degradation:
 | [8] |
where rHNTA2-is the rate of HNTA2– degradation (mole per liter per hour; it has a negative value indicating consumption), kbiom is the maximum specific rate of substrate utilization (mol/g cells/h), Xbiom is the biomass (g cells/L), Ks is the half-saturation constant for the substrate NTA (mol/L), Ka is the half-saturation constant for the electron acceptor O2 (mol/L), and C denotes concentration of species (mol/L). Equation [8] compared to the general form, Eq. [6], has only one mechanism, one product term Xbiom, and two Monod terms.
View this table:
[in this window]
[in a new window]
|
TABLE 3. List of concentrations for injection and background waters (Parkhurst and Appelo, 1999). A dash (–) is placed in the table for injection concentrations of Biomass, CoNTA(ads), and Co(ads), indicating not subject to transport. (ads = adsorption of surface species.)
|
|
The rate of biomass production is dependent on the rate of substrate utilization and a first-order decay rate for the biomass:
 | [9] |
where rbiom is the rate of cell growth (gram cells per liter per hour), y is the microbial yield coefficient (gram cells per mole NTA), and b is the first-order biomass decay coefficient (per hour). The parameter values for Eq. [8] and [9] are listed in Table 4
. Equation [9] consists of two mechanisms using Eq. [6].
Kinetic Sorption
The rate expressions for kinetic sorption reactions involving Co2+ and CoNTA– are given by
 | [10] |
where i is either Co2+ or CoNTA– (mole per liter), Si is the sorbed concentration (mole per gram sediment), km is the mass-transfer coefficient (per hour), and Kd is the distribution coefficient (liter per gram). Values for the coefficients are given in Table 5
. The values of Kd were defined to give retardation coefficients of 20 and 3 for Co2+ and CoNTA–, respectively. Two mechanisms involved in the sorption kinetics are each represented by a product term, using Eq. [6].
Comparison of Results
The evolution of aqueous and immobile constituents at the outlet of the column is shown in Fig. 2
. In the experiment, two pore volumes of water with NTA and cobalt were introduced into the column over the first 20 h, followed by 5.5 pore volumes of background water over the next 55 h. At 10 h, HNTA2– begins to appear at the column outlet, along with a rise in the pH (Fig. 2a). If NTA and cobalt were conservative and dispersion negligible, the graph would show square pulses that increase at 10 h and decrease at 30 h. However, the movement of the NTA and cobalt is retarded (relative to conservative movement) by the sorption reactions. The peak in NTA and cobalt concentrations occurs in the CoNTA complex between 30 and 40 h (Fig. 3c and 3d
), while the peak in Co2+ concentration is even more retarded by its sorption reaction and does not show up until near the end of the experiment. The sorbed CoNTA– concentration peaks between 30 and 40 h and slightly lags behind the peak in the dissolved concentration of the CoNTA– complex. Initially, no NTA is present in the column, and the biomass decreases slightly over the first 10 h because of the first-order decay of the biomass (Fig. 2b). As the NTA moves through the column, the biomass increases as the NTA substrate becomes available. After the peak of NTA has moved through the column, biomass concentrations level off and then begin to decrease because of decay. The Kd for cobalt sorption relates to a greater retardation coefficient than the Kd for CoNTA– sorption, and the sorbed concentration of Co2+ appears to be still increasing at the end of the experiment. The TOUGHREACT simulation results generally agree well with those of Bio-CORE2D. Some minor differences can be observed at the peak concentrations of dissolved and sorbed CoNTA– due to different numerical dispersion behaviors. TOUGHREACT uses an integral finite difference for spatial discretization, while Bio-CORE2D uses a finite element method. TOUGHREACT may have slightly larger numerical dispersion. From this test problem, the added processes of biodegradation and kinetic sorption to TOUGHREACT have been verified.
|
Denitrification and Sulfate Reduction
|
|---|
Problem Statement
To test the applicability of TOUGHREACT to reactive transport of denitrification and sulfate reduction, the column experiments of von Gunten and Zobrist (1993) were modeled. The experimental data were first used by Zysset et al. (1994) for their macroscopic model involving the transport of dissolved substances in groundwater-biofilm systems. Later, Doussan et al. (1997) used the same experimental data set for their model testing. They treated the diffusion-dominated microscopic transport processes within the biofilm by using a mass-transfer coefficient in the macroscopic transport equations. These experiments were designed to simulate infiltration of an organically polluted river into an aquifer. Thus, synthetic river water, including an organic substrate (lactate) and electron acceptors of oxygen, nitrate, and sulfate, were injected into columns filled with river sediments. Results of two column experiments were reported (von Gunten and Zobrist, 1993). The two columns have the same size of diameter, 5 cm, and length, 29 cm. The model aquifer was a sand fraction (0.125–0.25 mm) from a fluvio glacial deposit. The sand contained only a small amount of organic material (<0.1%). The columns were kept at ambient temperature. Column 1 was inoculated at the inflow with an aqueous extract from a humus soil. Column 2 was inoculated with 1 g of material from the first 1 to 2 cm of Column 1. In the first column, only nitrate is added to the solution as an electron acceptor. In the second column, oxygen, nitrate, and sulfate are used as electron acceptors.
Biodegradation Kinetics and Parameters
As described by the previous investigators, three major microbially mediated reactions are involved in the experiments. Three electron acceptors are reduced, while dissolved organic matter (DOM), using lactate (C3H5O3–) or propionate (C3H5O2–) in the experiment, is oxidized as follows:
 | [11] |
 | [12] |
 | [13] |
The bacterial growth rates due to three different electron acceptors are given in Eq. [14], [15], and [16], respectively. Denitrification is inhibited by oxygen, and sulfate reduction is inhibited by both oxygen and nitrate. The rate parameters for Eq. [14–16] are given in Table 6
.
 | [14] |
 | [15] |
 | [16] |
The overall biomass growth rate is expressed as
 | [17] |
where Xb is biomass concentration (milligram per liter), b is decay constant, and a value of 5.787 x 10–7 1/s was used according to Doussan et al. (1997). In this example, biomass is assumed not subject to transport. As mentioned in the section "Multiregion Model" above, most of the bacteria are fixed on the solid phase within geologic media.
Denitrification
In the first column, only nitrate is added to the solution as an electron acceptor. As already mentioned, Doussan et al. (1997) modeled the reactive transport of denitrification in Column 1. They treated the diffusion-dominated microscopic transport processes within the biofilm by a mass-transfer coefficient, which was calibrated by the experimental data. Maximum specific growth rates and half-saturation constants were also calibrated. These calibrated parameters were also used for their Column 2 simulation. These parameter estimates were obtained with biomass growth. A yield coefficient of 0.4 was used in their simulations. In the present study, a porous medium model for the denitrification is first used. Then, the general multiregion model is applied.
Porous Medium Model
The porous one-dimensional flow system has a cross-sectional area of 1 m2 and 0.1 m length, divided into 50 grid blocks of 0.002 spacing. The experiments were started by inoculating the water in the columns with a small amount of bacteria, but no biomass concentration was reported in von Gunten and Zobrist (1993). Here, a starting biomass concentration of 0.15 mg/L calibrated by Zysset et al. (1994) was used. Initial concentrations for all other species are set equal to a very small value of 10–10 mg/L (practically zero, because TOUGHREACT uses log10 calculations for a better convergence). Constant concentrations at the boundary were used in the simulation, which are given in Table 7
. A diffusion coefficient of 1 x 10–9 m2/s was used. TOUGHREACT uses an upstream weighting scheme, resulting in a numerical dispersivity of 
n = x/2 = 0.001 m, which equals the actual hydrodynamic dispersivity used by Zysset et al. (1994). Hydrodynamic dispersion is an important solute transport mechanism that arises from an interplay between nonuniform advection and molecular diffusion. In geologic media, velocities of fluid parcels are spatially variable due to heterogeneities on multiple scales, all the way from the pore scale to basin scale. The process is often represented by a Fickian diffusion analog (convection–dispersion equation), which has fundamental flaws and limitations, as has been demonstrated in numerous studies in the hydrogeology literature of the last 20 yr (Xu and Pruess, 2001). Therefore, the dispersion process is not considered in TOUGHREACT. Alternatively one has to use a finer grid and multiregion model to resolve this process.
The results of the porous medium model, together with the experimental data, are presented in Fig. 3. The simulated curve after 7 d is slightly below the observed data. The 14-d curve is much lower than the data, indicating that strong bioreactions occurred within the porous model but not within the actual experiment. In this model, bacteria are present throughout the pore space (a 0.4 porosity). After the porous medium model was found unable to reproduce the experiment, the multiregion model as illustrated in Fig. 1 was chosen.
Multiregion Model
The physical parameters used for the three regions are given in Table 8
. For the first trial run, 15% of the bulk porosity (0.4) was assumed for the immobile region. In the multiregion simulations, a porosity of 1 was assigned for the mobile region and a porosity of 0.5 was assigned for the immobile region consisting of stagnant water and biomass. Changes in porosity caused by bacterial growth are not currently considered. A 0.05 porosity was assumed for the solid particle region. In "Interacting with the Solid Particle" below, a mineral composition is assigned for the solid particle region to illustrate the coupling of microbially mediated redox reactions with mineral dissolution and precipitation. The distance (Table 8, d = d1 + d2) between the mobile and immobile regions is assumed to be 2 x 10–5 m, a fraction of the sand particle sizes range from 1.252 x 10–4 to 2.52 x 10–4 m (von Gunten and Zobrist, 1993). The interficial area (A) was calibrated by matching the measured nitrate concentrations. The distance between the immobile and solid regions (3 x 10–5 m) was set slightly larger than the mobile–immobile regions. The same A was used for immobile and solid regions. In fact, diffusive flux between regions is proportional to D0
A/d, where D0 is a diffusion coefficient (10–9 m2/s) and is tortuosity. The between two regions is calculated according to values in two regions (see Table 8) by = (d11 + d22)/(d1 + d2). Values of d, D0, , and A all affect the diffusive flux. Therefore, errors in d, D0, and will be passed along to A because only the interficial area is calibrated. The other conditions and parameters are unchanged from the previous one-dimensional porous medium simulation.
The concentration profile of nitrate at different times obtained with the multiregion model is presented in Fig. 4
. Starting from a 15% immobile region, the 7-d curve agrees well the observed data. However, the 14-d curve is above the corresponding data (less reactive) but improved from the porous model. This may result from bacterial growth increasing volume of the immobile region. Currently, we simply let the run stop at 7 d and then restart by increasing the immobile region to 25% of the bulk porosity. We kept the diffusional length the same. Dynamic changes of the diffusional length together with immobile region volume and the interfacial area will be addressed in the future. The 14-d curve of 25% immobile region captures the observed data (Fig. 4). Only one data point at x = 0.0167 was off the curve, possibly because of other reasons (such as measurement errors). The simulated concentrations of dissolved organic carbon (DOC) and biomass obtained are presented in Fig. 5
. The DOC decreases slowly here and is not a limiting factor for the bacterial growth. The growth is very slow until 7 d; it then accelerates dramatically.
View larger version (17K):
[in this window]
[in a new window]
|
FIG. 4. Nitrate concentrations obtained with the multiregion model after 7 and 14 d, together with measured data.
|
|
View larger version (13K):
[in this window]
[in a new window]
|
FIG. 5. Concentrations of dissolved organic carbon (DOC) and biomass obtained with the multiregion model after 7 and 14 d.
|
|
The matches were adjusted with the interficial area A between mobile and immobile regions. The calibrated values are 38 m2 per cubic bulk medium for the initial period and 75 m2 for the late period. The match and parameter calibration suggest that TOUGHREACT is not only a useful interpretative tool for biogeochemical experiments but also can produce insight into processes and parameters of microscopic diffusion and their interplay with biogeochemical reactions.
Sulfate Reduction
The evolution of sulfate in the Column 2 experiment of von Gunten and Zobrist (1993) was simulated together with that of oxygen and nitrate. Organic matter is successively oxidized by oxygen, nitrate, and sulfate. Denitrification is inhibited by oxygen, whereas sulfate reduction is inhibited by both oxygen and nitrate. The multiregion model uses 25% of the bulk porosity (0.4) for the immobile region throughout the simulation. The increased immobile volume was attributed (i) to a longer simulation time (35 over 14 d) and (2) to more electron acceptors (3 over 1 d). An interfacial area of 75 m2 per cubic meter bulk medium, calibrated from the simulation of 7 to 14 d of Column 1 is used for the simulation of Column 2. The two columns have different injection flow rates (1.83 m/d for Column 1, and 0.37 m/d for Column 2). All other physical parameters are assumed to be the same as Column 1 because the sediment material is the same in both columns. Simulated concentrations of sulfate captured the measured data well (Fig. 6
). According to von Gunten and Zobrist (1993), the spatial resolution of the experiments did not allow for clear separation of oxygen respiration from denitrification. Virtually complete nitrate reduction took place within the first several centimeters of the flow distance. The simulation results are consistent with the nitrate observations.
View larger version (12K):
[in this window]
[in a new window]
|
FIG. 6. The simulated concentration profiles (lines) of sulfate, nitrate, and oxygen at steady state (35 d), together with measured data of sulfate along Column 2.
|
|
The simulated concentrations of DOC and biomass (in the immobile region) obtained are presented in Fig. 7
. The DOC concentrations decrease with the distance due to biodegradation. Bacteria grow significantly close to the inlet owing to oxygen respiration and denitrification. Then a second growth peak occurs at about 0.05 m due to the start of sulfate reduction.
Interacting with the Solid Particle
To illustrate the possible participation of minerals in the solid particle region, and the interplay between biodegradation and geochemical reactions, two reactive minerals, gypsum and goethite, were assumed to be initially present in the region. SO4 dissolved from gypsum diffuses to the immobile region and then is reduced to HS. In the solid region, Fe generated from goethite dissolution and HS diffused back from the immobile region may precipitate as pyrite. Ca from gypsum and HCO3 from organic matter biodegradation may form calcite. The later two minerals were specified as potential secondary minerals. A simple rate law was used for mineral dissolution and precipitation:
 | [18] |
where r is the kinetic rate (positive values indicate dissolution, and negative values precipitation), k is the rate constant (moles per unit mineral surface area and unit time), A is the specific reactive surface area per kilogram H2O, K is the equilibrium constant for the mineral–water reaction written for the destruction of one mole of mineral, and Q is the reaction quotient. The parameters used for the rate law Eq. [18] are presented in Table 9
. Rate constants for goethite and pyrite were increased by five orders of magnitude from the values in Palandri and Kharaka (2004), who compiled and fitted many of the experimental data reported by a large number of investigators. The rate increase was explicitly for taking into account bacterial catalysis. The dissolution-rate constants for gypsum and calcite were assigned to be two times those of goethite and pyrite. The assigned geochemical reaction rates could be comparable to sulfate bioreduction. Note that kinetic parameters, volume fractions and reactive surfaces in Table 9 may not reflect the field conditions, but the purpose of considering mineral dissolution and precipitation is to demonstrate the TOUGHREACT capabilities for geochemical processes. It should also be mentioned that the abundant iron in sediments is typically present as Fe(III), and it inhibits sulfate reduction. This inhibition effects is not considered here. However, in the pH range of natural waters, Fe(III) is extremely insoluble, and dissolved iron exists mainly as Fe(II) (von Gunten and Zobrist, 1993).
The biogeochemical simulation of the multiregion model is based on sulfate reduction of Column 2, using 25% of the bulk porosity for the immobile region and an interficial area of 75 m2 per cubic meter medium. An arbitrary increased initial biomass concentration of 60 mg/L is used to allow for significant biodegradation from the start. Initial concentrations for all other species are set equal to a very small value. Only DOC is assumed to be available in the injection water, with a concentration of 10 mg/L. The only electron acceptor SO4 results from the leaching of the materials in the solid region (gypsum dissolution).
Simulation results are presented in Fig. 8
and 9
. SO4 concentrations increase along the distance (Fig. 8a) because of gypsum dissolution (Fig. 9b; Ca concentrations take on the same pattern). Significant SO4 concentration gradients occurs between solid and immobile regions because of (i) consumption by sulfate reduction in the immobile region and (ii) larger distance of solid–immobile (3 x 10–5 over 2 x 10–5 m) and smaller (see Table 9). SO4 concentration in the immobile region is slightly higher than that in the mobile region. HS concentrations are higher in the immobile region due to sulfate reduction there (Fig. 8b), but its concentrations are reduced in the solid region resulting from the sink of pyrite precipitation (Fig. 9d). HCO3 concentrations increase (Fig. 8c) and DOC decrease (Fig. 8d) along the distance because of biodegradation. The difference of the latter two species among the three regions is very small.
View larger version (28K):
[in this window]
[in a new window]
|
FIG. 8. Simulated concentrations of dissolved components in the three regions after 13 d for the biogeochemical simulation. DOC = dissolved organic carbon.
|
|
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 9. (a) Simulated biomass concentrations in the immobile region and (b, c, and d) mineral dissolution and precipitation in the solid region at different times.
|
|
Bacterial growth rates over the distance (Fig. 9a) increase close to the inlet because of the increase in SO4 concentration, but later downstream decrease due to the low DOC concentrations. Biomass concentrations increase almost linearly with time. Returning to the solid region, pyrite precipitation pattern is similar to that of goethite dissolution. The former volume is slightly larger than the latter because of the slightly larger mole volume (23.94 over 20.82 cm3/mol). Calcite was not found in the simulation for the shorter column distance, but it could form along flow path if the distance was long enough because the HCO3 concentrations steadily increase (Fig. 8c). Gypsum dissolution promotes bio-reduction of sulfate. Then the biogeneration of HS reacts with Fe from goethite dissolution, resulting in pyrite precipitation. Geochemical processes in the solid region and microbiological processes in the immobile region influence each other via microdiffusion.
|
Conclusions
|
|---|
Aqueous and sorption kinetics and microbiological processes have been incorporated into TOUGHREACT. This added capability was tested by one-dimensional reactive transport simulation that involved kinetic biodegradation and sorption. Simulation results agree very well with those obtained by other simulators. The flexible space discretization, based on an integral finite difference approach, can use irregular gridding to model biogeologic structures. TOUGHREACT can deal directly with general multiregion and multiple interacting continua models for saturated and unsaturated porous and fractured media.
The applicability of this enhanced multiregion model for reactive transport of denitrification and sulfate reduction was evaluated by comparison with column experiments and model calibration. The matches with measured nitrate and sulfate concentrations were adjusted with the interfacial area between mobile and immobile regions. The values of 38 m2 per cubic bulk medium for the initial period and 75 m2 for the late period were calibrated. The match and parameter calibration suggest that TOUGHREACT not only can be a useful interpretative tool for biogeochemical experiments but also can produce insight into the processes and parameters of microscopic diffusion and their interplay with biogeochemical reactions.
The interaction of an immobile region with solid particles consisting initially of reactive minerals gypsum and goethite, was modeled. Gypsum dissolution promotes bioreduction of sulfate, after which the biogeneration of HS reacts with Fe from goethite dissolution, resulting in pyrite precipitation. It has been shown that geochemical processes in the solid region and microbiological processes in the immobile region influence each other via microdiffusion.
The resulting biogeochemical-transport simulation capabilities may be useful for many subsurface problems, including acidic mine drainage remediation, organic matter decomposition, oil and gas maturation, sulfite reduction in oil fields, and effective environmental remediation of groundwater contamination. However, evaluation of the model against experimental data at larger than a column, and in multidimensional flow fields, is recommended before applying to field-scale problems.
|
ACKNOWLEDGMENTS
|
|---|
The author thanks Nicolas Spycher and Karsten Pruess for reviews of the manuscript and constructive discussions. The author appreciates Guoxiang Zhang for providing the numerical data from the Bio-CORE2D simulation. This work was supported by the Laboratory Directed Research and Development Program of the Ernest Orlando Lawrence Berkeley National Laboratory, under U.S. Department of Energy Contract No. DE-AC02-05CH11231.
|
REFERENCES
|
|---|
- Clement, T.P., B.S. Hooker, and R.S. Skeen. 1996. Macroscopic models for predicting changes in saturated porous media properties caused by microbial growth. Ground Water 34(5):934–942.[CrossRef][Web of Science]
- Curtis, G.P. 2003. Comparison of approaches for simulating reactive solute transport involving organic degradation reactions by multiple terminal electron acceptors. Comput. Geosci. 29:319–329.[CrossRef]
- Doussan, C., G. Poitevin, E. Ledoux, and M. Detay. 1997. River bank filtration: Modelling of the changes in water chemistry with emphasis on nitrogen species. J. Contam. Hydrol. 25:129–156.[CrossRef]
- Gwo, J.P., P.M. Jardine, G.V. Wilson, and G.T. Yeh. 1995. A multiple-pore-region concept to modeling mass transfer in subsurface media. J. Hydrol. 164:217–237.[CrossRef]
- Harvey, R.W., L.H. George, R.L. Smith, and D.R. Leblanc. 1989. Transport of microspheres and indigenous bacteria through a sandy aquifer: Results of a natural and forced gradient tracer experiment. Environ. Sci. Technol. 23:51–56.
- Narasimhan, T.N., and P.A. Witherspoon. 1976. An integrated finite difference method for analyzing fluid flow in porous media. Water Resour. Res. 12:57–64.[CrossRef]
- Palandri, J., and Y.K. Kharaka. 2004. A compilation of rate parameters of water–mineral interaction kinetics for application to geochemical modeling. USGS Open File Report 2004-1068. USGS, Denver, CO.
- 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, and inverse geochemical calculations, USGS Water-Resources Investigations Report 99-4259. USGS, Denver, CO.
- Pruess, K., and T.N. Narasimhan. 1985. A practical method for modeling fluid and heat flow in fractured porous media. Soc. Petrol. Eng. J. 25:14–26.
- Pruess, K., C. Oldenburg, and G. Moridis. 1999. TOUGH2 user's guide, version 2.0. LBL-43134. Lawrence Berkeley Lab., Berkeley, CA.
- imnek, J., and D.L. Suarez. 1994. Two-dimensional transport model for variably saturated porous media with major ion chemistry. Water Resour. Res. 30:1115–1133.[CrossRef]
- imnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272:14–35.[CrossRef]
- Sonnenthal, E., A. Ito, N. Spycher, M. Yui, J. Apps, Y. Sugita, M. Conrad, and S. Kawakami. 2005. Approaches to modeling coupled thermal, hydrological, and chemical processes in the Drift Scale Heater Test at Yucca Mountain. Int. J. Rock Mechanics Mining Sci. 42:698–719.[CrossRef]
- Spycher, N.F., E.L. Sonnenthal, and J.A. Apps. 2003. Fluid flow and reactive transport around potential nuclear waste emplacement tunnels at Yucca Mountain, Nevada. J. Contam. Hydrol. 62–63:653–673.
- Steefel, C.I. 1997. Incorporating intra-aqueous kinetics into reactive transport models. Paper presented at Reactive Transport Workshop, Richland, WA. October 2007.
- Steefel, C.I., and K.T.B. MacQuarrie. 1996. Approaches to modeling of reactive transport in porous media. p. 83–129. In P.C. Lichtner, C.I. Steefel, and E.H. Oelkers (ed.) Reactive transport in porous media. Reviews in Mineralogy 34. Mineralogical Soc. of America, Chantilly, VA.
- Tebes-Stevens, C., A.J. Valocchi, J.M. van Briesen, and B.E. Rittmann. 1998. Multicomponent transport with coupled geochemical microbiological reactions: Model description and example simulations. J. Hydrol. 209:8–26.[CrossRef]
- Thullner, M., M.H. Schroth, J. Zeyer, and W. Kinzelbach. 2004. Modeling of a microbial growth experiment with bioclogging in a two-dimensional saturated porous media flow field. J. Contam. Hydrol. 70:37–62.[CrossRef][Web of Science][Medline]
- van Beelen, P., and K. Fleuren Kemila. 1989. Enumeration of anaerobic and oligotrophic bacteria in subsoils and sediments. J. Contam. Hydrol. 4:275–284.[CrossRef]
- von Gunten, U., and J. Zobrist. 1993. Biogeochemical changes in groundwater infiltration systems: Column studies. Geochim. Cosmochim. Acta 57:3895–3906.[CrossRef]
- Wunderly, M.D., D.W. Blowes, E.O. Frind, and C.J. Ptacek. 1996. Sulfide mineral oxidation and subsequent reactive transport of oxidation products in mine tailings impoundments: A numerical model. Water Resour. Res. 32:3173–3187.[CrossRef]
- Xu, T., K. Pruess, and G. Brimhall. 1999. An improved equilibrium-kinetics speciation algorithm for redox reactions in variably saturated flow systems. Comput. Geosci. 25(6):655–666.[CrossRef]
- Xu, T., and K. Pruess. 2001. Modeling multiphase non-isothermal fluid flow and reactive geochemical transport in variably saturated fractured rocks: 1. Methodology. Am. J. Sci. 301:16–33.
- Xu, T., E.L. Sonnenthal, N. Spycher, and K. Pruess. 2006. TOUGHREACT: A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media. Comput. Geosci. 32(2):145–165.[CrossRef]
- Yeh, G.T., and V.S. Tripathi. 1991. A model for simulating transport of reactive multispecies components: Model development and demonstration. Water Resour. Res. 27:3075–3094.[CrossRef]
- Zhang, G. 2001. Non-isothermal hydrobiogeochemical models in porous media. Ph.D. diss., Univ. of La Coruña, Spain.
- Zysset, A., E. Stauffer, and T. Dracos. 1994. Modelling of reactive groundwater transport governed by biodegradation. Water Resour. Res. 30:2423–2434.[CrossRef]
This article has been cited by other articles:
|
|
|
|
 
S. Finsterle, C. Doughty, M.B. Kowalsky, G. J. Moridis, L. Pan, T. Xu, Y. Zhang, and K. Pruess
Advanced Vadose Zone Simulations Using TOUGH
Vadose Zone J.,
May 27, 2008;
7(2):
601 - 609.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
H.-H. Liu and T. H. Illangasekare
Preface: Recent Advances in Modeling Multiphase Flow and Transport with the TOUGH Family of Codes
Vadose Zone J.,
February 25, 2008;
7(1):
284 - 286.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
