CHROTRAN: A mathematical and computational model for in situ heavy metal remediation in heterogeneous aquifers

TL;DR

CHROTRAN is a comprehensive 3D reactive transport model for in situ heavy metal groundwater remediation, incorporating complex biogeochemical processes and heterogeneous flow, demonstrated through example simulations.

Contribution

The paper introduces CHROTRAN, a novel computational model that integrates multiple biogeochemical processes for in situ heavy metal remediation in heterogeneous aquifers.

Findings

Model handles multiple chemical species and injection points.

Simulations demonstrate unique remediation strategies.

Software supports complex flow and reactive transport coupling.

Abstract

Groundwater contamination by heavy metals is a critical environmental problem for which in situ remediation is frequently the only viable treatment option. For such interventions, a three-dimensional reactive transport model of relevant biogeochemical processes is invaluable. To this end, we developed a model, CHROTRAN, for in situ treatment, 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 as donor-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 many chemical 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. One simulation suggests an unorthodox approach to remediation of Cr(VI) contamination.

Figures

Equations

$$\frac{\partial C}{\partial t}\propto \frac{\partial D}{\partial t}\propto B\frac{C}{K_C+C}\frac{D}{K_D+D},$$

$$\frac{\partial C}{\partial t}\propto \frac{\partial D}{\partial t}\propto B\frac{K_{C^{\prime }}}{K_{C^{\prime }}+C}\frac{D}{K_D+D}.$$

$$\frac{\partial C}{\partial t}\propto B\frac{C}{K_C+C}.$$

$$\frac{d}{dt}\left({\rho }_w\theta \right)+\nabla \cdot \mathbf{q}=q_M(\mathbf{x},t),$$

$$\mathbf{q}=-\nabla \left({\rho }_{w}K(\mathbf{x},\mathbf{t})h(\mathbf{x},t)\right),$$

$$K(\mathbf{x},t)=K(\mathbf{x},0)\frac{\theta (\mathbf{x},t)}{{\theta }_0},$$

$$\theta (\mathbf{x},t)={\theta }_0-\frac{B(\mathbf{x},t)}{{\rho }_B},$$

$$\mathcal{T}\{c\}\equiv -\mathbf{q}\cdot \nabla c+\nabla \cdot (\theta \mathbf{D}(\mathbf{q})\nabla c),$$

$${\mu }_B={\lambda }_{B_1}B\stackrel{e^{-}\mathrm{donor}}{\overbrace{\frac{D}{K_D+D}}}\underset{\mathrm{crowding}}{\underbrace{{\left(\frac{K_B}{K_B+B}\right)}^{\alpha }}}\stackrel{\mathrm{inhibition}}{\overbrace{\frac{K_I}{K_I+I}}}\left[\frac{\mathrm{mol}}{{\mathrm{L}}_{\mathrm{b}}\mathrm{s}}\right].$$

$${\mu }_{CD}={\Gamma }_{CD}CD\left[\frac{\mathrm{mol}}{{\mathrm{L}}_{\mathrm{b}}\mathrm{s}}\right].$$

$$\frac{\partial \theta C}{\partial t}$$

$$\frac{\partial \theta D_m}{\partial t}$$

$$\frac{\partial \theta I}{\partial t}$$

$$\frac{\partial \theta X}{\partial t}$$

$$\frac{\partial B}{\partial t}$$

$$\frac{\partial D_i}{\partial t}$$

Full text

Keywords: in situ remediation, abiotic reduction, bioremediation model, subsurface flow, reactive transport

1 Introduction

Heavy metals, including chromium, arsenic, copper, nickel, selenium, technetium, uranium, and zinc, are widespread and 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, the chemical reduction of certain metals is known to reduce their mobility (Violante et al., 2010). This has inspired efforts to 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 been noted as a viable bioremediation technique by recent critical reviews (Hashim et al., 2011; Wu et al., 2010). Bioprecipitation, a process by which microbiological exudates react with metals to produce an insoluble compound, has been widely observed (Malik, 2004; Van Roy et al., 2006; Radhika et al., 2006) and has been noted by Wu et al. (2010) as a remediation method. Bio-stimulants have also recently been shown to effectively reduce chromium through abiotic redox 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 a multi-dimensional/multi-component numerical model of groundwater flow, contaminant transport, and biogeochemical processes to evaluate different remediation strategies under varying field conditions. The model should be capable of 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 is not vast. Early work focused on in situ bioremediation of toxic organic compounds through oxidation. A thorough mathematical and 2D numerical study representative of this approach is due to Chiang et al. (1991), who presented a three-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, which was 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 $\mathrm{CH_{4}}$. Travis (1993) presented a more complicated, unsaturated three-dimensional model, which introduced Monod dependence on nutrients, and the potential for two electron donors, with one inhibiting the other. This approach was further elaborated 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). Another complex oxidation model was developed by Suk et al. (2000), which explicitly modeled both mobile and immobile biomass, 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 studies have 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 externally applied electron donor, and the implied equations have a similar form to those devised by Chiang et al. (1991): linear in biomass, Monod in contaminant, and Monod in electron donor. For clarity, this is expressed symbolically as:

