Attribution: You must attribute the work in the manner specified by the author or licensor (but no in any way that suggests that they endorse you or your use of the work).Noncommercial ‒ you may not use this work for commercial purpose.No Derivative works ‒ You may not alter, transform, or build upon this work.Sharing ‒ Individual scientists are hereby granted permission, without fees or further requests to GSA, to use a single figure, a single table, and/or a brief paragraph of text in other subsequent works and to make unlimited photocopies of items in this journal for noncommercial use in classrooms to further education and science.

Abstract

Biologically mediated reductive dissolution and precipitation of metals and radionuclides play key roles in their subsurface transport. Physical and chemical properties of natural aquifer systems, such as reactive iron-oxide surface area and hydraulic conductivity, are often highly heterogeneous in complex ways that can exert significant control on transport, natural attenuation, and active remediation processes. Typically, however, few data on the detailed distribution of these properties are available for incorporation into predictive models. In this study, we integrate field-scale geophysical, hydrologic, and geochemical data from a well-characterized site with the results of laboratory batch-reaction studies to formulate two-dimensional numerical models of reactive transport in a heterogeneous granular aquifer. The models incorporate several levels of coupling, including effects of ferrous iron sorption onto (and associated reduction of reactive surface area of) ferric iron surfaces, microbial growth and transport dynamics, and cross-correlation between hydraulic conductivity and initial ferric iron surface area. These models are then used to evaluate the impacts of physical and chemical heterogeneity on transport of trace levels of uranium under natural conditions, as well as the effectiveness of uranium reduction and immobilization upon introduction of a soluble electron donor (a potential biostimulation remedial strategy).

BACKGROUND AND APPROACH

Microbial Reduction of Metals and Radionuclides

Subsurface contaminants of greatest concern at a number of U.S. Department of Energy (DOE) facilities include metals (e.g., lead, chromium, and mercury) and radionuclides (e.g., uranium, strontium, cesium, technetium, and plutonium). These elements are long-lived and relatively mobile in the subsurface environment and pose potential risks to humans and the natural environment (Lawrence Berkeley National Laboratory, 2003). In situ bioreme-diation is currently being studied as a potential technology for remediation of DOE waste sites. In the more common application of in situ bio-remediation to organic chemicals, microbes can degrade contaminants into daughter products while using the contaminant as a carbon source. Biotransformation of metals and radionuclides, however, must rely on a somewhat different process, since these elemental materials cannot be degraded into any less complex form and do not serve as sources of carbon. Most current research focuses on the potential for metals and radionuclides to be used in microbial respiratory processes as electron donors or acceptors. This process results in changes in valence state of the contaminants and often has significant implications for their mobility in the subsurface.

In this paper, we focus on the microbially mediated reduction of uranium from the soluble U(VI) valence state to the relatively insoluble U(IV) valence state. Metal-reducing micro-organisms obtain energy for cell growth and metabolism through oxidation-reduction reactions, which involve the transfer of electrons from one material (the electron donor) to another (the electron acceptor). It has been shown in laboratory and field studies that naturally occurring microorganisms can utilize dissolved uranium as an electron acceptor, thereby changing the valence state of the uranium and causing the precipitation of U(IV) minerals (Lovley et al., 1991; Lovley and Phillips, 1992; Anderson et al., 2003; Istok et al., 2004). Associated decreases in the concentration of uranium in the mobile aqueous phase reduce risks associated with transport in groundwater to environmental receptors (for example, springs or groundwater-fed streams).

Because uranium has a reduction potential similar to that of iron (and because iron occurs at higher concentration in most aquifers), uranium reduction typically occurs as a cotrans-formation in conjunction with iron reduction as the primary form of microbial respiration. Other substances common in groundwater that have higher reduction potentials than iron (such as oxygen and nitrate) must be removed from the aquifer system before iron and uranium reduction can proceed. Many shallow aquifers do not have sufficient natural carbon and nutrients to support microbial activity necessary to create anaerobic conditions, so in these cases, it is necessary to inject a soluble carbon source to stimulate the necessary level of microbial activity.

Effects of Physical and Chemical Heterogeneity

