Environ. Sci. Technol. 2004, 38, 2673-2679
Modeling Expected Solute Concentration in Randomly Heterogeneous Flow Systems with Multicomponent Reactions MARIA E. MALMSTRO ¨ M , * ,† GEORGIA DESTOUNI,‡ AND PHILIPPE MARTINET§ Industrial Ecology, Chemical Engineering and Technology, Royal Institute of Technology (KTH), SE-100 44 Stockholm, Sweden, Physical Geography and Quaternary Geology, Stockholm University (SU), SE-10691 Stockholm, Sweden, and Land and Water Resources Engineering, Royal Institute of Technology (KTH), SE-100 44 Stockholm, Sweden
Many environmental problems require assessment of extensive reaction systems within natural subsurface flow systems exhibiting large physical and biogeochemical heterogeneity. We present an approach to couple stochastic advective-reactive modeling of physical solute transport (LaSAR) with the geochemical model PHREEQC for modeling solute concentrations in systems with variable flow velocity and multicomponent reactions. PHREEQC allows for general and flexible quantification of a multitude of linear and nonlinear geochemical processes, while LaSAR efficiently handles field-scale solute spreading in stochastic heterogeneous flow fields. The combined LaSARPHREEQC approach requires very modest computational efforts, thereby allowing a large number of reactive transport problems to be readily assessed and facilitating handling of quantifiable uncertainty in environmental model applications. Computational efficiency and explicit handling of field-scale dispersion without introduction of excessive fluid mixing that may impair model results are general advantages of the LaSAR compared with alternative solute transport modeling approaches. The LaSAR-PHREEQC approach is restricted to steady or unidirectional flow fields, and our specific application examples are limited to homogeneous reaction systems without local or transverse dispersion-diffusion, although these are not general methodological limitations. As a comprehensive application example, we simulate the spreading of acid mine drainage in a groundwater focusing on Zn2+ and including relevant, major-component geochemistry. Model results show that Zn2+ may be substantially attenuated by both sorption and precipitation, with flow heterogeneity greatly affecting expected solute concentrations downstream of the mine waste deposit in both cases.
Introduction Irregular and essentially random flow heterogeneity in subsurface flow systems implies that there is uncertainty in * Corresponding author phone: +46-8-7908745; fax: +46-87905034; e-mail:
[email protected]. † Chemical Engineering and Technology, Royal Institute of Technology. ‡ Stockholm University. § Land and Water Resources Engineering, Royal Institute of Technology. 10.1021/es030029d CCC: $27.50 Published on Web 03/31/2004
2004 American Chemical Society
model input and in resulting predictions of solute concentrations and mass fluxes, which is rationally handled by a stochastic modeling approach (e.g., refs 1-6). In such highly variable and uncertain subsurface flow systems, environmentally relevant reactive transport problems also require assessment of extensive reaction systems and complex geochemical conditions. One such environmentally relevant transport problem, which will be specifically exemplified in this paper, is the fate of acidic, metal laden drainage in the surroundings of sulfidic mining waste deposits, i.e., the spreading, retention, and attenuation of acid mine drainage (AMD) in the groundwater downstream of the AMD source. Other environmentally important problems that require assessment of relatively comprehensive reaction systems include, for example, the spreading of organic pollutants in groundwater, geochemistry in and around landfills, engineered-barrier treatment of pollutants, and in-situ remediation of polluted sites. The LaSAR (Lagrangian Stochastic Advection-Reaction) approach to modeling advective-reactive solute transport in an ensemble of streamtubes, with spatially/temporally and randomly variable characteristics along and/or among them, has been found generally useful for a variety of heterogeneous flow systems and a wide range of mass transfer, retention, and attenuation processes, including microbially mediated processes (see, e.g., refs 1-15). The majority of previous LaSAR modeling, however, has been limited to relatively simple, single- or few-component reaction systems, for which mass balance equations can be solved analytically (e.g., refs 1-6, 8, 9, 11, 13, 16). Where such closed-form solutions have been sought, reaction quantification has often involved linearization and/or restriction to low solute concentrations, allowing for the assumption of constant values of geochemical master variables (e.g. pH, Eh, and ionic strength). More complex, yet still simplified and specifically predefined transport-reaction problems have also been handled by applying the LaSAR approach to computationally reasonable, numerical coupled flow and transport models (e.g., ref 12). Most environmentally relevant reactive transport problems, however, require more general and flexible geochemical modeling tools able to handle considerably more complex and extensive reaction systems and geochemical conditions. Yabusaki et al. (10) took a step toward such increased geochemical generality and flexibity, by coupling a LaSAR type of transport modeling approach to their own developed geochemical reaction path module. Furthermore, Tompson et al. (15) applied the GIMRT reaction path model (17) coupled with a LaSAR type of transport modeling approach to the specific problem of groundwater contamination from underground nuclear tests. In this study, we take a further step to couple the wellestablished, general, flexible, and widely used geochemical model PHREEQC (18) with a LaSAR type of transport modeling approach. On the one side, the PHREEQC model is able to flexibly address the complexity of a wide range of multicomponent, linear and nonlinear, equilibrium and kinetic geochemical reactions, and their feedback mechanisms on geochemical master variables, however, with relatively limited physical transport representation (e.g., refs 19-22). On the other side, the LaSAR transport modeling approach is a powerful tool for handling highly irregular and essentially random subsurface flow heterogeneity, however, so far applied only to relatively limited (e.g., refs 1-6, 8, 9, 11, 13, 16), or not widely accessible and/or generally flexible numerical reaction models (10, 12, 15). The particular, concrete, and novel coupling presented here of the geochemicalVOL. 38, NO. 9, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2673
FIGURE 1. Schematic representation of the solute transport conceptualization according to the stochastic, Lagrangian formulation. Solutes undergo reactive transport along individual streamlines between the injection plane (IP) and the control plane (CP). Due to flow heterogeneity, the streamlines have randomly variable advective travel (residence) time along and among them. Insert: Schematic representation of the PHREEQC conceptual model for reactive transport along an individual streamline. ly widely used PHREEQC model with the LaSAR transport modeling approach represents an important step across disciplinary boundaries, yielding a combined modeling tool that is considerably extended and more useful for environmental applications than each of the PHREEQC and previous LaSAR based model tools separately. The specific objectives of this study are to (1) develop the conceptual and computational methodology for coupling the LaSAR and PHREEQC modeling approaches for flow systems where local geochemical equilibrium can be applied as a first step toward the development of a more general equilibrium and kinetic reaction system representation; (2) test the model performance against the analytical solution for a simplified, coupled reaction-transport problem; and (3) investigate the model applicability and the effect of flow heterogeneity on model results for an environmentally relevant, comprehensive reaction system.
Modeling Methodology Stochastic Description of the Reactive Transport Problem. Figure 1 shows a schematic representation of the reactive transport of solute in a heterogeneous groundwater system according to a stochastic Lagrangian travel time formulation (3). The solute enters the groundwater system (transport domain) of interest through an injection plane (IP) that is normal to the mean water flow direction. The injection plane may represent, for example, the soil surface, or the boundary between the unsaturated and the saturated subsurface zone, or that between a waste deposit and the surrounding, natural geological formation. We consider a case where the injection of the solute is continuous for t > 0 and uniform over the IP. In a probabilistic sense, the solute can take any of several paths, each conceptualized as a streamline, to a control plane (CP) that is situated normal to the mean flow direction at an arbitrary distance x ) L from the IP. The CP represents our environmental point of interest, e.g. a compliance boundary, a subsurface-surface water interface, or a boundary to a sensitive domain of the groundwater system. Along each of the streamlines, (bio)geochemical reactions may occur. Processes such as biochemical degradation (organic pollutants) or radioactive decay (unstable isotopes), which destroy or transform the pollutant as well as interaction of the solute with the solid phase, e.g. sorption and precipitation, cause attenuation and retention of the solute along the streamline. The governing equation for reactive transport (e.g., ref 7), according to the Lagrangian formulation is
∂C ∂C + +R)0 ∂t ∂τ 2674
9
(1)
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 9, 2004
where C denotes total solute concentration in the aqueous phase, t denotes the time, τ denotes the travel (or residence) time of water and inert tracer, and R represents the rate of removal of solute from the aqueous phase due to (bio)geochemical or other processes. The specific water discharge rate, q, relates the travel time and length coordinates such that τ ) (xθ/q), with θ being the water content and (q/θ) the inert tracer velocity. The mean water and inert tracer travel time between IP and CP, τj, can, for instance, be estimated from the total volumetric flow rate of water, Q, and the flow volume in the system, V, as τj ) (V/Q). The flow volume may be estimated as V ) ALθ h from the average cross sectional area, A, the length, L, along the mean flow direction and the average water content, θ h , of the groundwater system of interest. Subsurface heterogeneity implies that each of the streamlines may have spatially and temporally variable characteristics along and among them. In order for the streamlines to also constitute pathlines for solute transport, the water flow must be either essentially steady or unidirectional. In reality, spatial subsurface variability may occur in physical as well as (bio)geochemical characteristics. In this study, however, we focus on the physical heterogeneity, causing spatially variable groundwater velocities, which in turn yield randomly variable solute travel times among the different streamlines, and assume that other subsurface properties, such as, e.g., temperature and chemical characteristics, are constant in space. We follow the common approach to neglect local and transverse dispersion within mobile water, thus reducing the full transport problem to coupled advection-reaction in a system of streamlines with randomly variable solute travel time; this simplification yields computational efficiency, but the LaSAR approach is not generally limited to such conditions (e.g. ref 14). Thus, macrodispersion is handled by the stochastic flow field, which is advantageous over the corresponding use of large dispersion coefficients in Eulerian approaches, which may impair model results by introducing unrealistic reactions due to implied extensive fluid mixing, the effect of which may be especially problematic for nonlinear reaction systems. Modeling of a Single Streamline. We use the PHREEQCcode (18) to quantify the geochemical processes and plugflow, advective transport occurring along an individual streamline. In the PHREEQC-model (version 2), the reactions can be conceptualized as occurring in a series of homogeneous, sequential cells, with identical volumes, ∆V, representing the spatial resolution of the system. Transport is modeled through stepwise movement of solute from one cell to the next in the series, where a step represents the temporal resolution of the model, ∆t. Reactions are decoupled from transport and quantified through mass-balance considerations in connection with geochemical equilibrium and kinetic expressions, which are evaluated for each cell between the transport steps. The single-streamline model is hence a 1-D reactive transport model, using a split-operator approach, in which advective transport is treated with explicit finite difference that is forward in time and space. For this application, we conceptualize a series of cells to comprise a streamtube, with total water volume VT from the IP to the CP, centered around one individual streamline, which by definition carries a constant fraction of the total water flow, ∆Q, from the IP to the CP (see insert to Figure 1). Advection couples the temporal and spatial scales of the model, as, for mass conservation reasons, ∆t ) ∆τ, where ∆τ is the water residence time in a cell, given by ∆τ ) (∆V/∆Q). The nth cell in a PHREEQC-simulation is hence associated with a streamline with a total travel time τn ) n∆τ from the IP to the CP, which is located at the downstream boundary of the cell. Similarly, the mth forward shift of cells corresponds to elapsed time tm ) m∆τ after solute injection.
Field-Scale, Average Solute Concentration. The expected field-scale resident concentration, quantified in the following, corresponds to the mean value of locally determined concentrations over the CP surface. This expected concentration, C h (t,x), is obtained as
C h (t,x) )
∫ C(t,τ;x)g(τ;x)dτ ∞
0
(2a)
in which C(t,τ;x)is the local concentration at time t and at CP position x along the mean flow direction associated with a streamline with water travel time τ(x), and g(τ;x) is the probability density function (pdf) of travel times τ(x) for all streamlines over the whole flow system (Figure 1). We interpret, as explained in the previous section, C(m,n) from a PHREEQC simulation as discrete evaluations of C(tm ) m∆τ,τn ) n∆τ). We obtain the expected field-scale concentration C h (t,x) by approximation of eq 2a as n)ntot
C h (t,x) ≈
∑ C(m,n)g(τ ;x)∆τ n
(2b)
n)1
The numerical summation in eq 2b was implemented in a simple MATLAB procedure. The separation of reaction process and transport into individual modules provides a convenient and computationally efficient model approach; however, it also implies limitations on the applicability of the combined model, such as that the flow field must be essentially steady or unidirectional and independent of the reaction system. Moreover, it is assumed that no mixing occurs between the streamtubes. In this study, we exemplify the use of eq 2 by a log-normal water travel time distribution, which has been indicated as relevant for groundwater solute transport in previous studies (e.g. refs 21 and 23)
g(τn;x) )
1 1 -1/2((lnτn-µ)/σ)2 e τn σx2π
(3)
where µ ) lnτ and σ are the average and the standard deviation of the log travel time, respectively, with µ ) lnτj σ2/2. However, any other type of pdf, for instance as obtained from field-scale tracer tests or simulated flow fields, can be implemented as easily. Note that for systems with processes at geochemical equilibrium, the PHREEQC output, C(m,n) is independent of the temporal and spatial scales of the problem, and, hence, the result of a single PHREEQC simulation can be used to assess several locations of the CP (x-values) or different temporal resolutions (∆τ-values). However, as the travel time to the CP and its statistics depend on the CP location (e.g., ref 23), the travel time pdf, g(τn, x), has to be evaluated separately for each CP scenario.
Illustrative Examples We apply the proposed modeling methodology for quantifying the spreading of Acid Mine Drainage (AMD) in an anoxic, heterogeneous groundwater surrounding a sulfidic mine waste deposit. Rather than quantifying the AMD source term, which has already been discussed in several previous publications (e.g., refs 25-29), we here aim at assessing the pollutant spreading, retention, and attenuation in the groundwater downstream of the AMD source, an assessment that is, for instance, necessary in environmental risk assessments. We assume that the oxygen flux, which would be dominated by vertical diffusion and thus be essentially normal to the mean groundwater flow direction, is negligible. This assumption is applicable for deeper parts of groundwater systems, where molecular oxygen has been consumed, for example, by overlying organic material or reduced minerals
TABLE 1. Impoundment Groundwater Composition Used in the Model Scenariosa component 2-
SO4 Fe(tot)b Zn2+ Ca2+
concn [M]
component
concn [M]
0.10 0.080 0.0060 0.0050
Mg2+
0.011 0.0020 4.9
Al pH
a Data based on monthly averages in Impoundment 1, Kristineberg, Sweden, prior to remediation, as reported by Ekstav and Qvarfort (30). b For this modeling purpose assumed to be Fe2+.
(compare e.g., ref 20). This AMD transport and natural attenuation problem provides a good example of a transportreaction system where it is necessary to assess multiple components, relatively high concentration of solutes, and nonlinear geochemical processes, a system that has so far been typically hard to quantify within the LaSAR modeling approach, using analytical or problem-specific numerical methods. We focus here on the transport and attenuation of Zn2+, a metal ion that commonly occurs at high concentration in AMD, and which, in part due to its single oxidation state, has a relatively simple geochemistry. We provide thereby an example case rather than a comprehensive model of a real system, to exemplify the application and potential usefulness of the suggested modeling methodology. Hence, we use a few simplifying, but not necessary, assumptions regarding dominant processes, process parametrization, and parameter values. In a more comprehensive assessment of a real system, the selection of the dominant transport and geochemical processes to be included in the model and their quantifications must be based on detailed field and laboratory observations. AMD Plume from a Sulfidic Mine Waste Deposit: Model Scenarios. We use a composition of the AMD that corresponds to the groundwater in mill tailings impoundment 1 in Kristineberg, northern Sweden, prior to remediation (Table 1, ref 30). For simplicity, we assume that reactions in the aqueous phase and between solutes and solid phases can be quantified through geochemical equilibrium expressions and that calcite is the only mineral providing buffer capacity to the groundwater system. The calcite content was arbitrarily set to 0.3 mol per liter of soil solution, corresponding to ∼0.5 vol % of the solids. We consider a partial mineral succession that is consistent with the one suggested by Blowes and Ptacek (31). At the plume front, calcite dissolves and buffers the pH around 7.9. Behind the calcite depletion front, the pH is buffered by dissolution of secondary mineral phases at consecutively lower pH levels. We expect the geochemistry of Zn2+ in the surroundings of the deposit to be dominated by aqueous speciation and retention/attenuation by sorption, precipitation, or coprecipitation. Zinc sorption onto ferric phases has been reported to be extensive above pH ∼5, with an increasing fraction of the Zn2+ being sorbed with increasing pH (e.g., refs 32 and 33). The pH dependence of Zn2+ sorption implies that sorption will be intimately coupled to the major ion chemistry in the AMD plume, which governs the pH of the soil solution. Zinc ions are fairly soluble, however, according to PHREEQC speciation calculations using the MINTEQ-database (34), several solid phases may become oversaturated and thus have the potential to precipitate in the plume. Although calculations of saturation indices indicate the most stable phases, these do often not occur or do not exert solubility control on dissolved metal ion concentrations in natural systems, due to kinetic limitations on precipitation/dissolution reactions. Of the zinc containing phases indicated to have a potential to precipitate, Nordstrom and Alpers (35) rate only smithsonite, ZnCO3(s), as “likely to control metal VOL. 38, NO. 9, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2675
concentrations” in AMD-environments. Extensive coprecipitation of Zn2+ in secondary ferric precipitates has been observed in AMD environments (36) but is not considered here. In this paper, we assess two possible scenarios where Zn2+ attenuation is dominated by either (A) sorption onto hydrous ferric oxide (HFO), assumed to be evenly distributed in the groundwater system with an arbitrarily set surface site concentration of 0.01 mol per liter of soil solution or (B) precipitation of secondary smithsonite. In both cases, we assume that the plume is spreading in an anoxic, but nonreducing groundwater environment, such that oxidation of ferrous iron, and subsequent formation of secondary ferric precipitates as well as reduction of sulfate, and subsequent formation of insoluble sulfide minerals can be omitted from the geochemical model. We assume an average specific groundwater discharge of q j ) (Q/A) ) 3.3 × 10-7m s-1, corresponding to a water flow of 1 × 10-7 m3 s-1 per square meter of IP/CP and an average water content of 0.3. We address the time evolution of expected solute concentration on a CP at 20 m distance from the IP (corresponding to τj ) 694 days (1.9 years)) and investigate an interval of standard deviation in the log travel time corresponding to relatively low (σ ) 0.1) to high (σ ) 1.0) degrees of subsurface flow heterogeneity. The PHREEQC implementation of the geochemical model scenarios, including 9 components, ∼50 species, and 4-5 solid phases, is detailed in Appendix 1 in the Supporting Information. Scenario A: Simplified Reaction System. To verify the numerical model against an analytical solution to the coupled reactive transport problem, we assess first a simplified reaction system for our scenario A, in which we only consider aqueous speciation of components and sorption at a constant pH of 6.5. As shown in Appendix 2 in the Supporting Information, this coupled reactive transport system can, after appropriate simplifications, be expressed as tot,aq tot,aq ∂CZn ∂NZn2+ ∂CZn 2+ 2+ + + )0 ∂t ∂τ ∂t
(4)
∂Ctot,aq ∂Ctot,aq + )0 ∂t ∂τ
(5)
tot,aq CZn C 2- + KZn(SO4)23-CSO42-2) (6) 2+ ) CZn2+(1 + KZnSO 4(aq) SO4
(
tot,aq CMg 2+ tot,aq + CSO 2- ) CSO 2- 1 + KMgSO 4 4 4(aq)1 + C SO42-
)
tot,aq CFe 2+ KFeSO4(aq) (7) 1 + CSO42-
RCZn2+ NZn2+ ) β + KsCZn2+
(8)
where N denotes the sorbed concentration of Zn2+, R and β denote empirical sorption isotherm constants, Ks is the sorption equilibrium constant, and Ki is the equilibrium constant for formation of aqueous species i. Eq 5 regards nonreactive solute transport; in our simplified reaction system, this applies to sulfate, magnesium, and iron. The problem defined by eqs 4-8 has, as shown in Appendix 3 in the Supporting Information, an analytical solution tot,aq tot,aq tot,aq Ctot,aq (t,τ) ) Cj,in + (Cj,0 - Cj,in )H(t - τ) j
(j ) SO42-, Mg2+, Fe2+) (9) tot,aq CZn 2+ (t,τ)
)
tot,aq Cin,Zn 2+
+
tot,aq (C0,Zn 2+
tot,aq - Cin,Zn 2+)H(τs(t) - τ) (10)
for initial and boundary conditions as Ctot,aq ) Ctot,aq for x > in 0, t ) 0 and Ctot,aq ) Ctot,aq for x ) 0, t > 0; and where H is 0 the Heaviside step function. In eq 10, τs ) [1 + (N0 - Nin)/ 2676
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 9, 2004
FIGURE 2. Modeled breakthrough curves (BTCs), C h (t,x), at L ) 20 m from the IP (τj ) 1.9 years) for the simplified model scenario A (sorption of Zn2+ onto HFO at pH 6.5) and different travel time 2+ variabilities. (a) SO24 and (b) Zn . Symbols show results from the analytical solution (eqs 9 and 10), and lines show results of the LaSAR-PHREEQC numerical model. Full line and triangles - deterministic model. Dotted line and diamonds - low travel time variability (σ ) 0.1). Dashed-dot line and circles - high travel time variability (σ ) 1.0). In panel b, results to the left are for a situation where Zn2+ is not attenuated, and results to the right are for sorption onto a HFO, with a surface site concentration of 0.01 M. tot,aq tot,aq -1t, where N ) RγCtot,aq /(β + K γCtot,aq ), (C0,Zn 0 s 2+ - Cin,Zn2+)] 0,Zn2+ 0,Zn2+ tot,aq tot,aq Nin ) RγCin,Zn2+/(β + KsγCin,Zn 2+), and where γ ) [1 + KZnSO4(aq) / / 2 -1 with C / 2C SO 2 + KZn(SO4)2 (C SO 2-) ] SO42 as obtained from 4 4 tot,aq / tot,aq / C0,SO42- ) C SO42-(1 + KMgSO4(aq)(C0,Mg2+/(1 + C SO + 2-) 4 tot,aq / KFeSO4(aq)C0,Fe /(1 + C )]. 2+ SO42Figure 2a,b shows the results of the analytical solution (symbols) along with the coupled LaSAR-PHREEQC model (lines) for different travel time variability. For comparison, the result of a scenario where we assume that Zn2+ is not attenuated is also shown in Figure 2b. For this case, Zn2+ arrives simultaneously as SO24 , which is modeled as a nonreactive tracer. Figures 2a,b show good agreement between the analytical solution and the numerical LaSAR-PHREEQC model, hence, verifying the numerical implementation of the latter, combined model. Moreover, the fact that the deterministic PHREEQC results in Figure 2 (solid lines) closely resemble step-functions indicates that numerical error within the PHREEQC model is small and does thus not contribute to the modeled spreading of the expected breakthrough curves (BTCs), C h (t,x). In general, numerical errors, including rounding and truncation errors, may occur in any numerical reactive-transport code such as PHREEQC. To ensure that numerical error does not considerably contribute to model results, the deterministic BTCs, i.e., BTCs obtained from only PHREEQC simulations, from successively finer model discretization may be compared. Finally, Figure 2 also demonstrates that consideration of flow heterogeneity, as expressed through the travel time variability, may considerably affect the spreading of the BTC. This effect is even more pronounced in the case where sorption attenuates Zn2+, as has previously been reported for sorption reactions in general (16).
FIGURE 3. Expected BTC at L ) 20 m from the IP (τj ) 1.9 years) from the LaSAR-PHREEQC model for the full model scenario A (sorption of Zn2+ onto HFO at pH as determined by major component chemistry) and different travel time variabilities. (a) Fe2+; (b) SO24 ; (c) pH (see text for interpretation); and (d) Zn2+. Full line - deterministic model. Dotted line - low travel time variability (σ ) 0.1). Dashed-dot line - high travel time variability (σ ) 1.0). Scenario A: Full Reaction System. Figure 3 shows the result of the LaSAR-PHREEQC model for the more comprehensive model scenario A and for different degrees of travel time variability. In this scenario, Mg2+ acts as a nonreactive component and shows similar behavior as SO24 in the simplified reaction system (BTC for Mg2+ not shown). At the plume front, calcite is unstable with respect to siderite, FeCO3(s), and gypsum, CaSO4 × 2H2O(s), and dissolved sulfate and ferrous iron are both attenuated through precipitation of these secondary phases
CaCO3(s) + Fe2+ + SO42- + 2H2O S FeCO3(s) + CaSO4 × 2H2O(s) (11) Iron is almost fully attenuated by the reaction in eq 11 until all calcite has been consumed, which delays the arrival of iron compared to magnesium (Figure 3a). As for magnesium, increased flow heterogeneity considerably increases the early arrival of iron and the overall spreading of its expected BTC. In the presence of calcite, the reaction in eq 11 is in our case limited by ferrous iron, due to the stoichiometry of the process and the lower concentration of ferrous iron than of excess sulfate in the impoundment water. Excess sulfate, CSO 2- ) 4 tot,imp. tot,imp. -1, is mobile and breaks CSO C ) 0.02 mol L 2Fe2+ 4 through at t ) τj, causing a steplike appearance of the sulfate BTC at low measure values of the travel time variability (Figure 3b). The small sudden increase in sulfate concentration at time around 5 years in the deterministic model is associated with the breakthrough of Zn2+ at this point, in conjunction with higher solubility of sulfate in the presence of Zn2+. After depletion of calcite up to the CP (at t ≈ 8.5 years), the sulfate concentration increases considerably and exceeds that of the impoundment water somewhat. This is due to redissolution of precipitated gypsum, which is kept at a low level by the high concentration of sulfate in the impoundment water and which continues throughout the simulation. The dramatic differences between the expected sulfate BTC for the simplified and full reaction systems (compare Figures 2a and 3b), with considerably delayed breakthrough of sulfate
in the full reaction system at low flow variability, are due to the reaction in eq 11 and demonstrate the necessity of considering the more comprehensive reaction system. Note that for higher flow variability, differences between the simplified and full reaction systems are smaller, indicating that effects of geochemical processes on the BTC can be masked in highly heterogeneous flow fields. Greater flow heterogeneity results in smoothing of the steps also for the full reaction system and increased early arrival of sulfate, an effect that may be considerable and that also must be accounted for in highly variable flow systems. For our model scenario, the unpolluted groundwater has a pH of ∼8.2. Upon arrival of acidic impoundment water (aluminum acidity is dominating at the modestly low pH), calcite buffers pH at around pH 7.9 but is simultaneously consumed in the favor of siderite (reaction in eq 11). After depletion of calcite, dissolution of siderite, FeCO3(s) + H + S Fe2+ + HCO3 , buffers pH at around pH 6.4 in each individual streamline. The expected BTC for pH shown in Figure 3c should be interpreted as related to the numerical average of locally measured pH values in samples extracted from the CP. As pH is the logarithm of an intensity variable, this average will be different from the pH value measured in a physical mixture of the samples. Nevertheless, for low travel time variability, the expected BTC shows the steplike feature associated with consecutive buffering by calcite and siderite, with increasing variability increasingly smoothing out these steps. For the highest considered variability (σ)1.0), the steps have totally disappeared, and the BTC is smooth. Figure 3d, finally, shows the expected BTC for Zn2+. The strong similarity to the corresponding BTC for the simplified reaction system (compare Figure 2b) is due to that major ion chemistry buffers pH above the sorption edge of Zn2+, which suppresses geochemical differences between the cases. Scenario B: Full Reaction System. Because the major ion chemistry is little affected by Zn2+, BTCs for Mg2+, Fe2+, SO24 , and pH are similar between model scenarios A (Figure 3a-c) and B (not shown). Figure 4a shows the expected BTC for Zn2+ for model scenario B, in which precipitation of VOL. 38, NO. 9, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2677
Discussion
FIGURE 4. LaSAR-PHREEQC model results for scenario B (precipitation of smithsonite) and different travel time variabilities. (a) Expected BTC for Zn2+ at L ) 20 m from the IP (τj ) 1.9 years); (b) Zn2+ concentration profile at 2, 5, and 8 years. Full line deterministic model. Dotted line - low travel time variability (σ ) 0.1). Dashed-dot line - high travel time variability (σ ) 1.0). Note that the vertical axis has been truncated and that the position of the control plane is marked by the line labeled “CP” in panel b. smithsonite is the Zn2+-attenuating process. The peak appearance of the deterministic BTC is due to that, although Zn2+ precipitates in the presence of calcite, Zn2+ + CaCO3(s) S ZnCO3(s) + Ca2+, smithsonite is unstable with respect to siderite and redissolves at the calcite depletion front
ZnCO3(s) + Fe2+ S Zn2+ + FeCO3(s)
(12)
thereby remobilizing Zn2+. Figure 4a shows a high sensitivity of the Zn2+ concentration to flow heterogeneity, which greatly lowers the maximum load (peak height) but increases the duration of increased load levels (peak width) even for small travel time variability. At the highest variability, the peak in fact disappears totally from the expected BTC. Moreover, increasing variability increases the early arrival of Zn2+ such that a considerable concentration occurs long before t ) τj. Figure 4b shows expected profiles for dissolved Zn2+ concentration vs distance from the IP at time 2, 5, and 8 years after injection. The profiles were obtained by moving the CP along the profile and evaluating eq 2b for discrete values of x. The broadening of the Zn2+ peak in the deterministic model (solid line) as the plume moves along the profile is due to the accumulating amount of Zn2+ that is attenuated before calcite is depleted in conjunction with a fixed Zn2+ concentration during the redissolution of precipitated smithsonite (equilibrium associated with the reaction in eq 12). Model results for σ ) 0.1 (dotted line) show that travel time variability causes further broadening of the Zn2+ peak and that this broadening increases as the peak moves along the profile. For the highest considered travel time variability (σ) 1.0; dashed-dot line), the profiles are smooth with a pronounced early arrival of Zn2+. 2678
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 9, 2004
We have coupled the general, flexible, and geochemically widely used PHREEQC model for multicomponent geochemical reactions with a LaSAR modeling approach to reactive solute transport in randomly heterogeneous flow systems. This coupling uses similar methodology to that presented by Yabusaki et al. (10) for combining LaSAR modeling with their own developed geochemical reaction path module. The present LaSAR- PHREEQC coupling, however, provides a potentially more generally accessible and flexible tool for quantification of extensive multicomponent reaction systems, with high concentrations of dissolved constituents and complex interactions between solutes and geochemical master variables, within randomly heterogeneous subsurface flow systems. In common with previous LaSAR model applications of e.g. Yabusaki et al. (10), also the present LaSAR-PHREEQC transport-reaction modeling tool is restricted to flow systems that are essentially steady or unidirectional and independent of the reaction system. In the coupled model, PHREEQC provides a general and flexible means to quantify (bio)geochemical processes in a detailed, process-based, building-from-first principles manner, and the LaSAR modeling approach provides a general and flexible means to account for the highly irregular and essentially random spatial heterogeneity of subsurface flow systems, in terms of relevant site- and problem-specific nonreactive solute travel time distributions. A general advantage of LaSAR, as compared to alternative transport modeling approaches, is the computational efficiency gained by decoupling advection dominated solute transport from the geochemical reactions, which is especially useful in modeling of reactive transport in heterogeneous subsurface systems. Moreover, the LaSAR transport modeling approach handles explicitly the commonly large subsurface heterogeneity effect on field-scale solute spreading, through the averaging of different streamtube compositions over the relevant physical travel time pdf. Alternative transport modeling approaches, such as the use of the Eulerian advection-dispersion equation, may instead handle physical fieldscale solute spreading by use of a large dispersion coefficient, implying mixing of large water volumes and possibly introducing unrealistic reaction effects, which may impair model results, especially for nonlinear reaction systems. Any (bio)geochemical reaction path module can in principle be coupled with a LaSAR type of transport modeling approach, following a similar methodology as Yabusaki et al. and others (see e.g., refs 10 and 15). A particular advantage of the LaSAR-PHREEQC coupling presented here, however, is the wide geochemical use of the component model PHREEQC, which may make the present coupled tool more generally flexible and accessible across disciplinary boundaries, for instance to geochemists, who have so far not commonly addressed and handled heterogeneous flow, or to contaminant hydrologists, who may widely handle heterogeneous flow systems but not commonly coupled with complex geochemical systems. For the specific objectives and application examples of this paper, we have assessed geochemically homogeneous systems using equilibrium geochemical process expressions.The LaSAR approach, however, is not generally limited to thermodynamic conditions (e.g. refs 1, 2, 5, 6, 10, 12); for the specific LaSAR-PHREEQC coupling presented here, consideration of (bio)geochemical kinetics is achieved through appropriate matching of the time step length in the PHREEQC module and the travel time discretization in the LaSAR module. Furthermore, assessment of variable chemistry between the streamtubes, such as e.g. describing physical and geochemical variability in a stratified aquifer, may readily be assessed through appropriate averaging over multiple, individually homogeneous streamtubes, but with different
reaction path chemistry between them. Assessment of deterministically known, variable chemistry along individual streamtubes, such as representing e.g. zones with different known reaction path chemistry, would require book-keeping and a combination of a sufficient number of different LaSARPHREEQC realizations. For a simplified transport-reaction problem, the LaSARPHREEQC model showed excellent agreement with an analytical solution, thereby verifying the proposed numerical solution. We applied the proposed LaSAR-PHREEQC model to example scenarios of AMD spreading, retention, and attenuation in the groundwater downstream from the AMD source. The two different investigated scenarios considered Zn2+ attenuation through either sorption onto HFO or precipitation as smithsonite. In both scenarios, complex interactions between solutes and solid phases in a system of moving reaction fronts governed major ion chemistry, which determine e.g. pH, the most important geochemical variable for sorption and solubility reactions, as well as the Zn2+ attenuation itself. The results for both example scenarios were promising, as we, after a minimum of model optimization regarding discretization choices, were able to derive expected BTCs for all components at the CP with very modest computational efforts. As an example, we also derived expected concentration profiles for Zn2+ between the IP and CP for different times. Model results indicated that Zn2+ may be substantially attenuated by both sorption and precipitation under the considered geochemical conditions, with flow heterogeneity greatly affecting expected solute concentrations in both cases. Based on the demonstrated usefulness of the coupled LaSAR-PHREEQC model in assessing the groundwater spreading and natural attenuation of AMD, as one example system, we find the approach highly promising for further testing on other complex, problem- and site-specific modeling of coupled transport-multicomponent-reaction systems. As PHREEQC is a standard tool for geochemical modelers and the handling of stochastic solute travel times is mathematically and numerically relatively straightforward in the presented coupled approach, we believe that the LaSARPHREEQC approach will not only extend the number of reactive transport problems that can be readily assessed but also facilitate the use of stochastic modeling and thereby the explicit handling of uncertainty in standard model applications as well as in geochemical research modeling tasks.
Acknowledgments
Dr. Sten Berglund, Land and Water Resources Engineering, Royal Institute of Technology, Stockholm, Sweden is acknowledged for valuable discussions on LaSAR modeling. Three anonymous reviewers and Dr. M. J. Small are acknowledged for comments during the review process. This research was partly funded by the Swedish Foundation of Strategic Research (MISTRA) through the research program “Mitigation of the environmental impact from mining waste” (MiMi). Georgia Destouni also acknowledges financial support from FORMAS, the Swedish Research Council for Environmental, Agricultural and Spatial Planning.
Supporting Information Available The PHREEQC implementation of base-case model scenarios and the mathematical definition as well as the analytical solution to the simplified reactive transport problem available as Appendices. This material is available free of charge via the Internet at http://pubs.acs.org.
Literature Cited (1) Cvetkovic, V.; Shapiro, A. M. Water Resour. Res. 1990, 26, 20572067. (2) Destouni, G.; Cvetkovic, V. Water Resour. Res. 1991, 27, 13151325. (3) Cvetkovic, V.; Dagan, G. J. Fluid Mech. 1994, 265, 189-215.
(4) Destouni, G.; Graham, W. D. Water Resour. Res. 1995, 31, 19351944. (5) Ginn, T. R.; Simmons, C. S.; Wood, B. D. Water Resour. Res. 1995, 31, 2689-2700. (6) Simmons, C. S.; Ginn, T. R.; Wood, B. D. Water Resour. Res. 1995, 31, 2675-2688. (7) Dagan, G.; Cvetkovic, V. Proc. R. Soc. London, Ser. A. 1996, 452, 285-301. (8) Cvetkovic, V.; Dagan, G. Proc. R. Soc. London, Ser. A 1996, 452, 303-328. (9) Eriksson, N.; Destouni, G. Water Resour. Res. 1997, 33, 471483. (10) Yabusaki, S. B.; Steefel, C. I.; Wood, B. D. J. Contam. Hydrol. 1998, 30, 299-331. (11) Cvetkovic, V. Phys. Fluids 2000, 12, 2279-2294. (12) Kaluarachchi, J. J.; Cvetkovic, V.; Berglund, S. J. Contam. Hydrol. 2000, 41, 335-365. (13) Painter, S.; Cvetkovic, V.; Turner, D. Ground Water 2001, 39, 326-338. (14) Fiori, A.; Berglund, S.; Cvetkovic, V.; Dagan, G. Water Resour. Res. 2002, 38, art. no. 1137. (15) Tompson, A. F. B.; Bruton, C. J.; Pawloski, G. A.; Smith, D. K.; Bourcier, W. L.; Shumaker, D. E.; Kersting, A. B.; Carle, S. F.; Maxwell, R. M. Environ. Geol. 2002, 42, 235-247. (16) Berglund, S.; Cvetkovic, V. Water Resour. Res. 1996, 32, 23-32. (17) Steefel, C. I.; Yabusaki, S. B. OS3D/GIMRT Software for Modeling Multicomponent- Multidimensional Reactive Transport: User Manual & Programmer’s guide; Pacific Northern National Laboratory Report PNL-11166: Richland, WA, U.S.A., 1996. (18) Parkhurst, D. L.; Appelo, C. A. J. User’s guide to PHREEQC (version 2); USGS Water-Resources investigations report 99-4259; USGS: Denver, CO, U.S.A., 1999. (19) Appelo C. A. J.; Verweij, E.; Scha¨fer, H. Appl. Geochem. 1998, 13, 257-268. (20) Brown, J. G.; Bassett R. L.; Glynn, P. D. J. Hydrol. 1998, 209, 225-250. (21) van Breukelen, B. M.; Appelo, C. A. J.; Olsthoorn, T. N. J. Hydrol. 1998, 209, 281-296. (22) Postma, D.; Appelo, C. A. J. Geochim. Cosmochim. Acta 2000, 64, 1237-1247. (23) Destouni, G.; Simic, E.; Graham, W. Water Resour. Res. 2001, 37, 2303-2308. (24) Demmy G.; Berglund, S.; Graham, W. Water Resour. Res. 1999, 35, 1965-1973. (25) Elberling, B.; Nicholson, R. V. Water Resour. Res. 1996, 32, 17731784. (26) Wunderly, M. D.; Blowes, D. W.; Frind, E. O.; Ptacek, C. J. Water Resour. Res. 1996, 32, 3173-3187. (27) Bain, J. G.; Blowes D. W.; Robertson W. D.; Frind E. O. J. Contam. Hydrol. 2000, 41, 23-47. (28) Mayer, K. U.; Frind, E. O.; Blowes, D. W. Water. Resour. Res. 2003, 38, art. no. 1174. (29) Salmon, S. U.; Malmstro¨m, M. E. Appl. Geochem. 2004, 19, 1-17. (30) Ekstav, A.; Qvarfort, U. Metallbalans Kristineberg; Kvarta¨rgeologiska avd., Uppsala Universitet: Uppsala, Sweden (in Swedish), 1989. (31) Blowes, D. W.; Ptacek, C. J. In Env. Geochem of Sulphide MineWastes; Jambor, J. L., Blowes, D. W., Eds.; Min. Assoc. Can.: Nepean, Canada, 1994; pp 271-292. (32) Dzombak, D. A.; Morel, F. M. M. Surface complexation modeling; John Wiley and Sons: New York, U.S.A., 1990. (33) Jo¨nsson, J.; Lo¨vgren, L. In Proceedings of Securing the Future: International Conference on Mining and the Environment, 2001, Skellefteå, Sweden, June 25-July 1 2001; pp 320-326. (34) Allison, J. D.; Brown, D. S.; Novo-Gradac, K. J. MINTEQA2/ PRODEFA2, A geochemical assessment model for environmental systems; Version 3.0 User’s Manual; EPA/600/3-91/021; U.S. Environmental Protection Agency: Athens, GA, 1994. (35) Nordstrom, D. K.; Alpers, C. N. In Environmental geochemistry of mineral deposits, Reviews in Economic Geology 6A; Plumlee, G. S., Logsdon, M. J., Eds.; Soc. Ec. Geologists: Chelsea, U.S.A., 1999; pp 133-160. (36) Lo¨vgren, L. Personal communication, 2002.
Received for review February 17, 2003. Revised manuscript received January 20, 2004. Accepted January 23, 2004. ES030029D VOL. 38, NO. 9, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2679