[TABLE]

where $C$, is the U(VI) concentration, $B$ is the biomass concentration, $D$ is the electron donor concentration, $K_{C}$ is the metal reduction Monod constant, and $K_{D}$ is the electron donor Monod constant. $K_{C}$ and $K_{D}$ respectively represent the concentration of $C$ and $D$ at which the reaction rate is halved. Recently, Molins et al. (2015) have published a numerical study of a column experiment with multiple species, all of whose dynamics are of the above form, but including an extra chemical inhibition factor. The models of Li et al. (2010, 2011) were implemented at field-scale in CrunchFlow (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 by Alam (2004), and by Shashidhar et al. (2007). Shashidhar et al. (2007) described the Cr(VI) degradation reaction slightly differently from Li et al. (2011):

[TABLE]

$K_{C^{\prime}}$ 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 but rather as an inhibitor. Interestingly, since the RHS of (2) is a proxy for the biomass growth reaction, $C$ consumption is modeled as proportional to biomass growth, but the biomass growth rate is modeled as independent of $C$. Biomass dynamics were governed by a growth term proportional to donor consumption and a first-order decay term, accounting for 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, and also included two enzymes (both created due to biomass growth, but one conserved, and one irreversibly consumed during bio-reduction). Neglecting the irreversibly-consumed enzyme and the mobile-immobile behavior, this model shares its electron donor and biomass dynamics with the model of Shashidhar et al. (2007). It differs significantly from other models that we are aware of by treating the Cr(VI) degradation reaction in this model as an incidental enzymatic process, and is governed by the following Monod equation:

[TABLE]

There is strong experimental support for this approach (e.g., Okeke, 2008), and this is arguably more defensible in a real complex geochemical system in which there are multiple competing donors and receptors, and given that there is 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 thus appears that the primary example of a bio-reduction model applicable to modeling a real-world remediation scheme is the CrunchFlow model of uranium treatment at the Rifle site, which was discussed above. We set out to develop a new model, dubbed chrotran, which is optimized for modeling bioremediation of Cr(VI), but of sufficient generality that it might be used for bioremediation of other metals, or for abiotic reduction, with ease. The key features of the model 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 donor and Cr(VI) is the contaminant, rather than the bio-mediated reaction previously posited. It is thus crucial to include 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) for modeling Cr(VI) degradation with (3), we implemented this general formulation as opposed to one which ties all 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 schemes is 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 pass through it. The model thus contains feedback between local biomass concentration and flow parameters such as porosity and hydraulic conductivity.

Biomass crowding

Similarly, if biomass becomes overly dense, this causes cell stress, which reduces the rate of further 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 well before they are consumed, additional chemicals may be injected to reduce biomass concentrations, and their reactive transport behavior is incorporated.

Multiple donor consumption pathways

The best model of electron donor consumption by biomass may be proportional 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 transport codes (Steefel et al., 2015). To accomplish our goal, we turned to pflotran (Lichtner et al., 2015), which is open source, has a modular structure, and a “reaction sandbox” interface (Hammond, 2015) that allows derivative versions with custom reaction behavior to be developed and compiled. In this fashion, no changes to the flow and transport part of pflotran are needed. We developed chrotran based on the existing pflotran code framework, taking advantage of the reaction sandbox interface to implement complex model features not included in basic microbial packages 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 the model. In Section 3, we discuss how to install and use the software. In Section 4, we present two numerical studies which illustrate chrotran and also suggest an interesting conclusion regarding Cr(VI) remediation. In Section 5, we briefly summarize what has been presented.

2 Model description