Spatial variability in aquifer properties (heterogeneity) can exert significant control on the effectiveness of in situ bioremediation processes. In the case of microbial uranium reduction, successful implementation requires that all of the following constituents occur together in the subsurface environment: contaminant (secondary electron acceptor), metal-reducing bacteria, electron donor, and primary electron acceptor. In most cases, the co-occurrence of these constituents will cause sufficient microbial growth and activity to reduce the concentration of any strong oxidants (e.g., oxygen and nitrate) to sufficiently low concentrations and render metal reduction thermodynamically favorable.

Each of these necessary constituents is typically heterogeneous in natural systems. In the case considered here, the contaminant (uranium in the oxidized U[VI] form) is strongly sorbed to naturally occurring iron oxides that are heterogeneously distributed in distinct laminations (see Fig. 1). Natural subsurface bacteria commonly attach to solid grain surfaces and may preferentially attach to fine-grained and/or iron-coated solids that are heterogeneously distributed, thereby creating a heterogeneous distribution of biomass. Heterogeneity in biomass can also arise from spatial variations in labile carbon concentration or other controls on microbial growth; these also can impact the type of organisms favored and therefore lead to variations in microbial community structure. The flow pathways of an injected carbon source (electron donor) are controlled by the physical heterogeneity of the aquifer, in particular the variability in hydraulic conductivity or permeability. Iron oxides often serve as the primary electron acceptor in the microbial respiration process, and, as has already been noted, are heterogeneously distributed. Perhaps most importantly, the heterogeneous distributions of each are controlled by related properties and processes, and can therefore be expected to correlate in important ways. For example, if electron acceptors, contaminants, and bacteria preferentially associate with fine-grained sediments, while the electron donor is preferentially delivered to coarse-grained (more permeable) sediments, this will likely have a negative impact on the success of the field-scale biostimulation. In laboratory microcosms, biologically reacting systems are usually well-mixed and uniform. But in the field, physical and chemical heterogeneity are the rule rather than the exception. Currently, the impacts of joint (correlated) heterogeneity in physical and biogeochemical aquifer properties on the effectiveness of in situ bioremediation are poorly understood. This is primarily because there are few sites where the spatial distribution of these properties has been well-characterized, and because numerical simulation of the multiple interacting processes is complex and computationally demanding.

Numerical Experimentation Approach

We present here the methods and results of a set of numerical experiments. The goal of these experiments is to quantify, under hypothetical but realistic model conditions, the impact of various forms of heterogeneity on the effectiveness of field-scale immobilization of uranium by biologically mediated reduction. The numerical experiments are based on hydrologic, geophysical, and biogeochemical data from an analog field research site (not actually contaminated with uranium) for which extensive characterization data are available. This field site in eastern Virginia, known as the South Oyster site, was used for research on subsurface microbial transport sponsored by the U.S. Department of Energy and has been described in several previous publications (e.g., DeFlaun et al., 1997; Balkwill et al., 2001; Chen et al., 2001; Hubbard et al., 2001; Scheibe et al., 2001; Zhang et al., 2001; Scheibe and Chien, 2003; Mailloux et al., 2003).

We developed a two-dimensional numerical model of field-scale groundwater flow and biogeochemically reactive transport; the model incorporates highly resolved heterogeneous property distributions generated using geostatistical methods. We then used this model to simulate a hypothetical uranium contamination and bioremediation scenario.

METHODS

Biogeochemical Reaction Model

The reaction model includes 27 equilibrium speciation reactions and 5 kinetic reactions, which are listed in 01Table 1. Equilibrium speciation reactions include the major groundwater components (Na-acetate, U[VI], Ca2+, CO32−), as well as the ferrous iron generated during Fe(III) reduction. The speciation system includes surface complexation of U(VI) and Fe(II) according to a generalized composite nonelectrostatic model as described in Davis et al. (1998). Inclusion of Fe(II) surface complexation is critical to modeling microbial Fe(III) oxide reduction because of the inhibitory effect of surface-bound Fe(II) on enzymatic electron transfer to Fe(III) oxides (Roden and Urrutia, 2002). In addition, knowledge of the aqueous-phase versus solid-phase partitioning of biogenic Fe(II) is required in order to estimate the impact of aqueous-phase thermodynamic limitations on Fe(III) oxide reduction rates.

