Environ. Sci. Technol. 2008, 42, 2457–2463
Beyond the Rayleigh Equation: Reactive Transport Modeling of Isotope Fractionation Effects to Improve Quantification of Biodegradation B O R I S M . V A N B R E U K E L E N * ,† A N D HENNING PROMMER‡ Department of Hydrology and Geo-Environmental Sciences, Faculty of Earth and Life Sciences, VU University Amsterdam, De Boelelaan 1085, NL-1081 HV Amsterdam, The Netherlands, and CSIRO Land and Water, Private Bag No. 5, Wembley WA 6913, Australia
Received August 09, 2007. Revised manuscript received November 29, 2007. Accepted November 30, 2007.
The Rayleigh equation is commonly applied to evaluate the extent of degradation at contaminated sites for which compoundspecific isotope analysis (CSIA) data are available. However, it was shown recently that (i) the Rayleigh equation systematically underestimates the extent of biodegradation in physically heterogeneous systems, while (ii) it overestimates biodegradation if sorption-based carbon isotope fractionation is relevant. This paper further explores these two isotope effects not captured by the Rayleigh equation by means of a numerical modeling approach. The reactive multicomponent transport simulations show that the systematic underestimation is considerably larger for fringe-controlled and Monod-type degradation reactions than for previously assumed redoxinsensitive first-order degradation kinetics, while for the nonsteady state front portion of plumes, the Rayleigh equation may falsely indicate the occurrence of and/or overestimate biodegradation. The latter anomaly results from carbon isotope fractionation during sorption. It occurs for both supplycontrolled degradation at the plume fringe and slow, reactioncontrolled degradation inside the plume core. The numerical model approach enables a more accurate interpretation of CSIA data and thereby improves the quantification of biodegradation processes.
Introduction To evaluate the performance of natural attenuation or enhanced remediation of organic pollutants at contaminated sites, it is essential to determine to what extent concentration reductions in time or space are caused by degradation, since other processes such as sorption and dilution do not contribute to the reduction of the total contaminant mass. For various low molecular weight organic pollutants such as BTEX, naphthalene, MTBE, and chlorinated aliphatic hydrocarbons, compound-specific isotope analysis (CSIA) can be used to distinguish degradation from other concentration* Corresponding author phone: +31-20-5987393; fax: +31-205989940; e-mail:
[email protected]. † VU University Amsterdam. ‡ CSIRO Land and Water. 10.1021/es071981j CCC: $40.75
Published on Web 02/26/2008
2008 American Chemical Society
attenuating processes (1, 2). Stable isotope fractionation analysis is based on the fact that molecules containing a heavy isotope (e.g., 13C, 2H) degrade at relatively slower rates than those composed exclusively of the more abundant light isotopes (e.g., 12C, 1H). This results in an increase in the isotopic ratio (δ13C, δ2H) of the residual contaminant. Since other nondegradative processes such as dilution, sorption, volatilization, and nonaqueous phase liquid (NAPL) dissolution incur only minor isotope fractionation (1–4), any observed increases in isotope ratio are generally considered to reflect mass destruction. The extent of biodegradation, as expressed by a change of the isotope ratio, is usually quantified by application of the Rayleigh equation (5, 6): RSt RS0
) f (R-1)
(1)
The equation describes the change in isotope ratio (R) of an element in a molecule (S) during degradation (from t ) 0 to t ) tmax) as a function of the remaining fraction f of the molecule and the process-specific kinetic isotope fractionation factor, R, in a closed and perfectly mixed reservoir. The Rayleigh equation, therefore, does not necessarily apply to open systems such as those encountered under field conditions in natural porous media (7). Abe and Hunkeler (8), for example, predicted through modeling that physical heterogeneity causes the Rayleigh equation to slightly, but systematically, underestimate (K OC , and consequently, Rsorption > 1. To simulate faster, K OC the isotope fractionation effect of sorption in our numerical model, separate Kd values were used for each isotope. The Kd value for the light isotope N is as follows:
K dN )
1-n N × foc × Fsolid × K OC n
(12)
where n is the porosity of the aquifer matrix (dimensionless), foc is the fraction of organic carbon (dimensionless), and Fsolid is the specific density of the solid medium (M L-3). The Kd value for the heavy isotope N* can be calculated by combining eqs 11 and 12. Model Set-Up and Parameter Values Selected. A growing contaminant plume was simulated by applying a fixed concentration boundary at the upgradient end of the model domain, where the lower 5 m represented the source zone (contaminant ) 1 mM, nondegrading tracer ) 1 mM, CO2 )
0 mM, δ13Corg ) -30‰: [12Corg] ) 9.892177 × 10-1 mM and [13Corg] ) 1.078236 × 10-2 mM) and the upper 7 m represented pristine groundwater (contaminant ) 0 mM, nondegrading tracer ) 0 mM, CO2 ) 0.26 mM () air-aqueous equilibrium at 25 °C)). The model domain has a length of 240 m and a height of 12 m. The model domain was discretized by 240 columns (∆x ) 1 m) and 96 layers (∆y ) 0.125 m). A groundwater flow velocity of 20 m/year was enforced (porosity ) 0.3; hydraulic conductivity ) 10 m/day; hydraulic gradient ) 0.00164 m/m). The total simulation time was 10 years and the selected iteration frequency between flow/ transport and reaction step was between 1 and 10 days. A longitudinal dispersivity (RL) of 1 m and a vertical transverse dispersivity (RV) of 0.01 m were selected. Those values can be considered as representative for field conditions and appropriate for the spatial scale of the simulation problem (18). An effective diffusion coefficient of 3 × 10-10 m2/s was used. The following values for degradation rate parameters were selected: koxic ) 0.1–100/yr; CO2 ) 0.5 mg/L (19); kanoxic ) 0.7/yr; KI ) 0.001 mg/L (19); Vmax ) 0.2 mM/yr; and Ks ) 4.7 mg/L (20). A kinetic carbon isotopic enrichment factor, ) (R - 1) × 1000, of -2.5‰ was applied, which corresponds to the average value for benzene degradation (2). To simulate the carbon isotope fractionation effect of benzene sorption, a Kd value of 2.07 was selected (Koc, 67 L/kg (9); foc, 0.005 () 0.5%); Fsolid, 2.65 kg/L; porosity, 0.3; Fbulk, 1.86 kg/L), and a carbon isotope fractionation factor, Rsorption of 1.00017 was adopted from Kopinke et al. (9). Only measurable effects are presented; therefore, model results are shown for concentrations equal or greater than 10-5 of the source concentration, which approximately represents the range between common source concentrations and the lower analytical limits of CSIA. Comparing the Extents of Biodegradation Predicted by the Rayleigh Equation with Corresponding Numerical Model Predictions. The numerically simulated extent of biodegradation, BM (%), was determined by scaling the contaminant concentrations with the corresponding simulated concentrations of the nondegrading tracer with similar sorption behavior and source concentration as the contaminant: BM(%) )
[Ccontaminant] × 100 [Ctracer]
(13)
The extent of biodegradation as calculated by means of the Rayleigh equation (1, 2), BRayleigh (%), was calculated by applying the equation to the (carbon) isotope ratios predicted by the model:
( (
BRayleigh(%) ) 1 -
Raquifer Rsource
) ) 1000⁄
× 100
(14)
where Rsource and Raquifer are the isotopic ratios (N*/N) of the source and within the aquifer, respectively. Abe and Hunkeler (8) calculated the underestimation (θ) of the Rayleigh equation relative to the true model-predicted extent of biodegradation as follows:
(
θ) 1-
)
BRayleigh(%) × 100 BM(%)
(15)
The shortcoming of this determination of the underestimation, as they also remarked, is that it becomes 0% for extensive degrees of biodegradation when both BM and BRayleigh approach 100%, even if the remaining fraction differs by orders of magnitude between BM and BRayleigh. Therefore, a more representative measure of underestimation that will be used in the current paper is
(
θ) 1-
)
(
)
kRayleigh ln(fRayleigh) × 100 ) 1 × 100 kM ln(fM)
(16)
This expression was also applied by Abe and Hunkeler (8) to calculate the underestimation of the (assumed) first-order rate constant, since -kt ) ln(f ), where f is the fraction of contaminant remaining due to biodegradation.
Results and Discussion Isotope Fractionation in Subsurface Systems. The Rayleigh equation describes the redistribution of the isotopes of an element in a molecule that is undergoing degradation in a fully mixed and closed reservoir. Groundwater systems, however, represent open systems where concentration and isotope ratio gradients develop as a result of both transport and degradation processes. Simulations (presented in the Supporting Information) show that the Rayleigh equation only predicts the extent of biodegradation perfectly (θ ) 0, eq 16) in the absence of dispersion and diffusion. As explained in more detail in the Supporting Information, hydrodynamic dispersion continuously fades isotope ratio gradients formed by degradation processes, causing the Rayleigh equation to systematically underestimate the extent of degradation, such that θ > 0. These simulations confirm the findings of Abe and Hunkeler (8) that larger physical heterogeneity, represented in our simulations by a higher hydrodynamic dispersion coefficient, increases θ, as the fading intensity increases. The simulations also confirm that higher degradation rates increase θ, because concentration and isotope ratio gradients become larger, enhancing dispersive fluxes, and fading of the expressed isotope effect. The type of degradation kinetics strongly influences the spatial distribution of θ in contaminant plumes as will be shown in the subsequent sections. Core Degradation Solely. Some compounds, in particular strongly halogenated compounds like perchloroethylene, only degrade under anaerobic, highly reducing conditions that prevail in the plume core, while their degradation is inhibited at the aerobic plume fringe. Figures 1a-c show that θ is low (∼ 4%) in the core of a plume, and comparable to values determined by Abe and Hunkeler (8), who used a model that neglected the redox-dependency of the degradation rate. However, for a “core-contaminant”, θ increases considerably toward the fringe (up to 20–30%), while the isotopic shift diminishes, and the extent of degradation decreases. Note, that even in the pristine groundwater (e.g., nondegrading tracer fraction