We consider flow and transport at aquifer scale. Conceptually, the aquifer is modeled as saturated, with incompressible water moving in accordance with Darcy’s law (chrotran can also simulate partially-saturated vadose-zone flow and transport). 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 the model is intended to be used for remedial design, every effort was made to simplify the formulation to use the smallest number of explanatory variables and parameters, and to keep the equations at a high level of abstraction, so they are not 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 with their own symbols:

Biomass, $B\>\mathrm{[mol\>L_{b}^{-1}]}$,

representing the concentration of all microbes and their associated extracellular material. The quantification 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 presentation of the model, so all governing equations have the same units. A mole of biomass should be understood as an equivalent mass: any quantity can be used, as long as one uses a consistent definition throughout the model. In the examples in this paper, we use the definition 1 mol $\equiv$ 1 g of biomass.

Aqueous contaminant, $C\>\mathrm{[mol\>L^{-1}]}$,

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 be

1. 1.

immobile, represented by $D_{i}\>\mathrm{[mol\>L_{b}^{-1}]}$, or 2. 2.

mobile, represented by $D_{m}\>\mathrm{[mol\>L^{-1}]}$,

with exchange of mass between the two states.

Nonlethal biomass-growth inhibitor, $I\>\mathrm{[mol\>L^{-1}]}$,

such as ethanol, which is modeled as a conservative species but acts to slow microbial growth.

Biocide, $X\>\mathrm{[mol\>L^{-1}]}$,

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 the formula $D=\frac{D_{i}}{\theta(\bm{x},t)}+D_{m}\>\mathrm{[mol\>L^{-1}]}$, where $\theta(\bm{x},t)$ [-] is the current porosity at $\bm{x}$. For simplicity, we assume that both the mobile and immobile donor participate equally in all reactions.

2.1 Flow and transport

2.1.1 Groundwater flow equations

Flow may be modeled using the balance of water mass given by

[TABLE]

with the water mass fluxes related to head via Darcy’s law:

[TABLE]

where $K(\bm{x})\ \mathrm{[m/s]}$ is the local hydraulic conductivity and $h(\bm{x},t)\ \mathrm{[m]}$ is the local hydraulic head, $\theta$ is the porosity and $q_{M}(\bm{x})\ \mathrm{[kg\ m^{-3}\ s^{-1}]}$ is the local mass injection rate into the system and $\rho_{w}\ \mathrm{[kg\ m^{-3}]}$ is the density of water. We note that, since chrotran is built on top of pflotran, it inherits all of pflotran’s groundwater flow modeling capabilities. This includes the ability to consider unsaturated and otherwise multiphase flow conditions, which are out-of-scope for the present discussion. Please see the pflotran user manual (Lichtner et al., 2015) for details on its complete capabilities.

The hydraulic conductivity is continually updated in accord with the relation

[TABLE]

where $\theta_{0}$ [-] is the spatially-uniform initial porosity, and $\theta(\bm{x},t)$ is calculated according to

[TABLE]

where $\rho_{B}\ \mathrm{[mol\ L^{-1}]}$ is the intrinsic biomass density. (Note that, using our proposed definition of 1 mol of biomass as 1 g of biomass, 1 $\mathrm{mol\ L^{-1}}=$ 1 $\mathrm{kg\ m^{-3}}$.)

2.1.2 Advective-dispersive transport operator

We define $\mathscr{T}\{\cdot\}$ to be an advective-dispersive transport operator, which characterizes the hydrodynamic effects on solute transport. For $c$, the concentration of an arbitrary mobile species,

[TABLE]

where $\bm{D}$ is a dispersion tensor that depends on the longitudinal and transverse dispersivities, molecular diffusion as well as the Darcy flux. For the work in this paper, we will only consider molecular diffusion and thereby we set $\bm{D}=D_{m}\bm{I}.$ Note that, while this is not shown explicitly for compactness, all symbols in this equation are functions of $\bm{x}$ and $t$.

2.2 Biogeochemical reactions