The first kinetic reaction describes the production of Fe(II) via microbial reduction of Fe(III) oxide according to dual Monod kinetics with respect to biomass and acetate and first-order kinetics with respect to “free” Fe(III) oxide surface sites, which are bioavailable, considering the free-energy limit:  
formula
where RFe(III) is the reaction rate for reaction R26 in 01Table 1; [BM] and [Ac] are the concentrations of iron-reducing bacteria biomass and acetate, respectively; Km,BM and Km,Ac,Fe(III) are the Monod half-saturation constants for biomass and acetate, respectively; and ([FesOH] + [FewOH] – [Fe(II)s]) represents the bioavailable Fe(III) oxide surface sites (the total number of strong [s] and weak [w] Fe[III] oxide surface sites minus the total sorbed Fe[II] concentration). Biomass concentration is measured in grams dry weight biomass per cubic decimeter of bulk sediment. This formulation implies that binding sites on the Fe(III) oxide surface occupied by sorbed Fe(II) are not available for microbial reduction, a phenomenon consistent with experimental studies of the geochemical and microbiological controls on this process (Roden and Zachara, 1996; Urrutia et al., 1998; Roden and Urrutia, 1999). Vmax,Fe(III) is the maximum surface-area–normalized rate of enzymatic reduction at saturating cell density; and f (▵Grxn) is a term accounting for the influence of aqueous-phase chemical conditions on the thermodynamic favorability of Fe(III) oxide reduction, defined by Liu et al. (2001) as:  
formula
 
formula
where R is the gas constant.
The total amounts of strong and weak binding sites on iron-oxide surfaces (formed by reactions R31 and R32 in 01Table 1) are computed according to:  
formula
and  
formula
where MW equals the molecular weight of oxide mineral (g mol−1), SA equals the specific surface area of oxide mineral (m2g−1), SSD equals the surface site density (mol m−2), fs is the ratio of strong to total surface sites, φ equals the porosity, and ρ equals the aquifer material bulk density.

The second kinetic reaction is the bacterial reduction of U(VI) according to dual Monod kinetics with respect to aqueous U(VI) and acetate, and first-order kinetics with respect to biomass:

where RU(VI) is the reaction rate for reaction R27 in 01Table 1; [BM] and [Ac] are the concentrations of iron-reducing bacteria biomass and acetate, respectively; Km,U(VI) and Km,Ac,U(VI) are the Monod half-saturation constants for U(VI) and acetate, respectively; and Vmax,U(VI) is the maximum rate of enzymatic reduction.

The third kinetic reaction describes the attachment to and detachment from solid grain surfaces of iron-reducing bacteria as a first-order process with rate coefficients kf (attachment) and kr (detachment), a formulation commonly employed in subsurface microbial transport models (Murphy and Ginn, 2000):  
formula
and  
formula
where [BMs] is the concentration of attached bacteria (cells/g soil).
The fourth and fifth kinetic reactions describe a generalized first-order cell death/maintenance process applied to both aqueous and attached iron-reducing bacteria biomass.  
formula
and  
formula
Where appropriate, reaction model parameters were assigned values from available literature. Stability constants (log K values) for the aqueous equilibrium speciation reactions listed in 01Table 1 were obtained from Nordstrom et al. (1990), with the exception of those for the Fe2+-CO32−-H2O system, which were obtained from Bruno et al. (1992). The total abundance of oxide surface sites was defined by the specified bulk Fe(III) oxide concentration (mol Fe dm−3); an assumed Fe(III) oxide surface area of 350 m2 g−1 (chosen to be representative of poorly crystalline goethite, a common Fe[III] oxide component of soils and sediments; Cornell and Schwertmann, 1996); and the standard mineral surface site density of 3.84 μmol m−2 recommended by Davis and Kent (1990). The generalized composite model for depicting U(VI) and Fe(II) sorption by Fe(III) oxide surfaces sites included both strong and weak sites (reactions R17–R24 in 01Table 1), whose relative abundance could be altered to fit experimental sorption data. However, since both U(VI) and Fe(II) sorption data could be adequately described using a single type of site, we assumed that equal quantities of strong and weak sites were present and used the same log K values for both types of sites. The half-saturation constant for acetate utilization by Fe(III) oxide reduction was set at 100 μM based on recent studies with Shewanella putrefaciens strain CN32 (Liu et al., 2001).

