CHROTRAN: A mathematical and computational model for in situ heavy metal remediation in heterogeneous aquifers
Scott K. Hansen, Sachin Pandey, Satish Karra, Velimir V. Vesselinov
CCHROTRAN : A mathematical and computational model for in situheavy metal remediation in heterogeneous aquifers ∗ Scott K. Hansen †1 , Sachin Pandey , Satish Karra , and Velimir V. Vesselinov Computational Earth Science Group (EES-16), Los Alamos National LaboratoryMarch 7, 2017
Abstract
Groundwater contamination by heavy metals is a critical environmental problem for which in situ remediationis frequently the only viable treatment option. For such interventions, a three-dimensional reactive transport modelof relevant biogeochemical processes is invaluable. To this end, we developed a model,
CHROTRAN , for in situtreatment, which includes full dynamics for five species: a heavy metal to be remediated, an electron donor, biomass,a nontoxic conservative bio-inhibitor, and a biocide. Direct abiotic reduction by donor-metal interaction as well asdonor-driven biomass growth and bio-reduction are modeled, along with crucial processes such as donor sorption,bio-fouling and biomass death. Our software implementation handles heterogeneous flow fields, arbitrarily manychemical species and amendment injection points, and features full coupling between flow and reactive transport.We describe installation and usage and present two example simulations demonstrating its unique capabilities. Onesimulation suggests an unorthodox approach to remediation of Cr(VI) contamination.
Keywords : in situ remediation, abiotic reduction, bioremediation model, subsurface flow, reactive transport ∗ Los Alamos National Laboratory technical report:
LA-UR-16-29041 † Corresponding author. [email protected] a r X i v : . [ q - b i o . O T ] M a r Introduction
Heavy metals, including chromium, arsenic, copper, nickel, selenium, technetium, uranium, and zinc, are widespreadand hazardous subsurface contaminants in groundwater aquifers (Appelo and Postma, 2004; Tchounwou et al., 2012).For many heavy metals, their most stable oxidation state is often the most toxic (Duruibe et al., 2007; Hashim et al.,2011), and this oxidation state is typically the highest that occurs under near-surface conditions. Additionally, thechemical reduction of certain metals is known to reduce their mobility (Violante et al., 2010). This has inspired effortsto manipulate in situ conditions to stimulate microbial growth and achieve biologically mediated metals reduction.This technique has been demonstrated, at least in some settings, for chromium, uranium and selenium (Lovley, 1993,1995), nickel (Zhan et al., 2012), technetium (Istok et al., 2004), and copper (Andreazza et al., 2010), and has beennoted as a viable bioremediation technique by recent critical reviews (Hashim et al., 2011; Wu et al., 2010). Biopre-cipitation, a process by which microbiological exudates react with metals to produce an insoluble compound, has beenwidely observed (Malik, 2004; Van Roy et al., 2006; Radhika et al., 2006) and has been noted by Wu et al. (2010) asa remediation method. Bio-stimulants have also recently been shown to effectively reduce chromium through abioticredox pathways (Chen et al., 2015; Hansen et al., 2016), and after fermentation for other metals (Hashim et al., 2011).Naturally, designing a remedial intervention using one of this family of techniques benefits greatly from the use of amulti-dimensional/multi-component numerical model of groundwater flow, contaminant transport, and biogeochemi-cal processes to evaluate different remediation strategies under varying field conditions. The model should be capableof capturing the transport behavior of electron donors, biomass, and other species, dominant biogeochemical reactions,and how these processes influence and are influenced by subsurface flow.Although the development on in situ bioreactive transport models goes back to at least the 1980s, the literature isnot vast. Early work focused on in situ bioremediation of toxic organic compounds through oxidation. A thoroughmathematical and 2D numerical study representative of this approach is due to Chiang et al. (1991), who presented athree-equation model involving a mobile electron donor (assumed to be the contaminant), mobile dissolved oxygen,and immobile biomass. The contaminant was assumed to be consumed only in the microbial growth reaction, whichwas linear in biomass, Monod (Monod, 1949) in electron donor, and Monod in electron acceptor. Wheeler et al. (1992)subsequently extended a reactive model of this sort to three dimensions to simulate biodegradation of CH . Travis(1993) presented a more complicated, unsaturated three-dimensional model, which introduced Monod dependence onnutrients, and the potential for two electron donors, with one inhibiting the other. This approach was further elabo-rated upon in a study of TCE degradation (Travis and Rosenberg, 1998) by accounting for living and dead microbes,microbial predators, and first-order kinetic sorption of all aqueous species (microbes were treated as mobile). Anothercomplex oxidation model was developed by Suk et al. (2000), which explicitly modeled both mobile and immobilebiomass, contained a decay network, and featured both anaerobic and aerobic oxidation, in competition.The development of models for metal reduction is comparatively more recent. For U(VI), field-scale modeling studieshave been performed on bio-reduction under anaerobic conditions at the Old Rifle Site in Colorado (Li et al., 2010,2011; Yabusaki et al., 2011). These conceptions treat the contaminant as the sole electron acceptor, with an externallyapplied electron donor, and the implied equations have a similar form to those devised by Chiang et al. (1991): linearin biomass, Monod in contaminant, and Monod in electron donor. For clarity, this is expressed symbolically as: ∂ C ∂ t ∝ ∂ D ∂ t ∝ B CK C + C DK D + D , (1)where C , is the U(VI) concentration, B is the biomass concentration, D is the electron donor concentration, K C is themetal reduction Monod constant, and K D is the electron donor Monod constant. K C and K D respectively representthe concentration of C and D at which the reaction rate is halved. Recently, Molins et al. (2015) have published anumerical study of a column experiment with multiple species, all of whose dynamics are of the above form, butincluding an extra chemical inhibition factor. The models of Li et al. (2010, 2011) were implemented at field-scale inCrunchFlow (Steefel et al., 2015), using its capability to represent single and multiple Monod formulations.Systems of governing reactive transport equations for enzymatic microbial Cr(VI) reduction have been presented byAlam (2004), and by Shashidhar et al. (2007). Shashidhar et al. (2007) described the Cr(VI) degradation reactionslightly differently from Li et al. (2011): ∂ C ∂ t ∝ ∂ D ∂ t ∝ B K C (cid:48) K C (cid:48) + C DK D + D . (2)2 C (cid:48) is the concentration of Cr(VI) at which the reaction rate is halved, which is similar to K C . However, although(2) appears superficially similar to (1), the C factor represents entirely different behavior: not as an energy source butrather as an inhibitor. Interestingly, since the RHS of (2) is a proxy for the biomass growth reaction, C consumptionis modeled as proportional to biomass growth, but the biomass growth rate is modeled as independent of C . Biomassdynamics were governed by a growth term proportional to donor consumption and a first-order decay term, accountingfor eventual biomass die-off. Other authors (e.g., Somasundaram et al., 2009) have used a similar approach. Alam(2004) presented a relatively complex model which included transport with both mobile and immobile biomass, andalso included two enzymes (both created due to biomass growth, but one conserved, and one irreversibly consumedduring bio-reduction). Neglecting the irreversibly-consumed enzyme and the mobile-immobile behavior, this modelshares its electron donor and biomass dynamics with the model of Shashidhar et al. (2007). It differs significantly fromother models that we are aware of by treating the Cr(VI) degradation reaction in this model as an incidental enzymaticprocess, and is governed by the following Monod equation: ∂ C ∂ t ∝ B CK C + C . (3)There is strong experimental support for this approach (e.g., Okeke, 2008), and this is arguably more defensible in areal complex geochemical system in which there are multiple competing donors and receptors, and given that thereis evidence for indirect reduction pathways, e.g., by metabolites (Priester et al., 2006). All of the models of Cr(VI)bio-reduction, discussed above, appear to be one-dimensional only.Our literature review did not reveal discussion of field-scale bio-reduction models for other heavy metal species. It thusappears that the primary example of a bio-reduction model applicable to modeling a real-world remediation scheme isthe CrunchFlow model of uranium treatment at the Rifle site, which was discussed above. We set out to develop a newmodel, dubbed CHROTRAN , which is optimized for modeling bioremediation of Cr(VI), but of sufficient generalitythat it might be used for bioremediation of other metals, or for abiotic reduction, with ease. The key features of themodel we developed are as follow:
Direct abiotic reaction between electron donor and contaminant
Recent experimental results (Chen et al., 2015;Hansen et al., 2016) have established a rapid direct redox reaction when molasses is used as an electron donorand Cr(VI) is the contaminant, rather than the bio-mediated reaction previously posited. It is thus crucial toinclude this behavior in a model aimed at remediation design.
Indirect Monod kinetics
On account of the evidence (Wang and Xiao, 1995; Okeke, 2008; Hansen et al., 2016) formodeling Cr(VI) degradation with (3), we implemented this general formulation as opposed to one which tiesall contaminant degradation to a single biomass growth equation.
Bio-fouling / Bio-clogging
It is well known in practice that one of the problems afflicting bioremediation schemesis build-up of biological material near the amendment injection point. This reduces the hydraulic conductivity,interfering with amendment injection, and may rapidly consume any amendment that does manage to passthrough it. The model thus contains feedback between local biomass concentration and flow parameters such asporosity and hydraulic conductivity.
Biomass crowding
Similarly, if biomass becomes overly dense, this causes cell stress, which reduces the rate offurther growth. Since clogging is enabled, this behavior was added as well.
Modeling of amendment additives
To address clogging or to attempt to spread electron donors farther from the wellbefore they are consumed, additional chemicals may be injected to reduce biomass concentrations, and theirreactive transport behavior is incorporated.
Multiple donor consumption pathways
The best model of electron donor consumption by biomass may be propor-tional to biomass concentration or biomass growth, and the model can handle any such combination.Building this functionality required custom programming beyond what is embedded in existing reactive transportcodes (Steefel et al., 2015). To accomplish our goal, we turned to
PFLOTRAN (Lichtner et al., 2015), which is opensource, has a modular structure, and a “reaction sandbox” interface (Hammond, 2015) that allows derivative versionswith custom reaction behavior to be developed and compiled. In this fashion, no changes to the flow and transportpart of
PFLOTRAN are needed. We developed
CHROTRAN based on the existing
PFLOTRAN code framework, taking3dvantage of the reaction sandbox interface to implement complex model features not included in basic microbialpackages and leverage other aspects of
PFLOTRAN , such as its high-performance computing capabilities.In Section 2, we present the mathematical details of
CHROTRAN and justify some of the decisions underlying themodel. In Section 3, we discuss how to install and use the software. In Section 4, we present two numerical studieswhich illustrate
CHROTRAN and also suggest an interesting conclusion regarding Cr(VI) remediation. In Section 5, webriefly summarize what has been presented.
We consider flow and transport at aquifer scale. Conceptually, the aquifer is modeled as saturated, with incompressiblewater moving in accordance with Darcy’s law (
CHROTRAN can also simulate partially-saturated vadose-zone flow andtransport). Two transport processes are considered, namely, advection with the Darcy flow and Fickian dispersion.Multiple reaction terms are then added in order to capture the complex chemical dynamics during remediation. As themodel is intended to be used for remedial design, every effort was made to simplify the formulation to use the smallestnumber of explanatory variables and parameters, and to keep the equations at a high level of abstraction, so they arenot tied to one particular set of chemical species.The following are the several species whose dynamics are captured by the system of reaction equations, each withtheir own symbols:
Biomass, B [ mol L − ] , representing the concentration of all microbes and their associated extracellular material. Thequantification of biomass as a “molar” rather than a mass concentration is unusual, and was done for two reasons:(i) to avoid hard-coding units in which biomass concentration is to be specified, and (ii) to simplify presentationof the model, so all governing equations have the same units. A mole of biomass should be understood as anequivalent mass: any quantity can be used, as long as one uses a consistent definition throughout the model. Inthe examples in this paper, we use the definition 1 mol ≡ Aqueous contaminant, C [ mol L − ] , which we here assume is a heavy metal ion in its oxidized state, such as Cr(VI)or U(VI). Electron donor, which is part of the chemical amendment, and may be1. immobile, represented by D i [ mol L − ] , or2. mobile, represented by D m [ mol L − ] ,with exchange of mass between the two states. Nonlethal biomass-growth inhibitor, I [ mol L − ] , such as ethanol, which is modeled as a conservative species butacts to slow microbial growth. Biocide, X [ mol L − ] , which reacts directly with biomass and is consumed.For convenience, we also define a total species aqueous concentration of the electron donor, D , according to theformula D = D i θ ( x , t ) + D m [ mol L − ] , where θ ( x , t ) [-] is the current porosity at x . For simplicity, we assume that boththe mobile and immobile donor participate equally in all reactions. Flow may be modeled using the balance of water mass given by ddt ( ρ w θ ) + ∇ · q = q M ( x , t ) , (4)4ith the water mass fluxes related to head via Darcy’s law: q = − ∇ ( ρ w K ( x, t ) h ( x , t )) , (5)where K ( x ) [ m / s ] is the local hydraulic conductivity and h ( x , t ) [ m ] is the local hydraulic head, θ is the porosityand q M ( x ) [ kg m − s − ] is the local mass injection rate into the system and ρ w [ kg m − ] is the density of water. Wenote that, since CHROTRAN is built on top of
PFLOTRAN , it inherits all of
PFLOTRAN ’s groundwater flow modelingcapabilities. This includes the ability to consider unsaturated and otherwise multiphase flow conditions, which areout-of-scope for the present discussion. Please see the
PFLOTRAN user manual (Lichtner et al., 2015) for details on itscomplete capabilities.The hydraulic conductivity is continually updated in accord with the relation K ( x , t ) = K ( x , ) θ ( x , t ) θ , (6)where θ [-] is the spatially-uniform initial porosity, and θ ( x , t ) is calculated according to θ ( x , t ) = θ − B ( x , t ) ρ B , (7)where ρ B [ mol L − ] is the intrinsic biomass density. (Note that, using our proposed definition of 1 mol of biomass as1 g of biomass, 1 mol L − = − .) We define T {·} to be an advective-dispersive transport operator, which characterizes the hydrodynamic effects onsolute transport. For c , the concentration of an arbitrary mobile species, T { c } ≡ − q · ∇ c + ∇ · ( θ D ( q ) ∇ c ) , (8)where D is a dispersion tensor that depends on the longitudinal and transverse dispersivities, molecular diffusion aswell as the Darcy flux. For the work in this paper, we will only consider molecular diffusion and thereby we set D = D m I . Note that, while this is not shown explicitly for compactness, all symbols in this equation are functions of x and t . We define one governing equation for each species, mobile or immobile, as well as two equations defining reactionrate expressions for algebraic convenience. The equations involve numerous parameters, whose symbols, units, andlong-form name in the
CHROTRAN input file are summarized in Table 1. The parameter symbols follow a schemein which the first letter encodes the physical interpretation of the parameter and the subscript specifies the governingequation in which they participate. A symbol beginning with Γ is a second-order mass action rate constant, withunits [ L mol − s ] . A symbol beginning with K is Monod or inhibition constant with units of concentration, mol L − or mol L − , and represents the concentration at which a process rate becomes 50% of its maximum rate, all otherparameters being equal. A symbol beginning with λ has units of s − and is interpreted as a pure first-order reactionrate constant. A symbol beginning with S is dimensionless, and represents a stoichiometric relationship between areaction rate and the consumption rate of a certain species. Before presenting the equations, it is useful to review allof the chemical processes that are incorporated into the model: Abiotic reduction
This is an aqueous-phase bimolecular reaction between the electron donor, D , and the contaminant, C . It is modeled with a classical second-order mass action rate law. Bio-reduction
This represents the removal of the contaminant, C by the biomass, B . The process is linear in B , andMonod in C , which is a common assumption for bio-mediated processes. Note that we are not assuming thatreduction of C is directly tied to any particular cell metabolic process. This form is sufficiently general that itcan capture other bio-remediation processes besides bio-reduction of heavy metals.5 iocide reaction This is an inter-phase bimolecular reaction between the biocide, X , and the biomass, B . It is modeledwith a classical second-order mass action rate law, with the added condition that B cannot fall below a specifiedminimum concentration B min . Biomass growth
The core biomass growth reaction irreversibly consumes electron donor, D to increase biomass, B . As a biologically catalyzed reaction, it is assumed to be linear in B and Monod in C . Two inhibitioneffects are assumed: a biomass crowding term, tunable with exponent α , attenuates growth rate as the biomassconcentration rises. The nonlethal inhibitor concentration, I , also reduces the reaction rate as its concentrationincreases. Mobile-immobile mass transfer (MIMT)
This is a process with first-order kinetics, which models sorptive retarda-tion of the electron donor.
Natural decay
This is an empirical process reflecting the idea that, if left unstimulated, both the amount of livingcells and the amount of extracellular material in the aquifer will ultimately return to their natural backgroundlevel (i.e. B min ). This is modeled as a first-order process. Respiration
This represents consumption of the electron donor for purposes of life maintenance, unrelated to biomassgrowth. This is described by a first-order rate law which is proportional to biomass concentration, B .The explanations of the operative processes and of parameter interpretation above help the descriptions of factors andterms in the governing equations presented below. The biomass growth reaction is linear in biomass concentration, has a Monod dependency on electron donor, a tunableinhibition factor due to biomass crowding, and a classic inhibition factor describing the impact of the nonlethal growthinhibitor (as indicated by comment braces): µ B = λ B B e − donor (cid:122) (cid:125)(cid:124) (cid:123) DK D + D (cid:18) K B K B + B (cid:19) α (cid:124) (cid:123)(cid:122) (cid:125) crowding inhibition (cid:122) (cid:125)(cid:124) (cid:123) K I K I + I (cid:20) molL b s (cid:21) . (9)The direct, abiotic reduction reaction is represented by a classic, second-order mass action law: µ CD = Γ CD CD (cid:20) molL b s (cid:21) . (10) The mobile components are all governed by the advection-dispersion operator, T , defined previously, and also affectedby extra terms implementing the chemical processes outlined earlier (as indicated by comment braces):6 θ C ∂ t = T { C } − bio − reduction (cid:122) (cid:125)(cid:124) (cid:123) λ C B CK C + C − S C µ CD (cid:124) (cid:123)(cid:122) (cid:125) abiotic reduction (cid:20) molL b s (cid:21) , (11) ∂θ D m ∂ t = T { D m } − biomass growth (cid:122) (cid:125)(cid:124) (cid:123) S D D m D µ B − λ D D m D B (cid:124) (cid:123)(cid:122) (cid:125) respiration − abiotic reduction (cid:122) (cid:125)(cid:124) (cid:123) S D µ CD D m D − λ D i θ D m + λ D m D i (cid:124) (cid:123)(cid:122) (cid:125) MIMT (cid:20) molL b s (cid:21) , (12) ∂θ I ∂ t = T { I } (cid:20) molL b s (cid:21) , (13) ∂θ X ∂ t = T { X } − biocide reaction (cid:122) (cid:125)(cid:124) (cid:123) Γ X BX (cid:20) molL b s (cid:21) . (14) The immobile component concentrations are affected only by the reactive processes outlined above: ∂ B ∂ t = biomass growth (cid:122)(cid:125)(cid:124)(cid:123) µ B − λ B ( B − B min ) (cid:124) (cid:123)(cid:122) (cid:125) natural decay − biocide reaction (cid:122) (cid:125)(cid:124) (cid:123) Γ B ( B − B min ) X (cid:20) molL b s (cid:21) , (15) ∂ D i ∂ t = − biomass growth (cid:122) (cid:125)(cid:124) (cid:123) S D D i D µ B − λ D D i D B (cid:124) (cid:123)(cid:122) (cid:125) respiration − abiotic reduction (cid:122) (cid:125)(cid:124) (cid:123) S D µ CD D i D + λ D i θ D m − λ D m D i (cid:124) (cid:123)(cid:122) (cid:125) MIMT (cid:20) molL b s (cid:21) . (16) CHROTRAN currently is undergoing an open source licensing process, and will be available as soon as possible asfreely downloadable Fortran 2003 source code. Until this process is complete, the software may be made available ona restricted basis by arrangement with the authors.
CHROTRAN must be compiled using the GFortran complier (freely available as part of the GNU Compiler Collec-tion). It is based on the open-source
PFLOTRAN code base, and the installation procedure is the essentially the sameas that required to build
PFLOTRAN from source, and
CHROTRAN requires all the libraries upon which
PFLOTRAN depends, including PETSc (Balay et al., 2016) and others. For installation of required libraries, the
PFLOTRAN in-stallation instructions are applicable, except that CHROTRAN , rather than
PFLOTRAN , should be cloned from itsrepository once all the dependencies have installed. To build CHROTRAN itself, navigate to
CHROTRAN executable will be called chrotran .) Available at http://documentation.pflotran.org/user_guide/how_to/installation/installation.html . Not officially released yet; anonymous access can be facilitated through editor for review purposes.
CHEMISTRY card (Figure 1) and the mathematical symbolsshown in Section 2. Symbol Units Name in
CHEMISTRY card α - EXPONENT B B min mol L − BACKGROUND CONC B Γ B L mol − s − MASS ACTION B Γ CD L mol − s − MASS ACTION CD Γ X L mol − s − MASS ACTION X K B mol L − INHIBITION B K C mol L − INHIBITION C K D mol L − MONOD D K I mol L − INHIBITION I λ B s − RATE B 1 λ B s − RATE B 2 λ C s − RATE C λ D s − RATE D λ D i s − RATE D IMMOB λ D m s − RATE D MOBIL ρ B mol L − DENSITY B S C - STOICHIOMETRIC C S D - STOICHIOMETRIC D 1 S D - STOICHIOMETRIC D 2 A CHROTRAN input file is of the same format as a
PFLOTRAN input file. Information on how to set up such a file isavailable in the
PFLOTRAN user manual (Lichtner et al., 2015). However, to use
CHROTRAN ’s additional functionality,a few of the input cards (top-level blocks, in
PFLOTRAN jargon) must contain some particular content. The required
CHEMISTRY card format is shown in Figure 1, with bold text being mandatory and standard-weight text being user-alterable. The required
SIMULATION , MATERIAL PROPERTY , and (initial)
CONSTRAINT card formats are shown in Figure2, again with bold text being mandatory and standard-weight text being user-alterable. Comments in the input file arepreceded by the character .In addition to these cards being properly formatted, there must exist a chemistry database at the (absolute or relative)path specified after the
DATABASE keyword in the
CHEMISTRY card, and it must, at a minimum contain the lines shownin Figure 3. The one exception to bold text being mandatory is that species names can be changed at will, as longas there is consistency between the
CHEMISTRY card and the chemistry database. For instance, one could change allinstances of the text
Cr(VI) in both of those locations to
U(VI) or all instances of the text chubbite to etibbuhc ,with no alteration in execution behavior (besides, obviously, the species names used in the output files).The chemistry database contains lines for five mobile species: water, plus the mobile species in the CHROTRAN kineticslisted in Section 2.3: C , D m , I , and X . The database also contains a line for a “dummy” mineral species, chubbite ,which does not correspond to any species previously mentioned. This species is treated as a mineral which is specifiedas inactive with respect precipitation/dissolution by setting its kinetic rate constant ( RATE CONSTANT ) to zero. Themineral is included as a surrogate for biomass and porous media volume in
CHROTRAN and is updated according toEquation 7 to track 1 − θ ( x , t ) . The initial volume fraction of chubbite thus defines the initial porosity. The formatof a chemistry database is discussed in more detail in the PFLOTRAN user manual.Once you have saved your input file as, e.g. test.in , it is easy to run the code from the console. Navigate to
CHROTRAN can save flow field velocities, concentrations of all species, permeabilities, and porosities at8ny specified times in an .h5 format file. This file format can be visualized natively using freely-available standalonetools such as VisIt and ParaView, and are also accessible from Python scripts by means of the h5py library and fromJulia scripts by means of the
HDF5 package.
The
PFLOTRAN software from which
CHROTRAN derives its numerical flow and reactive transport solvers has gonethrough extensive quality assurance testing, has been benchmarked against other reactive transport solvers, and isused inside and outside the U.S. Department of Energy for mission-critical analytical work. The new bio-reactivetransport model that is constitutive of
CHROTRAN is not available in any other software, so direct benchmarking is notpossible. However, extensive quality testing has been performed by the developers. We have validated through batchand multidimensional simulations that
CHROTRAN does satisfy the governing equations we present for chemistry andpermeability, and also that it gives plausible, physically consistent results for a wide range of scenarios. CHROTRAN remediation case studies
To demonstrate the capabilities of our software, we present two example studies, which together illustrate the interac-tions of all the types of chemical species it permits to be modeled, along with its treatment of bio-clogging. The inputfiles for these two examples can be found in the chrotran examples directory of the
CHROTRAN repository.
This study concerns the co-injection of molasses (electron donor, D ) and ethanol (nonlethal bio-inhibitor, I ) into asingle well drilled in a heterogeneous aquifer with an appreciable background Cr(VI) concentration. The competitionbetween direct abiotic reduction of Cr(VI) by molasses and bio-reduction of Cr(VI), which exists since both reductionpathways consume the electron donor, along with the impact of suppressing the biomass growth is explored. The basicparameters used are those shown in Figures 1 and 2, with changes as indicated below.Four related simulations are performed on the same 100 ×
100 m two-dimensional heterogeneous hydraulic con-ductivity field, with geometric mean conductivity K g = − m s − , a multi-Gaussian correlation structure with ex-ponential semivariogram with correlation length of 4 m and σ K =
2. Each simulation takes place over a span of500 days and begins with ε = − initial concentrations of all species, except B ( x , ) = B min = − mol L − and C = . × − mol L − (1000 ppb Cr(VI)). In all cases, there are no flow boundaries at the north and south ofthe domain ( y = y =
100 m), and constant head boundaries are imposed at the west and east of the domain( x = x =
100 m) such that there is a drop of head of 0.28 m between these faces. A single injection wellexists at ( x , y ) = (
25 m ,
50 m ) . For the first 10 days of the simulation, there is no injection into the well. From 10 d to30 d, injection is performed at the well with constant volumetric flow rate 272.55 m d − with species concentrationsdiscussed below. From 30 d to 500 d, there is again no injection at the well. A very large (arbitrary) ρ B is assumed, soas to eliminate the effect of biomass clogging from this simulation.The four simulations differ in their chemistry only. Two direct abiotic reduction rates are considered: Γ CD = − s − and Γ CD = − s − , as are two different ethanol concentrations in the injection fluid: I = − and I = ε mol L − , in all four possible combinations. The injection fluid chemistry always has Cr(VI) concentration equalto the initial concentration ( C = . × − mol L − ), ensuring that no chromium disappearance is due to dilution,and molasses concentration D = × − mol L − .Concentrations of Cr(VI) for each scenario are shown in Figure 4. It is apparent that little persistent reduction dueto biomass alone occurs, although ethanol co-injection does increase biomass footprint, which has a noticeable andpersistent effect. By contrast, the rapid abiotic reaction between Cr(VI) and a constituent of molasses has more impact.This is attributable to the fact that molasses has a large reducing capacity, background concentrations of Cr(VI) arerelatively low, and it has a retardation factor of around 150 (obtained from Shashidhar et al. (2006)), meaning that it hasthe potential to form a persistent permeable reactive barrier around the well. The better performance in the presence of9 HEMISTRY PRIMARY_SPECIES molasses Cr(VI) ethanol biocide END IMMOBILE_SPECIES biomass molasses_im END MINERALS chubbite
END REACTION_SANDBOX CHROTRAN_PARAMETERS NAME_D_MOBILE molasses NAME_D_IMMOBILE molasses_im NAME_C Cr(VI) NAME_B biomass NAME_I ethanol NAME_X biocide NAME_BIOMINERAL chubbite EXPONENT_B
BACKGROUND_CONC_B
MASS_ACTION_B
MASS_ACTION_CD
MASS_ACTION_X
RATE_B_1
RATE_B_2
RATE_C
RATE_D
RATE_D_IMMOB
RATE_D_MOBIL
INHIBITION_B
INHIBITION_C
MONOD_D
INHIBITION_I
DENSITY_B
STOICHIOMETRIC_C
STOICHIOMETRIC_D_1
STOICHIOMETRIC_D_2
END END MINERAL_KINETICS chubbite RATE_CONSTANT 0.d0 END END UPDATE_POROSITY MINIMUM_POROSITY 1.d-4 UPDATE_PERMEABILITY UPDATE_MINERAL_SURFACE_AREA DATABASE ./chem.dat
OUTPUT ALL FREE_ION TOTAL END
LOG_FORMULATION
END
Figure 1: Example
CHEMISTRY card for
CHROTRAN input file. Bold text should not be altered. However, additionalspecies may be added to the
PRIMARY SPECIES , IMMOBILE SPECIES , MINERALS , and
MINERAL KINETICS blocks, ifdesired. Additional sandboxes can also be used in the
REACTION SANDBOX block.10
IMULATION SIMULATION_TYPE SUBSURFACE PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE RICHARDS
END SUBSURFACE_TRANSPORT transport
GLOBAL_IMPLICIT NUMERICAL_JACOBIAN END END END MATERIAL_PROPERTY soil1 ID
1 TORTUOSITY 0.1d0
PERMEABILITY
DATASET Permeability
END PERMEABILITY_POWER 1.0 PERMEABILITY_CRITICAL_POROSITY 0.0
PERMEABILITY_MIN_SCALE_FACTOR 1.d-4
CHARACTERISTIC_CURVES cc1
END
CONSTRAINT initial
CONCENTRATIONS molasses ethanol biocide
Cr(VI)
END
IMMOBILE biomass molasses_im
END
MINERALS chubbite
END END
Figure 2: Additional cards that require particular content in order for
CHROTRAN to work properly. In the
SIMULATION card, the
NUMERICAL JACOBIAN option must be specified. In the
MATERIAL PROPERTY card, the OP-TION
PERMEABILITY MIN SCALE FACTOR 1.d4 option should be set. In
CONSTRAINT cards, species that are notpresent should have small, but non-zero concentrations assigned. The concentration of
NAME B ( biomass , here) shouldequal BACKGROUND CONC B in the
CHEMISTRY CARD . Finally, the initial porosity of the system is set by assigning thevolume fraction of
NAME BIOMINERAL ( chubbite , here). In general, bold text is required. However, other options maybe specified, if desired. 'temperature points' 8 0. 25. 60. 100. 150. 200. 250. 300. 'H2O' 3.0 0.0 18.0153 'Cr(VI)' 0. 0. 0. 'molasses' 0. 0. 0. 'ethanol' 0. 0. 0. 'biocide' 0. 0. 0. 'null' 0 0 0 'null' 1 0. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 'null' 0. 1 1. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0. 'chubbite' 1.0 0 0. 0. 0. 0. 0. 0. 0. 0. 1.0 'null' 0. 1 0. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0. 'null' 1 0. '0' 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. Figure 3: Minimal
CHROTRAN chemistry database. The text shown here should not be removed, however additionalspecies may be added, if desired. See
PFLOTRAN user manual for details on the database format.11 o ethanol co-injection Ethanol co-injection N o d i r e c t a b i o t i c r e du c t i o n D i r e c t a b i o t i c r e du c t i o n Figure 4: Maps of Cr(VI) concentrations [ppb] in the aquifer 470 d after injection ceased in each of the four scenariosdiscussed in Section 4.1. Injection well location is denoted by a black X .ethanol is attributable to the fact that ethanol co-injection prevented consumption of molasses by the biomass duringthe injection phase, and so molasses persists over a larger area. CHROTRAN has the capability to model hydraulic conductivity reduction due to bio-fouling and the use of biocideas a remediation strategy. To illustrate model capabilities, we perform a simulation of constant-head injection intoa homogeneous aquifer in which the injection fluid amended initially with acetate ( D = − M) for the first 400 d,which is subsequently replaced with dithionite ( X = . CHROTRAN input file is the same as in the study outlined in Section 4.1 (this is to say, as shown in Figures 1and 2), but with different
CHROTRAN parameter values, as shown in Table 2. We here make the reasonable (Ritmann,2004, p. 361) assumption that biomass has the same density as water (recall that we everywhere use the interpretationthat 1 mol of biomass is defined as 1 g of biomass).The simulation is performed on a 50 m square homogeneous hydraulic conductivity field, with constant hydraulicconductivity K = . × − m s − . Each simulation takes place over a span of 500 days, and begins with ε = − CHROTRAN parameter values used in the bio-fouling example in Section 4.2.Symbol Value Units α B min − mol L − Γ B . × − L mol − s − Γ CD − s − Γ X . × − L mol − s − λ B − s − λ B − s − λ C − s − λ D − λ D i − λ D m − s − K B × mol L − K C − M K D − M K I ρ B mol L − S C . × − - S D − - S D . × − -initial concentrations of all species, except B ( x , ) = B min = − mol L − and C = . × − mol L − . No flowboundaries are imposed at the north and south of the domain ( y = y =
50 m), and constant head boundariesare imposed at the west and east of the domain (head 0.28 m at x = x =
50 m). A single injectionwell exists at ( x , y ) = (
25 m ,
25 m ) , and constant head of 0.28 m is imposed at its location.A sequence of quiver plots representing the velocity field at nine points in time, superimposed on the intensity ofbiomass concentration are shown in Figure 5. During the first 400 d of the simulation, biomass concentration growsin the vicinity of the well, until hydraulic conductivity drops to zero at the well until no influx occurs there; onlyambient flow is apparent, flowing around the impermeable barrier near the well. At this point, the biomass has becomeuseless for bioremediation, as contaminated aquifer water no longer travels through it. However, at 400 d, dithionite isintroduced into the injection fluid and effectively eliminates biomass in the vicinity of the well. The region containingdithionite is relatively sterile and grows outwards until the biomass concentration approaches background, and theinitial flow regime is recovered at 416 d. Because initial and final conditions are the same, this cycle may be performedindefinitely. For modeling in situ remediation of aqueous groundwater contaminants by injection of aqueous amendment, we rec-ognized the importance of mathematical formulations and numerical codes that can represent three-dimensional fluidflow and multi-species contaminant transport in heterogeneous aquifers with arbitrary injection regimes. For the par-ticularly important case of heavy metal remediation, a number of contaminant-remediation processes (pathways) aresusceptible to a unified modeling framework: bio-reduction, bio-precipitation, and direct reduction by the chemicalamendment. There have previously existed no general tools appropriate for modeling such interventions. With thisbackground in mind, we developed a mathematical model that describes the reactive transport dynamics of an amend-ment (containing any combination of electron donor, non-lethal bio-inhibitor, and biocide) with biomass and aqueousheavy metal contaminant. We also implemented the mathematical model in a novel computational framework, called
CHROTRAN , that is based on the open-source code
PFLOTRAN . PFLOTRAN ’s modularity and the reaction sandboxcapability allowed us to implement the model easily without making any changes to the flow and transport part of13igure 5: A sequence of snapshots of the cell-center groundwater seepage velocity fields and biomass concentrationdistributions in the example of Section 4.2. Velocity magnitude is indicated by arrow length and direction by arroworientation; the arrow tails are located at the cell center; biomass concentration [g m − ] is indicated by green intensityin each superimposed map. The initial condition snapshot is shown in the upper left corner, with time increasing inthe clockwise direction, until the initial condition is reached again at 416 d. The same scale is used in each snapshot.14 FLOTRAN . CHROTRAN can harness the high-performance computing capabilities of
PFLOTRAN which allows forsimulations of complex models with large number of computational cells and degrees of freedom. We described ourcomputer implementation and explained how to use
CHROTRAN to solve practical problems.We also considered two demonstration studies related to chromium remediation. The presented synthetic problemswere formulated to be consistent with real-world groundwater contamination problems and illustrate the capabilityof
CHROTRAN to aid in the engineering design process. In one of the studies, we discovered that, contrary to muchexisting theory, Cr(VI) reduction was maximized by injecting molasses and suppressing biomass growth to maximizethe direct, abiotic reduction reaction. In the other, we showed the feasibility of pulsed injection of bio-stimulant andbiocide to alleviate bio-fouling in the context of ongoing bioremediation.We observe that because of the abstraction of our model and its parametric flexibility, the
CHROTRAN equations canbe used to model other reactive transport behaviors besides the heavy metal bio-reduction that we have focused upon,including basic advection-dispersion-reaction interaction (between C and D , in the absence of B ). The bio-reductionmodel captures any biodegradation that can be represented using a Monod equation, as long as the contaminant rep-resented by C is non-sorbing, and it does not explicitly require the contaminant to be reduced. This potentially allowsfor modeling the biodegradation of a wide range of organic contaminants, which include but are not limited to hydro-carbons, chlorinated solvents, pesticides, and volatile organic compounds. Acknowledgments
The authors acknowledge the support of the LANL Environmental Programs.
References
Alam, M. (2004).
Bioreduction of hexavalent chromium: flow-through column experiments and reactive transportmodeling . Ph. D. thesis, Washington State University.Andreazza, R., S. Pieniz, L. Wolf, M. K. Lee, F. A. O. Camargo, and B. C. Okeke (2010). Characterization of copperbioreduction and biosorption by a highly copper resistant bacterium isolated from copper-contaminated vineyardsoil.
Science of the Total Environment 408 (7), 1501–1507.Appelo, C. A. J. and D. Postma (2004).
Geochemistry, groundwater and pollution . CRC press.Balay, S., S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. Gropp,D. Karpeyev, D. Kaushik, M. Knepley, L. Curfman McInnes, K. Rupp, B. Smith, S. Zampini, H. Zhang, andH. Zhang (2016). PETSc users manual. Technical Report ANL-95/11 Rev 3.7, Argonne National Laboratory.Chen, Z.-F., Y.-S. Zhao, J.-W. Zhang, and J. Bai (2015). Mechanism and kinetics of hexavalent chromium chemicalreduction with sugarcane molasses.
Water, Air, & Soil Pollution 226 (11), 1–9.Chiang, C. Y., C. N. Dawson, and M. F. Wheeler (1991). Modeling of in-situ biorestoration of organic compounds ingroundwater.
Transport in Porous Media 6 (5-6), 667–702.Duruibe, J., M. Ogwuegbu, J. Egwurugwu, et al. (2007). Heavy metal pollution and human biotoxic effects.
Interna-tional Journal of Physical Sciences 2 (5), 112–118.Hammond, G. E. (2015). PFLOTRAN: Recent developments facilitating massively-parallel reactive biogeochemicaltransport. In . American Geophysical Union Fall Meeting.Hansen, S. K., H. Boukhalfa, S. Karra, D. Wang, and V. V. Vesselinov (2016). Chromium (VI) reduction in acetate-and molasses-amended natural media: experimental results and model development. Technical Report LA-UR-16-27473, Los Alamos National Laboratory.Hashim, M. A., S. Mukhopadhyay, J. N. Sahu, and B. Sengupta (2011). Remediation technologies for heavy metalcontaminated groundwater.
Journal of Environmental Management 92 (10), 2355–2388.15stok, J. D., J. M. Senko, L. R. Krumholz, D. Watson, M. A. Bogle, A. Peacock, Y. J. Chang, and D. C. White (2004).In Situ Bioreduction of Technetium and Uranium in a Nitrate-Contaminated Aquifer.
Environmental Science andTechnology 38 (2), 468–475.Li, L., N. Gawande, M. B. Kowalsky, C. I. Steefel, and S. S. Hubbard (2011). Physicochemical heterogeneity controlson uranium bioreduction rates at the field scale.
Environmental Science & Technology 45 (23), 9959–66.Li, L., C. I. Steefel, M. B. Kowalsky, A. Englert, and S. S. Hubbard (2010). Effects of physical and geochemicalheterogeneities on mineral transformation and biomass accumulation during biostimulation experiments at Rifle,Colorado.
Journal of Contaminant Hydrology 112 (1-4), 45–63.Lichtner, P. C., G. E. Hammond, C. Lu, S. Karra, G. Bisht, B. Andre, R. T. Mills, and J. Kumar (2015). PFLOTRANuser manual: A massively parallel reactive flow and transport model for describing surface and subsurface processes.Technical Report LA-UR-15-20403, Los Alamos National Laboratory.Lovley, D. R. (1993). Dissimilatory metal reduction.
Annual Reviews in Microbiology 47 (1), 263–290.Lovley, D. R. (1995). Bioremediation of organic and metal contaminants with dissimilatory metal reduction.
Journalof industrial microbiology 14 (2), 85–93.Malik, A. (2004). Metal bioremediation through growing cells.
Environment international 30 (2), 261–278.Molins, S., J. Greskowiak, C. Wanner, and K. U. Mayer (2015). A benchmark for microbially mediated chromiumreduction under denitrifying conditions in a biostimulation column experiment.
Computational Geosciences 19 (3),479–496.Monod, J. (1949). The growth of bacterial cultures.
Annual Reviews in Microbiology 3 (1), 371–394.Okeke, B. C. (2008). Bioremoval of hexavalent chromium from water by a salt tolerant bacterium, Exiguobacteriumsp. GS1.
Journal of Industrial Microbiology and Biotechnology 35 (12), 1571–1579.Priester, J. H., S. G. Olson, S. M. Webb, M. P. Neu, L. E. Hersman, and P. A. Holden (2006). Enhanced exopolymerproduction and chromium stabilization in Pseudomonas putida unsaturated biofilms.
Applied and EnvironmentalMicrobiology 72 (3), 1988–1996.Radhika, V., S. Subramanian, and K. A. Natarajan (2006). Bioremediation of zinc using Desulfotomaculum nigrifi-cans: Bioprecipitation and characterization studies.
Water Research 40 (19), 3628–3636.Ritmann, B. E. (2004). Biofilms in the water industry. In M. Ghannoum and G. O’Toole (Eds.),
Microbial Biofilms .Washington, DC: ASM Press.Shashidhar, T., S. M. Bhallamudi, and L. Philip (2007). Development and validation of a model of bio-barriers forremediation of cr(vi) contaminated aquifers using laboratory column experiments.
Journal of Hazardous Materi-als 145 (3), 437–452.Shashidhar, T., L. Philip, and S. M. Bhallamudi (2006). Bench-scale column experiments to study the containment ofcr(vi) in confined aquifers by bio-transformation.
Journal of Hazardous Materials 131 (13), 200 – 209.Somasundaram, V., L. Philip, and S. M. Bhallamudi (2009). Experimental and mathematical modeling studies onCr(VI) reduction by CRB, SRB and IRB, individually and in combination.
Journal of Hazardous Materials 172 (2-3), 606–617.Steefel, C. I., C. A. J. Appelo, B. Arora, D. Jacques, T. Kalbacher, O. Kolditz, V. Lagneau, P. C. Lichtner, K. U. Mayer,J. C. L. Meeussen, S. Molins, D. Moulton, H. Shao, J. ˇSim˚unek, N. Spycher, S. B. Yabusaki, and G. T. Yeh (2015).Reactive transport codes for subsurface environmental simulation.
Computational Geosciences 19 (3), 445–478.Suk, H., K.-K. Lee, and C. H. Lee (2000). Biologically reactive multispecies transport in sanitary landfill.
Journal ofEnvironmental Engineering 126 (May), 419–427.Tchounwou, P. B., C. G. Yedjou, A. K. Patlolla, and D. J. Sutton (2012). Heavy metal toxicity and the environment.In
Molecular, clinical and environmental toxicology , pp. 133–164. Springer.16ravis, B. J. (1993). Numerical simulation of in situ bioremediation. In
Transactions on Ecology and the Environment ,Volume 17, pp. 91–98.Travis, B. J. and N. D. Rosenberg (1998). Modeling in situ bioremediation of TCE at Savannah River: Effects ofproduct toxicity and microbial interactions on TCE degradation.
Environmental Science and Technology 31 (11),3093–3102.Van Roy, S., K. Vanbroekhoven, W. Dejonghe, and L. Diels (2006). Immobilization of heavy metals in the saturatedzone by sorption and in situ bioprecipitation processes.
Hydrometallurgy 83 (1), 195–203.Violante, A., V. Cozzolino, L. Perelomov, A. G. Caporale, and M. Pigna (2010). Mobility and bioavailability of heavymetals and metalloids in soil environments.
Journal of soil science and plant nutrition 10 (3), 268–292.Wang, Y.-T. and C. Xiao (1995). Factors affecting hexavalent chromium reduction in pure cultures of bacteria.
WaterResearch 29 (1), 2467–2474.Wheeler, M. F., K. R. Roberson, and A. Chilakapati (1992). Three-dimensional bioremediation modeling in heteroge-neous porous media. In
University of Colorado 9th International Conference on Computational Methods in WaterResources , Denver, CO, pp. 1–19.Wu, G., H. Kang, X. Zhang, H. Shao, L. Chu, and C. Ruan (2010). A critical review on the bio-removal of hazardousheavy metals from contaminated soils: Issues, progress, eco-environmental concerns and opportunities.
Journal ofHazardous Materials 174 (1-3), 1–8.Yabusaki, S. B., Y. Fang, K. H. Williams, C. J. Murray, A. L. Ward, R. D. Dayvault, S. R. Waichler, D. R. Newcomer,F. A. Spane, and P. E. Long (2011). Variably saturated flow and multicomponent biogeochemical reactive transportmodeling of a uranium bioremediation field experiment.
Journal of Contaminant Hydrology 126 (34), 271 – 290.Zhan, G., D. Li, and L. Zhang (2012). Aerobic bioreduction of nickel(II) to elemental nickel with concomitantbiomineralization.