We define one governing equation for each species, mobile or immobile, as well as two equations defining reaction rate expressions for algebraic convenience. The equations involve numerous parameters, whose symbols, units, and long-form name in the chrotran input file are summarized in Table 1. The parameter symbols follow a scheme in which the first letter encodes the physical interpretation of the parameter and the subscript specifies the governing equation in which they participate. A symbol beginning with $\Gamma$ is a second-order mass action rate constant, with units $[\mathrm{L\ mol^{-1}\ s}]$. A symbol beginning with $K$ is Monod or inhibition constant with units of concentration, $\mathrm{mol\ L_{b}^{-1}}$ or $\mathrm{mol\ L^{-1}}$, and represents the concentration at which a process rate becomes 50% of its maximum rate, all other parameters being equal. A symbol beginning with $\lambda$ has units of $\mathrm{s^{-1}}$ and is interpreted as a pure first-order reaction rate constant. A symbol beginning with $S$ is dimensionless, and represents a stoichiometric relationship between a reaction rate and the consumption rate of a certain species. Before presenting the equations, it is useful to review all of 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$, and Monod in $C$, which is a common assumption for bio-mediated processes. Note that we are not assuming that reduction of $C$ is directly tied to any particular cell metabolic process. This form is sufficiently general that it can capture other bio-remediation processes besides bio-reduction of heavy metals.

Biocide reaction