Other model parameters were empirically estimated based on laboratory batch studies using representative sediments collected from an excavation near the South Oyster site. Stability constants for complexation of Fe2+ by Fe(III) oxide surfaces were derived from sorption isotherms computed with data from batch Fe(III) oxide reduction experiments. The Vmax,Fe(III) and Km,BM parameters for Fe(III) oxide reduction were obtained from model fits to experiments in which initial reduction rates were determined as a function of Fe(III)-reducing bacterial cell density. The abundance of free Fe(III) oxide surface sites available for microbial reduction at any given time during the simulation was equal to the total surface density (sum of the total concentrations of strong sites and weak sites) minus the concentration of Fe(II) sorbed to strong and weak oxide surface sites. Thus, rates of Fe(III) oxide reduction were dynamically controlled by the changing abundance of sorbed Fe(II). This same strategy was used successfully in a previous model of synthetic goethite reduction (Roden and Urrutia, 1999), although a different Fe(II) sorption model (Urrutia et al., 1998) was used in the previous study. Currently, use of sorbed Fe(II) abundance as computed by some type of equilibrium speciation expression appears to be the most straightforward and tractable way to model the influence of surface-bound Fe(II) accumulation on Fe(III) oxide reduction activity. The initial sorbed iron-reducing bacteria cell biomass was set equal to 2.0 x 10-5 M, which corresponds approximately to 105 cells per mL bulk sediment assuming the standard formulation of cell material as C5H7O2N and a single cell mass of 5.0 x 10-12 g (Burgos et al., 2002). Subsequent production of Fe(III)-reducing bacterial biomass was computed based on calculated rates of Fe(III) reduction, using yield coefficients for growth of Geobacter metallireducens with Fe(III)-citrate (5–100 mM) as an electron acceptor and excess acetate (20 mM) as a carbon and energy source in batch culture. Cell biomass was determined by measurement of protein content using the bicinchoninic acid (BCA) method (Smith et al., 1985) (Pierce Chemicals) after NaOH digestion. Protein content was converted to dry weight biomass assuming that cells were 50% protein by mass (Madigan and Martinko, 2006). The amount of Fe(III) reduced was determined by diluting 0.25-mL culture samples into 5 mL of 0.5 M HCl, followed by determination of Fe(II) by the ferrozine method (Lovley and Phillips, 1986). Experimental data for growth yield experiments are shown in Figure 2. The results of the yield coefficient determination also agree well with recently published data for growth of G. sulfurreducens with acetate and Fe(III)-citrate in continuous culture (Esteve-Núñez et al., 2005).

The values of Vmax and Km,U(VI) for U(VI) reduction were obtained from U(VI) reduction kinetics studies with washed acetate/fumarate–grown G.sulfurreducens cells under nongrowth conditions in PIPES (Roden and Scheibe, 2005). Km,Ac,U(VI) was assumed to be 200 μM.

Because acetate concentration was present in excess in our study, U(VI) reduction rate did not have a strong dependence on the concentration of acetate.

Equilibrium speciation constants are summarized in 01Table 1, and kinetic reaction model parameters are summarized in 02Table 2.

Aquifer Properties Model

The aquifer properties model is based on data and observations from the South Oyster bacterial transport research site (e.g., Balkwill et al., 2001). Extensive data from the field site were available for use in development of a two-dimensional high-resolution stochastic model of aquifer properties. Data used here include (1) lithologic logs from sediment cores, including visual estimates of sediment type and percent sand, (2) permeability measurements taken on core samples and using a down-hole borehole flow meter, (3) grain-size distributions measured on sediment samples, (4) measurements of extractable iron (Fe[II] and Fe[III]) from core samples, and (5) ground-penetrating radar velocity and attenuation fields based on cross-borehole tomographic measurements. Figure 3 shows vertical logs of several data types collected from one of the three boreholes within the simulated transect. The geophysical data are the most spatially extensive and highly resolved, and (when integrated with borehole data) provide a basis for estimating aquifer properties at locations not directly interrogated by boreholes. Details on the geophysical methods and their application to inference of aquifer properties at the South Oyster site are provided by Chen et al. (2001), Hubbard et al. (2001), and Chen et al. (2004).

Three-dimensional simulations of heterogeneous property distributions have been generated for the South Oyster site by previous investigations (e.g., Scheibe and Chien, 2003). However, the computational demands associated with simulation of reactive transport of uranium and multiple other species over time periods of decades are very high. In addition, the most dense data (from hydrogeophysical methods) are available on two-dimensional transects along series of wells. Accordingly, we chose for this study to limit our simulations to two-dimensional vertical cross sections. We recognize that this decision likely impacts the character of flow and transport in the simulated systems and limits the validity of conclusions drawn. Nevertheless, incorporation of joint physical and chemical heterogeneity into multicomponent reactive transport models at high resolution represents a significant advance, and the results obtained here are likely to motivate investment in computational capabilities necessary to solve similar systems in three dimensions.

Three alternative methods of inferring the spatial distribution of reactive iron (hydr)oxides (Fe[III]) were used in this study.

Method 1: (Chen, 2004). The first method (for conciseness, referred to here as Chen2004) uses geostatistical methods to simulate Fe(III) distributions conditional to direct measurements at boreholes as hard data and geophysical observations as soft data. This method has been described in Chen et al. (2004) and will not be described further here. Chen et al. (2004) did not infer corresponding hydraulic conductivity distributions because the goal of their study was the estimation of geochemical parameters. For the simulated Fe(III) fields based on the methods of Chen et al. (2004), hydraulic conductivity (K) distributions were simulated using the method described in the following based on the corresponding lithology distributions of Chen et al. (2004).

Method 2: (Geophys2). In the current study, we developed an alternative approach (referred to here as “Geophys2”) to jointly simulate both Fe(III) and hydraulic conductivity (K) distributions based on empirical correlations to lithology distributions inferred from geophysical and borehole observations. A two-step process was used to simulate aquifer hydrologic and biogeochemical property distributions. First, geo-statistical simulations of binary sand/mud distributions were generated using the sequential indicator simulation (SISIM) method of GSLIB (Deutsch and Journel, 1998). Then, simulations of hydraulic conductivity and Fe(III) content were created using sequential Gaussian simulation (SGSIM) with locally varying means, where the local mean for a given cell was determined by the corresponding mud/sand simulation. The vertical variogram of the binary indicator of sand and mud (Fig. 4) was estimated from the sand and mud data identified during geologic logging of three cores collected within the simulated transect (taken from the same boreholes in which geophysical surveys were conducted). Because insufficient data existed to estimate the horizontal variogram, a 10:1 geometric anisotropy model was assumed (ratio of horizontal to vertical variogram range) based on prior experience at the site. The variogram sill was assumed to be isotropic. SISIM simulations of mud/sand distributions were conditioned on observations of ground-penetrating radar (GPR) attenuation as soft data. The soft data were derived by splitting the GPR attenuation data into five classes and determining the frequency of sand and mud in each class at boreholes (where data were co-located). The GPR attenuation data and the probability of sand being present given the GPR attenuation data are shown in Figure 5. The quantitative definition of the five classes is given in 03Table 3. Three indicator simulations representing a low, medium, and high proportion of mud are shown in Figure 6. Because the geophysical data strongly constrain the lithology, the proportion of mud varies only modestly across these three realizations from 12.4% to 17.8%. 04Table 4 provides summary statistics of GPR attenuation, hydraulic conductivity, and Fe(III) concentration data from boreholes by facies type. Hydraulic conductivity simulations were generated using SGSIM with the locally varying mean option. The locally varying means for the mud and sand facies were based on the median values of core permeability data as assigned to the two facies, 0.00000073 and 0.00278864 cm/s, respectively. In this manner, the simulations were conditioned indirectly on the geophysical soft data (through the mud/sand distributions), as well as directly to the core permeability measurements. Vertical vario-grams were calculated and modeled (Fig. 4) for the permeability residuals from the locally varying means (data value minus the locally varying mean), and again a 10:1 horizontal to vertical anisotropy ratio was assumed to define the horizontal variogram range. The three simulations of hydraulic conductivity that correspond to the indicator simulations shown in Figure 6 are shown in Figure 7. Simulations of Fe(III) were also generated using SGSIM with the locally varying mean option. The simulations were conditioned on the Fe(III) concentrations observed at the borehole sampling locations and the simulations of sand and mud based on the GPR attenuation data. The borehole data indicate that the Fe(III) concentrations are higher within sand units below the approximate level of the peat-rich zone, at an elevation of about 3 meters below mean sea level (see Figs. 3 and 8). Therefore, we used a two-part calibration, with a single mean for mud units wherever they were found, and separate calibrations for sand units above and below the −3 m elevation level. The three Fe(III) simulations paired with the indicator simulations discussed earlier are shown in Figure 9.