This is an inter-phase bimolecular reaction between the biocide, $X$, and the biomass, $B$. It is modeled with a classical second-order mass action rate law, with the added condition that $B$ cannot fall below a specified minimum concentration $B_{\mathrm{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 inhibition effects are assumed: a biomass crowding term, tunable with exponent $\alpha$, attenuates growth rate as the biomass concentration rises. The nonlethal inhibitor concentration, $I$, also reduces the reaction rate as its concentration increases.

Mobile-immobile mass transfer (MIMT)

This is a process with first-order kinetics, which models sorptive retardation of the electron donor.

Natural decay

This is an empirical process reflecting the idea that, if left unstimulated, both the amount of living cells and the amount of extracellular material in the aquifer will ultimately return to their natural background level (i.e. $B_{\mathrm{min}}$). This is modeled as a first-order process.

Respiration

This represents consumption of the electron donor for purposes of life maintenance, unrelated to biomass growth. 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 and terms in the governing equations presented below.

2.3 Reactive transport equations

2.3.1 Definitions of convenience reaction variables

The biomass growth reaction is linear in biomass concentration, has a Monod dependency on electron donor, a tunable inhibition factor due to biomass crowding, and a classic inhibition factor describing the impact of the nonlethal growth inhibitor (as indicated by comment braces):

[TABLE]

The direct, abiotic reduction reaction is represented by a classic, second-order mass action law:

[TABLE]

2.3.2 Partial differential equations for mobile chemical components

The mobile components are all governed by the advection-dispersion operator, $\mathscr{T}$, defined previously, and also affected by extra terms implementing the chemical processes outlined earlier (as indicated by comment braces):

[TABLE]

2.3.3 Partial differential equations for immobile chemical components

The immobile component concentrations are affected only by the reactive processes outlined above:

[TABLE]

3 Software implementation

3.1 Software availability and installation

chrotran currently is undergoing an open source licensing process, and will be available as soon as possible as freely downloadable Fortran 2003 source code. Until this process is complete, the software may be made available on a restricted basis by arrangement with the authors.

chrotran must be compiled using the GFortran complier (freely available as part of the GNU Compiler Collection). It is based on the open-source pflotran code base, and the installation procedure is the essentially the same as 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 installation instructions333Available at http://documentation.pflotran.org/user_guide/how_to/installation/installation.html. are applicable, except that chrotran, rather than pflotran, should be cloned from its repository444Not officially released yet; anonymous access can be facilitated through editor for review purposes. once all the dependencies have installed. To build chrotran itself, navigate to <path of cloned repository>/src/pflotran and type make chrotran. (The chrotran executable will be called chrotran.)

3.2 Specifying and running a simulation

A chrotran input file is of the same format as a pflotran input file. Information on how to set up such a file is available 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 Figure 2, again with bold text being mandatory and standard-weight text being user-alterable. Comments in the input file are preceded 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 shown in Figure 3. The one exception to bold text being mandatory is that species names can be changed at will, as long as there is consistency between the CHEMISTRY card and the chemistry database. For instance, one could change all instances 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 kinetics listed 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 specified as inactive with respect precipitation/dissolution by setting its kinetic rate constant (RATE_CONSTANT) to zero. The mineral is included as a surrogate for biomass and porous media volume in chrotran and is updated according to Equation 7 to track $1-\theta(\bm{x},t)$. The initial volume fraction of chubbite thus defines the initial porosity. The format of 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 <path of cloned repository>/src/pflotran, and type chrotran -pflotranin <path to input file>/test.in. The output of the simulation will be saved in the same directory as the input file. Depending on the options specified in the input file, chrotran can save flow field velocities, concentrations of all species, permeabilities, and porosities at any specified times in an .h5 format file. This file format can be visualized natively using freely-available standalone tools such as VisIt and ParaView, and are also accessible from Python scripts by means of the h5py library and from Julia scripts by means of the HDF5 package.

3.3 Quality assurance

The pflotran software from which chrotran derives its numerical flow and reactive transport solvers has gone through extensive quality assurance testing, has been benchmarked against other reactive transport solvers, and is used inside and outside the U.S. Department of Energy for mission-critical analytical work. The new bio-reactive transport model that is constitutive of chrotran is not available in any other software, so direct benchmarking is not possible. However, extensive quality testing has been performed by the developers. We have validated through batch and multidimensional simulations that chrotran does satisfy the governing equations we present for chemistry and permeability, and also that it gives plausible, physically consistent results for a wide range of scenarios.

4 chrotran remediation case studies

To demonstrate the capabilities of our software, we present two example studies, which together illustrate the interactions of all the types of chemical species it permits to be modeled, along with its treatment of bio-clogging. The input files for these two examples can be found in the chrotran_examples directory of the chrotran repository.

4.1 Remediation of Cr(VI) by molasses and ethanol co-injection

This study concerns the co-injection of molasses (electron donor, $D$) and ethanol (nonlethal bio-inhibitor, $I$) into a single well drilled in a heterogeneous aquifer with an appreciable background Cr(VI) concentration. The competition between direct abiotic reduction of Cr(VI) by molasses and bio-reduction of Cr(VI), which exists since both reduction pathways consume the electron donor, along with the impact of suppressing the biomass growth is explored. The basic parameters used are those shown in Figures 1 and 2, with changes as indicated below.

Four related simulations are performed on the same 100 $\times$ 100 m two-dimensional heterogeneous hydraulic conductivity field, with geometric mean conductivity $K_{g}=10^{-4}\ \mathrm{m\ s^{-1}}$, a multi-Gaussian correlation structure with exponential semivariogram with correlation length of 4 m and $\sigma^{2}_{\ln K}=2$. Each simulation takes place over a span of 500 days and begins with $\epsilon=10^{-20}$ initial concentrations of all species, except $B(\bm{x},0)=B_{min}=10^{-10}\ \mathrm{mol\ L_{b}^{-1}}$ and $C=1.923\times 10^{-5}\ \mathrm{mol\ L^{-1}}$ (1000 ppb Cr(VI)). In all cases, there are no flow boundaries at the north and south of the domain ($y=0$ m and $y=100$ m), and constant head boundaries are imposed at the west and east of the domain ($x=0$ m and $x=100$ m) such that there is a drop of head of 0.28 m between these faces. A single injection well exists at $(x,y)=(25\;\mathrm{m},50\;\mathrm{m})$. For the first 10 days of the simulation, there is no injection into the well. From 10 d to 30 d, injection is performed at the well with constant volumetric flow rate 272.55 $\mathrm{m^{3}d^{-1}}$ with species concentrations discussed below. From 30 d to 500 d, there is again no injection at the well. A very large (arbitrary) $\rho_{B}$ is assumed, so as 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: $\Gamma_{CD}=1\ \mathrm{L\ mol^{-1}\ s^{-1}}$ and $\Gamma_{CD}=0\ \mathrm{L\ mol^{-1}\ s^{-1}}$, as are two different ethanol concentrations in the injection fluid: $I=1\ \mathrm{mol\ L^{-1}}$ and $I=\epsilon\ \mathrm{mol\ L^{-1}}$, in all four possible combinations. The injection fluid chemistry always has Cr(VI) concentration equal to the initial concentration ($C=1.923\times 10^{-5}\ \mathrm{mol\ L^{-1}}$), ensuring that no chromium disappearance is due to dilution, and molasses concentration $D=1\times 10^{-2}\ \mathrm{mol\ L^{-1}}$.

Concentrations of Cr(VI) for each scenario are shown in Figure 4. It is apparent that little persistent reduction due to biomass alone occurs, although ethanol co-injection does increase biomass footprint, which has a noticeable and persistent 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) are relatively low, and it has a retardation factor of around 150 (obtained from Shashidhar et al. (2006)), meaning that it has the potential to form a persistent permeable reactive barrier around the well. The better performance in the presence of ethanol is attributable to the fact that ethanol co-injection prevented consumption of molasses by the biomass during the injection phase, and so molasses persists over a larger area.

4.2 Biomass clogging/unclogging due to acetate/dithionite injection

chrotran has the capability to model hydraulic conductivity reduction due to bio-fouling and the use of biocide as a remediation strategy. To illustrate model capabilities, we perform a simulation of constant-head injection into a homogeneous aquifer in which the injection fluid amended initially with acetate ($D=10^{-2}$ M) for the first 400 d, which is subsequently replaced with dithionite ($X=3.5$ M) for the remainder of the simulation. The basic structure of the chrotran input file is the same as in the study outlined in Section 4.1 (this is to say, as shown in Figures 1 and 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 interpretation that 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 hydraulic conductivity $K=9.8\times 10^{-5}\ \mathrm{m\ s^{-1}}$. Each simulation takes place over a span of 500 days, and begins with $\epsilon=10^{-20}$ initial concentrations of all species, except $B(\bm{x},0)=B_{min}=10^{-10}\ \mathrm{mol\ L_{b}^{-1}}$ and $C=1.923\times 10^{-5}\ \mathrm{mol\ L^{-1}}$. No flow boundaries are imposed at the north and south of the domain ($y=0$ m and $y=50$ m), and constant head boundaries are imposed at the west and east of the domain (head 0.28 m at $x=0$ m and head 0 m at $x=50$ m). A single injection well exists at $(x,y)=(25\;\mathrm{m},25\;\mathrm{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 of biomass concentration are shown in Figure 5. During the first 400 d of the simulation, biomass concentration grows in the vicinity of the well, until hydraulic conductivity drops to zero at the well until no influx occurs there; only ambient flow is apparent, flowing around the impermeable barrier near the well. At this point, the biomass has become useless for bioremediation, as contaminated aquifer water no longer travels through it. However, at 400 d, dithionite is introduced into the injection fluid and effectively eliminates biomass in the vicinity of the well. The region containing dithionite is relatively sterile and grows outwards until the biomass concentration approaches background, and the initial flow regime is recovered at 416 d. Because initial and final conditions are the same, this cycle may be performed indefinitely.

5 Summary and conclusions

For modeling in situ remediation of aqueous groundwater contaminants by injection of aqueous amendment, we recognized the importance of mathematical formulations and numerical codes that can represent three-dimensional fluid flow and multi-species contaminant transport in heterogeneous aquifers with arbitrary injection regimes. For the particularly important case of heavy metal remediation, a number of contaminant-remediation processes (pathways) are susceptible to a unified modeling framework: bio-reduction, bio-precipitation, and direct reduction by the chemical amendment. There have previously existed no general tools appropriate for modeling such interventions. With this background in mind, we developed a mathematical model that describes the reactive transport dynamics of an amendment (containing any combination of electron donor, non-lethal bio-inhibitor, and biocide) with biomass and aqueous heavy 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 sandbox capability allowed us to implement the model easily without making any changes to the flow and transport part of pflotran. chrotran can harness the high-performance computing capabilities of pflotran which allows for simulations of complex models with large number of computational cells and degrees of freedom. We described our computer implementation and explained how to use chrotran to solve practical problems.

We also considered two demonstration studies related to chromium remediation. The presented synthetic problems were formulated to be consistent with real-world groundwater contamination problems and illustrate the capability of chrotran to aid in the engineering design process. In one of the studies, we discovered that, contrary to much existing theory, Cr(VI) reduction was maximized by injecting molasses and suppressing biomass growth to maximize the direct, abiotic reduction reaction. In the other, we showed the feasibility of pulsed injection of bio-stimulant and biocide 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 can be 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-reduction model captures any biodegradation that can be represented using a Monod equation, as long as the contaminant represented by $C$ is non-sorbing, and it does not explicitly require the contaminant to be reduced. This potentially allows for modeling the biodegradation of a wide range of organic contaminants, which include but are not limited to hydrocarbons, chlorinated solvents, pesticides, and volatile organic compounds.

Acknowledgments

The authors acknowledge the support of the LANL Environmental Programs.