Method 3: (CoreData). A third set of correlated hydraulic conductivity and Fe(III) simulations (referred to here as “CoreData”) was generated based on core data alone, with no soft geophysical data constraint. The vertical variogram of the normal scores of the core hydraulic conductivity was estimated and modeled, and a 10:1 horizontal to vertical anisotropy ratio was assumed as in the previous methods. Three simulations of the hydraulic conductivity conditioned on core data alone are shown in Figure 10. The three simulations chosen for display correspond to the three simulations displayed in Figure 7 in that they were generated using the same random seed and following the same random path. Comparison of the corresponding simulations in both figures indicates that there are similarities between them based on the use of the same hard core data and the same random path. The observed differences between the two simulations result from the influence of the soft geophysical data on the hydraulic conductivity simulations shown in Figure 7. Correlated simulations of Fe(III) were generated using SGSIM with the collocated cokriging algorithm, using the suite of previously generated hydraulic conductivity simulations as soft data. That algorithm requires the specification of the correlation coefficient between the normal scores of the hard and soft data. We used a correlation coefficient of 0.6, which is approximately the correlation observed between the hydraulic conductivity and Fe(III) simulations that were generated conditioned on the indicator sand and mud simulations. The vertical variogram of the normal scores of the Fe(III) core data shows a pronounced trend due to the difference in Fe(III) concentrations seen in the upper and lower sand. The fitted vertical variogram has a range of 2 m, so the horizontal to vertical anisotropy ratio was assumed in this case to be 5:1, such that the horizontal range is 10 m, consistent with that used for the rest of the simulations. The three realizations of Fe(III) corresponding to the hydraulic conductivity simulations in Figure 10 are shown in Figure 11.

For each of the three methods, several realizations were generated and used as input to the flow and reactive transport model described in the following section. For method 1 (Chen et al., 2004), ten realizations (pairs of simulated Fe[III] and lithology fields) were selected randomly from a suite of 10,000 realizations previously generated by Chen et al. (2004). For method 2 (Geophys2), eleven realizations (pairs of Fe[III] and K fields) were selected by sorting a suite of 100 randomly generated realizations according to the proportion of mud facies and choosing realizations corresponding to the following percentiles: 2.5, 10, 20, 30, 40, 50, 60, 70, 80, 90, and 97.5. For method 3 (CoreData), eleven realizations corresponding to those from method 2 (in the sense that they used the same random seed and random path in the simulation algorithm) were selected.

Flow and Reactive Transport Model

The model domain (simulated aquifer cross section) is 11.9 m long and 5.4 m thick. This is representative of the scale of transport experiments conducted at the South Oyster site, and approximates the dimensions of a cross section along the centerline of the experimental flow cell parallel to the direction of groundwater flow. For the simulations here, a horizontal head gradient of 0.03 was imposed through constant head boundaries at both ends of the domain, and no-flow boundaries were imposed along the top and bottom of the domain (representing the phreatic water table and an underlying low-permeability aquitard unit, respectively). The imposed boundary conditions induce macroscopically uniform flow (from left to right in the orientation used for figures and animations). Groundwater flow and reactive transport were simulated using the HYDROGEOCHEM code (Yeh et al., 2004). Groundwater flow was simulated as steady state using the permeability distributions generated as described herein and the imposed boundary conditions. A spatially uniform porosity value of 0.34 was assigned based on inference from tracer tests conducted at the site (Scheibe et al., 2001). Bulk density of the aquifer solids was assumed to be 1.75 g/cm3. Relatively small values of longitudinal and transverse dispersivity (1.0 and 0.1 m, respectively) were assumed that are representative of local-scale microdispersivity in sandy porous media and consistent with the fact that physical heterogeneity at scales larger than 0.5 m is explicitly quantified in the model.

The reactive transport model was used to simulate the following hypothetical scenario: Step 1: The model domain was contaminated by aqueous-phase uranium (dissolved U[VI]) flowing into the model domain from an unspecified upgradient source at a uniform concentration of 30 μM. This uranium sorbed to natural iron-oxide surfaces in the model domain, gradually coming to a quasi-equilibrium state after ∼22 yr.

Step 2: Once the quasi-equilibrium state was reached in the model, a simulated bioremediation scenario was initiated by injecting a 10 mM solution of sodium acetate (soluble electron donor) at the upgradient boundary of the model domain. In the model, this is represented as a constant-concentration boundary at the left end of the model domain. The resulting microbial growth, acetate utilization, iron and uranium reduction, and transport of mobile species (bacteria, acetate, Fe[II], and U[VI]) was simulated over a 200 d period.

RESULTS AND DISCUSSION

The simulated results are presented as static and animated visualizations of the two-dimensional distributions of the state variables, and as summary information in terms of cumulative percent of uranium reduced as a function of time in each suite of simulations. Figures presented in this section are all in reference to the median mud content realization from the Geophys2 ensemble (center panel of Fig. 6). For purposes of visual comparison, animations are provided for one realization from each of the three suites (a randomly selected realization from the Chen et al. [2004] ensemble and realization #12 from both the Geophys2 and CoreData ensembles).

Figure 12 shows the simulated distribution of uranium (sorbed and aqueous phase) after the initial 22 yr loading period. A corresponding animation illustrating the development of this spatial distribution over the 22 yr period of uranium loading is provided in Animation 1. Comparable animations from the selected Chen2004 and CoreData realizations are provided in Animations 2 and 3, respectively. Note that the aqueous-phase (dissolved) U(VI) is distributed fairly uniformly at approximately the input concentration, whereas the solid-phase (sorbed) U(VI) is heterogeneously distributed and strongly associated with areas of high Fe(III) content (compare to Fig. 9, middle panel).

Figure 13 (upper left) shows the simulated distribution of sodium acetate (electron donor) after 200 d of biostimulation. Preferential delivery of the electron donor to areas of high hydraulic conductivity is evident; significant portions of the transect remain devoid of electron donor even 200 d after injection. In this model system, since hydraulic conductivity and initial Fe(III) content are positively correlated, the electron donor is delivered to portions of the aquifer suitable for iron-reducing bacterial activity. Figure 13 (upper right) shows the simulated distribution of aqueous Fe(II)–reduced iron generated by microbial iron reduction. The distribution is closely related to that of solid-phase Fe(III) (lower-left panel of Fig. 13), because the Fe(III) serves both as a source of reduced Fe(II) and as sites for sorption of reduced Fe(II) (lower-right panel of Fig. 13). Comparison of the lower-left panel of Figure 13 with the center panel of Figure 9 provides an indication of the amount of Fe(III) used during the 200 d period of bio-stimulation. A corresponding animation showing the development of the distributions shown in Figure 13 over the 200 d simulation period is provided in Animation 4. Comparable animations from the selected Chen2004 and CoreData realizations are provided in Animations 5 and 6, respectively.

Figure 14 shows the simulated spatial distributions of generated iron-reducing bacteria biomass at the end of the 200 d biostimulation, both solid-associated (lower-left panel) and aqueous (upper-left panel). Most of the biomass is associated with grain surfaces, and in particular Fe(III) solids. Cell concentrations are spatially heterogeneous (patchy), varying over two orders of magnitude within the simulated transect. The right-hand panels of Figure 14 show the simulated distributions of aqueous (top) and sorbed (bottom) U(VI) at the end of the 200 d biostimulation. The microbial activity as simulated is effective at reducing the majority of the contaminant within the time period simulated, but there remain small patches of moderate concentration, particularly of sorbed U(VI). A corresponding animation showing the development of the distributions shown in Figure 14 over the 200 d simulation period is provided in Animation 7. Comparable animations from the selected Chen2004 and CoreData realizations are provided in Animations 8 and 9, respectively.

Figure 15 (left-hand panel) shows a comparison of the percent of total U(VI) reduced as a function of time using the Chen2004 and Geo-phys2 ensembles of realizations. Each is presented in terms of the mean and upper and lower bounds over the set of realizations considered. Although there are some differences, the results are overall quite similar for the two ensembles, indicating that the two methods of incorporating the geophysical observations led to consistent predictions. For these two cases, the range of model predictions (representing uncertainty associated with modeled heterogeneity) is modest but significant. Because the scale of the numerical experiments is smaller than that of a typical field-scale remediation effort, the magnitudes of effective bioremediation time are probably not meaningful. However, the time to 90% remediation varies from ∼25 to 150 or more days, a ratio of six in effective treatment time. A potential variation factor of six in predicted treatment time would be significant in the context of design and cost estimation for a field bioremediation plan. We also note that at this particular site, there was observed a positive correlation between hydraulic conductivity and Fe(III) content. A positive correlation is favorable for bioremediation, since it leads to more effective delivery of electron donor to zones where contaminants preferentially reside. In many cases, however, there exists an inverse or negative correlation (e.g., Tompson et al., 1996). A negative correlation would render the system even more sensitive to the spatial distribution of physical and biogeochemical properties, and would likely greatly increase both the required time for bioremediation and the uncertainty in model predictions.

The right-hand panel of Figure 15 shows the results for individual realizations within the Geophys2 ensemble in decreasing order (top to bottom in the legend) of total proportion of mud facies. It can be seen from this plot that the proportion of mud facies is not the primary control on the time to 90% reduction. For example, the median case (14.8% mud) takes significantly longer to achieve 90% reduction than all the other realizations (∼150 d), whereas the 60th percentile case (14.5% mud) reaches 90% reduction most quickly of all the realizations (∼25 d). By comparing the animations of U(VI) reduction and the K and Fe(III) distributions associated with each realization, it is evident that the relatively slow remediation of realization 12 (the median case) is associated with the large low-permeability zone near the bottom left of the cross section (see Fig. 7, center panel). The 60th percentile case (not shown) contains few low-permeability zones outside the main mud layer. The smaller low-permeability zones near the bottom are more effectively penetrated during the U(VI) loading stage, and also contain more Fe(III), such that the initial concentration of U(VI) is higher than in the larger main mud layer. During the relatively rapid remediation stage, the bioremediation effectiveness is lower in these hydraulically isolated zones.

Figure 16 shows the percent of total U(VI) reduced as a function of time for the third ensemble of realizations (CoreData), conditioned to core data only. Comparison to Figure 15 shows that the predictions are very different, both in terms of the mean and the range of predicted treatment times. The mean curve for the CoreData case falls outside the upper bound of the Geophys2 case, and the range of predictions is relatively much smaller for the CoreData case. This indicates that incorporation of the geophysical observations as soft data in the geostatistical simulation methods significantly impacts model predictions. Although one might expect that more conditioning data would lead to a decrease in predictive uncertainty, here the opposite effect was observed. We believe this reflects the lack of ability of the geostatistical model to adequately represent the heterogeneous structure with only weak conditioning (i.e., large conceptual model error); the use of extensive and highly resolved geophysical observations provides much stronger conditioning and significantly improves the ability of the geostatistical methods to represent the complex heterogeneous structure of the aquifer. In this case, the geophysical observations were able to identify smaller low-permeability zones that were not inferred from core data alone. As discussed already, it appears to be the smaller, more isolated low-permeability regions that inhibit rapid U(VI) reduction, whereas the large continuous mud layer, which is represented in the core observations, does not play as significant a role because of its lower Fe(III) content and lower initial U(VI) content.

CONCLUSIONS

Physical and biogeochemical heterogeneity exerts significant control on the processes involved in the transport and fate of metals and radionuclides in the subsurface environment. In the example considered here, we demonstrated how diverse data (biogeochemical, hydrogeologic, and geophysical) can be integrated into a quantitative model of reactive transport in heterogeneous media. The numerical simulations presented here provide insight into the role that physical and chemical heterogeneity plays in biogeochemical reaction processes in a particular environment. There is a continued need for joint measurements of spatially distributed biological, chemical, and physical properties of subsurface media in a variety of environmental settings. Parallel laboratory and field experiments under well-characterized conditions are also needed to develop the understanding required for proper representation of heterogeneous (and often correlated) physical and biogeochemical aquifer properties in predictive models of reactive transport and bioremediation.

This research was supported by the U.S. Department of Energy Natural and Accelerated Bioremediation Research (NABIR) Program. Access to the South Oyster research site was granted by The Nature Conservancy, Virginia Coast Reserve.