- Email: [email protected]
Implementation of biofilm permeability models for mineral reactions in saturated porous media
ARTICLE IN PRESS
Computers & Geosciences 31 (2005) 968–977 www.elsevier.com/locate/cageo
Implementation of biofilm permeability models for mineral reactions in saturated porous media$ Vicky L. Freedman, K. Prasad Saripalli, Diana H. Bacon, Philip D. Meyer Pacific Northwest National Laboratory, P.O. Box 999, MSIN K9-36, Richland, Washington 99352, USA Received 8 July 2004; received in revised form 14 February 2005; accepted 22 February 2005
Abstract An approach based on continuous biofilm models is proposed for modeling permeability changes due to mineral precipitation and dissolution in saturated porous media. In contrast to the biofilm approach, implementation of the film depositional models within a reactive transport code requires a time-dependent calculation of the mineral films in the pore space. Two different methods for this calculation are investigated. The first method assumes a direct relationship between changes in mineral radii (i.e., surface area) and changes in the pore space. In the second method, an effective change in pore radii is calculated based on the relationship between permeability and grain size. Porous media permeability is determined by coupling the film permeability models (Mualem and Childs and Collis–George) to a volumetric model that incorporates both mineral density and reactive surface area. Results from single mineral dissolution and single mineral precipitation simulations provide reasonable estimates of permeability, though they predict smaller permeability changes relative to the Kozeny and Carmen model. However, a comparison of experimental and simulated data show that the Mualem film model is the only one that can replicate the oscillations in permeability that occur as a result of simultaneous dissolution and precipitation reactions occurring within the porous media. r 2005 Published by Elsevier Ltd. Keywords: Permeability; Geochemical models; Hydrogeology; Numerical models
1. Introduction There is an increasing interest in reactive transport models that integrate porosity and permeability changes that occur as a result of mineral precipitation and dissolution in a porous medium. Such models are needed for accurately predicting contaminant fate and transport through porous media undergoing hydraulic property $ Code on server at http://www.iamg.org/CGEditor/index.htm Corresponding author. Tel.: +1 509 372 4067; fax: +1 509 372 6089. E-mail address: [email protected] (V.L. Freedman).
0098-3004/$ - see front matter r 2005 Published by Elsevier Ltd. doi:10.1016/j.cageo.2005.02.003
changes. Several models exist in the literature. In the reactive transport literature, for example, a discrete approach is usually taken to describe changes in permeability. Typically, the surface area of a mineral is used as input to the model, or calculated from the given volume fraction and initial sphere radius for the mineral (Ortoleva et al., 1987; Sanford and Konikow, 1989; Carnahan, 1990; Steefel and Lasaga, 1994; LeGallo et al., 1998). When mineral-fluid reactions occur, either the mineral radius or its surface area is recalculated accordingly. In the biological reactive solute transport literature, a continuous biofilm approach has been developed for describing hydraulic property changes due to biomass
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
accumulation in porous media. In this method, it is assumed that the solid particles that make up the porous media are uniformly covered with a biofilm. This approach assumes that the porosity and permeability reduction is thereby uniform within the pore space (Taylor et al., 1990; Cunningham et al., 1991). This uniform change in pore space may also occur for mineral reactions. For example, uniform dissolution and precipitation is consistent with surface complexation models that formulate reaction rates based on the surface concentration of specific chemical groups (Yamamoto et al., 1996; Duckworth and Martin, 2003). Observations of glass corrosion experiments for low-level vitrified waste have also shown that as the glass dissolves over time, the secondary mineral precipitates form as a uniform coating over the glass grains (McGrail et al., 1999). The goal of this work is to develop a film depositional model applicable to mineral precipitation and dissolution reactions. The model is developed to introduce an alternative permeability model that can be coupled to a reactive transport simulator for cases where uniform precipitation and dissolution may occur. Although Taylor et al. (1990) also presented a film-based porosity model, early testing demonstrated that unless the film model accounted for mineral density, it could not be adapted for use in a reactive transport model. Hence, this paper focuses on the permeability model development, and describes the film model theory and the main program functions. The model is implemented utilizing two separate methods for calculating the relative change in pore radii. An example problem involving simultaneous dissolution and precipitation is demonstrated. Model results are compared to laboratory experiments to ensure that some degree of generality can be attached to model conclusions and interpretations.
2. Film model theory The cut-and-random-rejoin approach for describing permeability assumes that porous media contains pores of different radii that are randomly distributed in space. The hydraulic conductivity of porous media depends on the number of interconnected pores and their geometric configuration. Permeability can be expressed in terms of the pore-size distribution function, which Taylor et al. (1990) define as a power function given as b r l1 f ðrÞ ¼ , (1) rmax rmax where f ðrÞ ¼ the volume of full pores of radii between r and r þ dr per unit volume of medium, r ¼ the pore radius ½L, rmax ¼ maximum pore radius ½L within the pore-size distribution ½L, and b ¼ a dimensionless
969
constant. The parameter l represents an index of the pore-size distribution, and whose range is from 0olp1, with larger values being associated with more homogeneous pore sizes. In the section that follows, Eq. (1) is used in the development of the permeability models, which closely follows the presentation given in Taylor et al. (1990). As outlined below, the development of a mineral film model differs from the biofilm model of Taylor et al. (1990) inasmuch as it accounts for a increase in the radii defined by the pore-size distribution.
3. Film model development The porosity ðnÞ of porous media can be expressed as a function of the pore-size distribution ðf ðrÞÞ and is given as Z rmax n¼ f ðrÞoðrÞ dr, (2) rmin
where r ¼ the pore radius ½L and o ¼ the weighting factor for radius r, which is assumed to be unity. Maximum and minimum pore radii are defined by rmax ½L and rmin ½L, respectively. Thus, utilizing the pore-size distribution function defined in (1), the initial diffusive porosity of a clean porous medium ðnÞ is " #1 b rmin l n¼ 1 . (3) l rmax In integral form, the Mualem (1976) and Childs and Collis–George (1950) permeabilities (km and kccg , respectively) are Z rmax 2 tn0:5 km ¼ rf ðrÞ dr , (4) 8 rmin kccg ¼
1 4
Z
rmax
r2 ½n yðrÞf ðrÞ dr,
(5)
rmin
where t ¼ tortuosity ½L, and yðrÞ ¼ the unsaturated moisture content. Substituting (1) into both permeability formulations yields ( " #)2 tb2 r2max n0:5 1 rmin 1þl km ¼ , (6) 1 8 rmax 1þl
kccg
" ( # b2 r2max 1 rmin 2þl 1 ¼ 1 4ln rmax 2þl 2 þ 2l " #) rmin 2þ2l . 1 rmax
ð7Þ
Eqs. (3), (6), and (7) describe a porous medium in the absence of a mineral film. When mineral precipitation occurs, it is assumed that a film of uniform thickness, Lf
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
970
½L, uniformly coats the internal surfaces in the porous media. Mineral film growth completely seals pores having a radius less than Lf , and reduces the radius of those pores larger than Lf . For example, pores with a new radius r, previously had a radius r þ Lf , and their cylindrical volume (assuming long tubular pores) is changed by the ratio of their volumes ½ðphÞr2 =½ðphÞðr þ Lf Þ2 , where h represents the length of the pore tube. If f p ðrÞ describes the pore-size distribution after mineral precipitation and dissolution occurs, then f ðrÞ and f p ðrÞ are related by f p ðrÞ ¼
r2 f ðr þ Lf Þ, ðr þ Lf Þ2
rmin Lf orormax Lf .
ð8Þ
By substituting f p ðrÞ for f ðrÞ in Eq. (2), the porosity resulting from mineral precipitation is Z rmax Lf r2 np ¼ f ðr þ Lf Þ dr, (9) ðr þ Lf Þ2 rmin where rmin ¼ max½rmin Lf ; 0. The permeabilities for mineral precipitation are obtained in the same manner. Substituting f p ðrÞ for f ðrÞ in (4), and (5) yields "Z #2 rmax Lf tn0:5 r3 p f ðr þ L Þ dr , (10) kmp ¼ f 8 ðr þ Lf Þ2 rmin
kccgp ¼
1 4
Z
rmax Lf
½np yr ðrÞ rmin
r4 f ðr þ Lf Þ dr. ðr þ Lf Þ2 (11)
For mineral dissolution, it is assumed that the internal surfaces of the porous media decrease by a uniform thickness. Thus, pores with a new radius r, previously had a radius r Lf , and their volume is increased by the ratio of volumes r2 =ðr Lf Þ2 . If f d ðrÞ describes the voids affected by a decrease in a uniform mineral film thickness, then f d ðrÞ and f ðrÞ are related by f d ðrÞ ¼
r2 f ðr Lf Þ, ðr Lf Þ2
rmin þ Lf orormax þ Lf .
ð12Þ
By substituting f d ðrÞ for f ðrÞ in Eq. (2), the permeabilities resulting from mineral dissolution are "Z #2 rmax þLf tn0:5 r3 d kmd ¼ f ðr Lf Þ dr , (13) 8 ðr Lf Þ2 rmin þLf
kccgd ¼
1 4
Z
rmax þLf
½nd yr ðrÞ rmin þLf
r4 f ðr Lf Þ dr. ðr Lf Þ2 (14)
The film-based permeability equations for Mualem [Eqs. (10) and (13)], and Childs and Collis–George [Eqs. (11) and (14)] are used for calculating permeability in the reactive transport simulator.
4. Program implementation The mineral film permeability model presented in this paper was first developed as a stand-alone program. Required input included porosity (n), tortuosity (t), initial pore radius (r), as well as the minimum and maximum radii of the pore-size distribution (rmin and rmax , respectively). In this work, it is assumed that tortuosity is constant. Also required as input to the program were the incremental changes in pore radii (Lf ). A critical component in the implementation of the film models in a reactive transport simulator was in determining a time-dependent method for the incremental changes in the pore space. For the stand-alone program, and in Taylor et al. (1990), it was assumed that the film thickness was static. In fact, the stand-alone program is capable of calculating both porosity and permeability based on film models that utilize a static film thickness. However, the stand-alone model does not consider mineral surface areas and densities, which is an important aspect for implementation with a reactive transport simulator. Hence, porous media porosity was determined using a volumetric model that considered both density and surface area. For permeability calculations, the mineral film thickness needed to dynamically change as a function of time-dependent reaction processes. To this end, mineral volume changes due to mineral precipitation were interpreted as a change in the film thickness that uniformly altered the radii of all pores in the pore-size distribution. For dissolution, the change in volume also occurred uniformly on the solid-phase surface. Thus, a negative change in the mineral film thickness was translated into an increase in the pore radii throughout the entire pore-size distribution. To calculate changes in hydraulic properties within the reactive transport simulator, porous media permeability was determined by coupling the film permeability models to a volumetric model that incorporated both mineral density and reactive surface area. Once new mineral volumes were determined, the film permeability models required the change in pore radii (Lf ) to be dynamically calculated from the precipitation and dissolution reactions. Two separate methods for determining the change in pore radii were used. In the first method, it was assumed that the cumulative pore radii changes were equal to the cumulative changes in mineral radii: X X t1 Ltf ¼ Rtm Rt1 (15) m þ Lf ,
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
where t ¼ the current timestep, t 1 ¼ the previous timestep, and Rm ¼ mineral radius ½L, which is not related to the pore radius (r) initially described in Eq. (1). Alternatively, an effective Lf was determined based on a relationship between permeability and grain size. Like the film models of Taylor et al. (1990), the Kozeny and Carmen permeability model is based on conduit flow, which assumes steady, laminar flow through bundles of capillary tubes. Because it is difficult in natural systems to determine the tube radius, the Kozeny and Carmen model was modified to relate permeability to the mineral grain size. An earlier version of the Kozeny and Carmen equation [in Pape et al., 2000] attempts to incorporate the effects of effective hydraulic pore radius and the length of parallel tubes k¼
nr2eff
. (16) 8t Although the effective hydraulic radius cannot be directly measured, it can be determined from the following relationship (Pape et al., 2000) reff ¼
2Rm tn . 3ð1 nÞ
(17)
By including an expression for the effective porosity, Eq. (17) accounts for the fact that only a fraction of the pores are interconnected and participate in the flow process. Thus, in the second method, the cumulative Lf ðeff Þ was determined as 2Rm tn t 2Rm tn t1 t Lf ðeff Þ ¼ þ Lft1 (18) ðeff Þ . 3ð1 nÞ 3ð1 nÞ Eq. (18) is also consistent with the Mualem permeability model in as much as it accounts for the tortuosity of the flow path. In Taylor et al. (1990), analytical solutions were developed to calculate permeability. In this work, numerical solutions for the integrals in Eqs. (10), (11), (13) and (14) were generated using the mathematical software Maple (by Maplesoft).
5. Reactive transport simulator The film permeability models were coupled to the reactive mass transport and flow simulator, Subsurface Transport Over Reactive Multiple phases (STORM) (Bacon et al., 2004). In STORM, flow of fluid phases was computed from Darcy’s law. The resulting flow fields were then used to solve conservation equations for reactive aqueous phase transport. Transport of phase components was computed from Fick’s law. Once flow and transport calculations were complete, porosity and permeability calculations were performed. For brevity’s sake, the governing equations describing flow, transport
971
and chemical reactions are not presented in this document, but are described in detail in Bacon et al. (2004). However, a brief description of the porosity and permeability models used in STORM is provided below. The two governing equations describing saturated flow and advective-dispersive transport are discretized into algebraic equations using the control volume method (Patankar, 1980). The nonlinear equation sets are converted into linear form utilizing iterative Newton–Raphson method. Details of the numerical solution methods can be found in Bacon et al. (2004). A primary assumption associated with the solute mass conservation equation is that solute concentrations are infinitely dilute. This assumption implies that fluid properties are independent of solute concentrations, which allows the solute conservation equations to be solved independently of the flow equation. Thus, once the flow equation is solved, the chemical reaction equations and the solute transport equations are solved simultaneously. Reactions that only involve aqueous species are assumed to be at equilibrium. Mineral precipitation and dissolution reactions are treated kinetically. Reaction rates in STORM are based on transition state theory and consider the forward rate and the equilibrium constant 2 3 , jnij j jnij j Y Y W j ¼ AmðjÞ kmðjÞ 4 ðgi ci Þ f j ðgi ci Þ f j ðK rðjÞ Þ1=f 5, nij o0
nij 40
(19) where W j ¼ the reaction rate for reaction j [mol/MT], Am ðjÞ ¼ the effective mineral reaction area ½L2 ], kj ¼ the rate constant [1/t], K r ðjÞ ¼ the equilibrium constant, gi ¼ the activity coefficient of species i, ci ¼ the concentration [moles/M], and f i ¼ the factor of the reaction order. The stoichiometric coefficient of species i in reaction j is represented by nij . The bars j j represent the absolute value and the sign of the reaction is positive for precipitation and negative for dissolution. Once the flow, transport and chemical reaction equations are solved, porosity and permeability updates are then calculated by the STORM texture solver, and are used as input values for flow in the time stepping routine. 5.1. Volumetric equations STORM calculates volume fractions by assuming that mineral grains are spherical in shape, and occupy a given volume based on the number of grains and the mineral radius. Assuming that the number of grains remains constant, changes in porosity are calculated based on a change in the mineral radius. STORM describes the relationship between the mineral volume fraction ðV m ½L3 =L3 Þ and the porosity
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
972
Table 1 Parameter values for numerical simulations
of the porous medium as n ¼ 1:0
Nm X
V m.
(20)
Property
Value
Initial porosity (n) Initial permeability (k) Maximum pore radius (rmax )
0.50 2:92 1010 m2 1:125 103 m (Mualem) 1:441 104 m (Childs and Collis–George) 0.0 m 3.0 1:86 102 m 5:84 1010 m
m¼1
where N m is the total number of minerals. Assuming a spherical geometry, the number of grains (nm ) present for each mineral is defined as nm ¼
3V m 4pR3m
(21)
and the effective reactive surface area (Am ) is Am ¼ 4pR2m nm .
(22)
By assuming that a dissolution/precipitation reaction is a surface-limited kinetic reaction, the change rate of the volume fraction is qV m qRm ¼ 4pR2m nm . qt qt
(23)
Minimum pore radius (rmin ) Pore-size distribution (l) Tortuosity (t) Kozeny and Carmen fitting coefficient ðC o Þ Calcite initial surface area (SA) Quartz initial surface area (SA) Darcy flux ðqÞ Domain length ðLÞ Grid discretization ðDzÞ Total simulation time (T)
1:2084 102 m2 =m3 2:5000 105 m2 =m3 2:5 105 m=s 2.0 m 0.1 m 5.0 days
5.2. Permeability equations Permeability in the STORM simulator is calculated by a relationship defined by Kozeny and Carmen. This permeability equation is empirical and requires the input of fitting parameters. For example, the Kozeny and Carmen ðk ½L2 Þ equation is given as (Bear, 1972) k ¼ Co
n3 , ð1 nÞ2 A2s
(24)
where As is the specific surface area and C o ¼ an empirical fitting coefficient.
6. Model testing Evaluation of the film depositional models was performed by comparing general trends in hydraulic property changes with the Kozeny and Carmen model [Eq. (24)] already implemented in STORM. By adjusting the fitting parameters, the initial permeabilities were approximately equal between the STORM and film models at the beginning of the simulations. To calculate changes in permeability, the models required n, l, rmin , and rmax as input. Because hydraulic conductance largely depends on the connectivity of larger pore sizes, rmax essentially determined the magnitude of permeability. Since the permeability was insensitive to the lower limit of integration, the minimum radius was assumed to have a value close to zero rmin ¼ 1 1015 . These calculations were independent of the initial mineral radius [Eq. (22)] required for input in STORM.
Although several example problems were run to test the depositional models, results from only a few representative simulations are presented in the sections that follow.
6.1. Precipitation A one-dimensional, saturated flow simulation involving the precipitation of calcite ðCaCO3 Þ at 25 1C was used to test the film models. Initially, the column was slightly saturated with respect to calcite. Input parameters appear in Table 1. For flow, a specified flux boundary condition was assigned at the top of the column, and a free gradient boundary at the bottom. An initial conditions boundary type was assigned at both the top and bottom boundaries for solute concentrations. This boundary type fixed the concentration to the initial concentration value at the node adjacent to the boundary surface, and was invariant with time (Bacon et al., 2004). At time greater than zero, a water flux entered the column from the top, which maintained super saturation with respect to calcite. As the calcite precipitated, incremental changes in the pore radii occurred, causing a porosity and permeability reduction. The precipitation of calcite in the column was simulated for 5 days. Results from this simulation are presented below. In Figs. 1 and 2, results for surface area changes and findings from the Mualem film permeability model are plotted with STORM’s implementation of the Kozeny
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
973
Distance from Outlet (m)
2.0 Kozeny and Carmen C and C-G (Lf) C and C-G (Lf-eef)
1.5
1.0
0.5
0.0 0.6
0.8
1.0
1.2
K/Ko
Fig. 1. Relative changes in reactive surface area for calcite precipitation simulation.
Fig. 3. Permeability reduction profiles for the Kozeny and Carmen (grain packing) and Childs and Collis–George film models.
Distance from Outlet (m)
2.0 Kozeny and Carmen Mualem (Lf) Mualem (Lf-eef)
1.6 1.2 0.8 0.4 0.0 0.6
0.8
1.0
1.2
K/Ko
Fig. 2. Permeability reduction profiles for the Kozeny and Carmen (grain packing) and Mualem film models.
and Carmen model, respectively. Results of the Mualem permeability calculations showed that even though surface area changes were the same for all three permeability models, the magnitude of the Kozeny and Carmen permeability reduction is greater than both of the implementations of the Mualem film model. At the top of the column, the relative change in surface area was only 13%, yet the permeability reduction was 31% for the Kozeny and Carmen model, and only 26% and 11% for the Mualem models ðLf and Lf ðeff Þ , respectively). The shape of the permeability reduction profiles was fairly similar for the Kozeny and Carmen and Mualem ðLf Þ model, but the effective ðLf ðeff Þ Þ implementation was less steep as the permeability reduction predicted over the length of the column was less than the other two models. Shown in Fig. 3 are predictions from the implementation of the Childs and Collis–George film permeability model with the Kozeny and Carmen equation. Similar to the Mualem permeability model, the Childs and Collis–George equation experienced a larger permeabil-
ity reduction due to calcite precipitation when Lf was calculated as in Eq. (15). However, this difference was quite small. At the top of the column, the Lf implementation of the Childs and Collis–George model predicted a 5% permeability reduction, whereas the Lf ðeff Þ implementation predicted a 4% reduction. Relative to the Kozeny and Carmen model, the Childs and Collis–George film models estimate a smaller permeability reduction. As is evident from the results shown in Figs. 2 and 3, predictions made by the Kozeny and Carmen model showed a more dramatic decrease in permeability throughout the column relative to all of the film model implementations. Similar behavior was demonstrated in the permeability changes for mineral dissolution presented in the following section. 6.2. Dissolution Inorganic reactions may cause an increase in the porosity and permeability of the porous media by increasing the pore radii due to solid phase dissolution. To test the film growth models in STORM, quartz dissolution simulations for a one-dimensional column were conducted at 100 1C. The column was undersaturated with respect to quartz, but otherwise had the same physical setup as the calcite precipitation test case (Table 1). Shown in Fig. 4 is the relative changes in surface area profile, and in Figs. 5 and 6 the permeability profiles for the Kozeny and Carmen model as well as the Mualem and Childs and Collis–George film models, respectively. Fig. 5 shows that the two lines representing the permeability enhancement profiles for each implementation of the Mualem model coincide, as the predictions were nearly identical. Although the maximum change in surface area was only 5% at the top of the column, both implementations predict only a 4% increase in the
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
974
Fig. 4. Relative changes in reactive surface area for quartz dissolution simulation.
6.3. Column experiment with simultaneous precipitation and dissolution
Distance from Outlet (m)
2.0 1.6 Kozeny and Carmen Mualem (Lf) Mualem (Lf-eef)
1.2 0.8 0.4 0.0 0.8
1.0
1.2
1.4
1.6
K/Ko
Fig. 5. Permeability enhancement profiles for the Kozeny and Carmen (grain packing) and Mualem film models.
Distance from Outlet (m)
2.0 1.6 1.2 Kozeny and Carmen C and C-G (Lf)
0.8
C and C-G (Lf-eff) 0.4 0.0 0.8
two lines representing the Childs and Collis–George permeability profiles coincide in Fig. 6, but relative to the Kozeny and Carmen model, the two implementations predict a smaller permeability enhancement (7%). Although all of the permeability profiles presented for both mineral precipitation and dissolution reactions yield plausible results, the film models estimated smaller permeability changes to the Kozeny and Carmen model. However, comparisons of film-based permeability estimates to the Kozeny and Carmen relation were only used for reference, not as an absolute scale of performance. To this end, permeability predictions for a column experiment are examined in the following section.
1.0
1.2 K/Ko
1.4
1.6
Fig. 6. Permeability enhancement profiles for the Kozeny and Carmen (grain packing) and Childs and Collis–George film models.
permeability at the top of the column compared to a 15% increase predicted by Kozeny and Carmen. A similar result occurred for the predictions made by the Childs and Collis–George models. Not only do the
Matrix acidizing is a process which uses acids to dissolve minerals near a wellbore to increase permeability and enhance oil recovery. However, some of the dissolution products can precipitate as the acid reacts with geologic materials. Commonly, iron hydroxides are observed to plug pore throats and decrease permeability. The simulations presented in this section are based on the experiments of Rege and Fogler (1989). In their investigation, Rege and Fogler (1989) acidized a limestone core using 0.75 M iron chloride with a pH of one. Dissolution of the calcite ðCaCO3 ) by the acid caused the solution pH to increase and ferric hydroxide ðFeðOHÞ3 ) to precipitate. Variations in the permeability occurred as a result of the simultaneous precipitation and dissolution reactions occurring within the core. Both the film models and the Kozeny and Carmen equation are compared to the experimental permeability ratio reported by Rege and Fogler (1989). The onedimensional column measures 0.089 m, and is discretized into 20 equally spaced nodes. Since the initial porosity and permeability of the experimental limestone core are not reported, the porosity was assumed to be equivalent to the value used in the previous test simulations reported in Table 1. The initial permeability of the medium was assumed to be 5:8177 1010 m2 . The kinetic reaction rate for FeðOHÞ3 ðsÞ was an order of magnitude lower ð1 107 mol=m2 s) than that of the CaCO3 ð1 106 mol=m2 s). The flow rate for the simulations is 1 108 m3 =s, and is equivalent to the one used in the acidizing matrix experiments of Rege and Fogler (1989). A reaction front in the simulations is created by injecting a 0.075 M solution of ferric chloride ðFeCl3 Þ with a pH of 1.7 into a one-dimensional column containing calcite. The solution in the column is initially undersaturated with respect to calcite and has a pH of 7.9. This setup differs from the pH and strength of concentration (1.0 and 0.75 M, respectively) of the
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
injected solution in the experiments. However, concentration limits imposed by both the gauss-seidel iterative solver in STORM, and the b-dot ion activity model (Helgeson, 1969), which is valid for ionic strengths p0:1 M, prevented an exact replication of the Rege and Fogler (1989) experiments. Hence, the simulations presented here only attempt to replicate the permeability trends reported by Rege and Fogler (1989), and not produce an exact replication of the chemistry of their system. 6.3.1. Comparison of results The changes in reactive surface area for both calcite and FeðOHÞ3 ðamÞ for the top of the column are found in Fig. 7a, and the calculated permeability responses for this work, as well as the one found by Rege and Fogler (1989) are shown in Fig. 7b. This includes the simulated permeability ratios for the Kozeny and Carmen, Mualem and Childs and Collis–George models at 0.0134 m from the top of the column. In all cases, acid is consumed through the dissolution of calcite, which is shown to decrease its surface area throughout the simulation. As the pH increases, the iron in solution precipitates as iron hydroxide. Because the rate of iron hydroxide precipitation exceeds the rate of calcite dissolution, the permeability initially decreased in the column. For the Mualem predictions, both the Lf and Lf ðeff Þ methods were used to calculate changes in the pore space. Although an attempt was made to use the Lf method for the Childs and Collis–George model, the permeability decreased to zero and convergence was not obtained. Hence, for the Childs and Collis–George permeability predictions shown in Fig. 7, only the Lf ðeff Þ method is shown. In the Rege and Fogler (1989) experiments, fluctuations in the permeability response are observed that are absent in all but the Mualem Lf simulation. Rege and Fogler (1989) explain that these oscillations were due to the deposition of the precipitate ahead of the dissolution front. Initially, the precipitate formed in flow paths with higher flow rates. As deposition of the precipitate occurred, there was an increase in the flow resistance, which caused flow to be diverted to flow paths with lower rates of flow. As the smaller flow paths lengthened, flow resistance decreased. This caused the sudden increase in permeability at 8 pore volumes. In addition to the diverted flow paths, the oscillations and excursions in the permeability profiles of the (Rege and Fogler, 1989) experiments and in the Mualem permeability profile, could also be a result of the cyclical nature of precipitation–dissolution reaction events in the medium. The rate equation for a typical precipitation–dissolution reaction shows that the rate of reaction is proportional to ð1 Q=KÞ, where Q is the activity product of aqueous species involved in reaction, and K is
975
the equilibrium constant for reaction. Although precipitate formation is thermodynamically favorable at or above saturation values (i.e., when Q=K41Þ, precipitation will not generally occur until a significant degree of super-saturation is reached. Deposition of the precipitated mass in the porous medium is a faster, or more sudden event compared to the slow dissolution process. Once such precipitation occurs, the pore fluids are depleted of the solute mass and the system becomes thermodynamically more favorable for dissolution, which starts a fresh precipitation–dissolution cycle. Since the permeability of the porous medium is sensitive to the thickness of the precipitate film [as shown in Eqs. (4) and (5)], it follows that the cyclic or oscillatory precipitation–dissolution behavior would likely lead to a similar oscillatory change in permeability as seen in Fig. 7. Using the Lf method, the Mualem model was able to capture this behavior, due to the sensitivity of the method to mineral volume changes. Not only was the Mualem model able to replicate the oscillations in permeability, but it also predicted a general increase in permeability at 7 pore volumes, similar to the increase in permeability observed at 8 pore volumes. Because of the differences in the concentration and pH of the injected fluid, the simulated permeability decrease occurred less rapidly than the one observed experimentally. With the Lf ðeff Þ method, however, the translation of the mineral volume changes to changes in pore radii were too small to capture this behavior. In fact, a monotonic permeability reduction was predicted by both the Mualem ðLf ðeff Þ ) and Childs and Collis–George ðLf ðeff Þ ) models because of the insensitivity of the method to the increases in the pore radii. A less dramatic, monotonic permeability reduction was predicted by the Kozeny and Carmen model, but the Kozeny and Carmen model also failed to capture the oscillatory behavior of the permeability in the system. In contrast to the Mualem and Kozeny and Carmen models, the Childs and Collis–George film model predicted dramatic changes in permeability. Between 4 and 9 pore volumes, the permeability decreased by nearly 100%. This decrease in permeability was not due to the lengthening of smaller flow paths, rather fluctuations in film thicknesses. Because the Childs and Collis–George film model does not include an explicit relationship to the porosity, changes in permeability are entirely dependent on film thickness. As the ferric hydroxide precipitated and removed iron from the system, the dissolution of calcite caused the net cumulative change in film thickness to initially increase by a small measure (1% between 0 and 2 pore volumes), but then monotonically decrease thereafter. Although a quantitative comparison could not be made, the comparison of experimental and simulated data demonstrated that the Mualem film model was
ARTICLE IN PRESS
Calcite Surface Area (m2/m3)
1.23E+06
4.00E+00
1.20E+06
3.00E+00
1.17E+06
2.00E+00
1.14E+06
1.00E+00
1.11E+06
Fe(OH)3 Surface Area (m2/m3)
V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
976
0.00E+00 0
2
4
8
6
Pore Volumes
(a)
Calcite
Fe(OH)3
1.0
K/Ko
0.8 0.6 0.4 0.2 0.0 0
2
4
6
8
Pore Volumes
(b)
Mualem (Lf)
C & C-G (Lf-eff)
Kozeny and Carmen
Rege & Fogler
Mualem (Lf-eff)
Fig. 7. (a) Changes in surface area for calcite and FeðOHÞ3 (am); and (b) Permeability response for ferric chloride/carbonate system for the experimental results of Rege and Fogler (1989) and the simulated results using the Kozeny and Carmen (grain packing), Mualem (film), and Childs and Collis–George (film) permeability models.
capable of predicting similar trends in permeabilities comparable to the one demonstrated by the experiment. By contrast, the lack of dependence on porosity exhibited by the Childs and Collis–George film model generated a large decrease in the permeability.
7. Summary and conclusions An approach based on the biofilm models of (Taylor et al., 1990) has been developed for modeling changes in permeability of saturated porous media due to inorganic
chemical reactions. Results of the film model implementation in the reactive transport code, STORM, demonstrate that the Mualem film permeability model can be successfully applied to mineral precipitation and dissolution reactions as demonstrated by the comparison of experimental and simulated data. The Mualem model is preferred over the model of Childs and Collis–George, because it directly relates porosity to estimates of permeability. The Mualem model, utilizing the Lf method for determining changes in pore radii, was successful in capturing the oscillatory behavior of permeability
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
observed experimentally. However, similar to the Childs and Collis–George model, the Mualem model predicted a monotonic decrease in permeability when using the Lðeff Þ method. This demonstrated that the manner in which film thickness is calculated is critical to the successful implementation of the film model within a reactive transport code. This result emphasizes an important difference in the work presented in this manuscript from the work of Taylor et al. (1990). In Taylor et al. (1990), the authors linked film thickness to permeability change as an indirect measure of the biofilm thickness. In this work, mineral film thickness was calculated from precipitation and dissolution reactions and linked to changes in permeability. Other methods of translating the mineral volume changes to estimates of film thickness may improve the film permeability model performance. Other performance enhancements may also include eliminating the assumption of a uniform film thickness, or even a different conceptual model for the shape of the pore space.
References Bacon, D.H., White, M.D., McGrail, B.P., 2004. Subsurface Transport Over Reactive Multiphases (STORM): a general, coupled, nonisothermal multiphase flow, reactive transport, and porous medium alteration simulator, version 3, user’s guide. Technical report PNNL-14783. Pacific Northwest National Laboratory, Richland, WA, p. 218. Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York, p. 764. Carnahan, C.L., 1990. Coupling of precipitation–dissolution reactions to mass diffusion via porosity changes. In: Melchior, D.C., Bassett, R.L. (Eds.), Chemical Modeling of Aqueous Systems II. American Chemical Society, Washington, DC, pp. 235–242 (Chapter 18). Childs, E.C., Collis–George, N., 1950. The permeability of porous materials. Proceedings of the Royal Society London A 201, 392–405. Cunningham, A.B., Characklis, W.G., Abedeen, F., Crawford, D., 1991. Influence of biofilm accumulation on porous media hydrodynamics. Environmental Science and Technology 25 (7), 1305–1311. Duckworth, O.W., Martin, S.T., 2003. Connections between surface complexation and geometric models of mineral
977
dissolution investigated for rhodochrosite. Geochimica et Cosmochimica Acta 67 (10), 1787–1801. Helgeson, H.C., 1969. Thermodynamics of hydrothermal systems at elevated temperatures and pressures. American Journal of Science 267, 729–804. LeGallo, Y., Bilstein, O., Brosse, E., 1998. Coupled reactionflow modeling of diagenetic changes in reservoir permeability, porosity and mineral composit. Journal of Hydrology 209 (1–4), 366–388. McGrail, B.P., Icenhower, J.P., Bacon, D.H., Vienna, J.D., Jiricka, A., Ebert, W.L., Martin, P.F., Schaef, H.T., O’Hara, M.J., Rodriguez, E.A., 1999. Waste form release data package for the 2001 immobilized low-activity waste performance assessment. Technical Report PNNL-13043, Rev. 1. Pacific Northwest National Laboratory, Richland, WA, p. 178. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research 12 (3), 513–522. Ortoleva, P., Chadam, J., Merino, E., Sen, A., 1987. Geochemical self-organization II: the reactive-infiltration instability. American Journal of Science 287 (10), 1008–1040. Pape, H., Clauser, C., Iffland, J., 2000. Variation of permeability with porosity in sandstone diagenesis interpreted with a fractal pore space model. Pure and Applied Geophysics 157 (4), 603–619. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing, Washington, DC, p. 197 Rege, S.D., Fogler, H.S., 1989. Competition among flow, dissolution, and precipitation in porous media. American Institute of Chemical Engineers Journal 35 (7), 1177–1185. Sanford, W.E., Konikow, L.F., 1989. Simulation of calcite dissolution and porosity changes in saltwater mixing zones in coastal aquifers. Water Resources Research 25 (4), 655–667. Steefel, C.I., Lasaga, A.C., 1994. A coupled model for transport of multiple chemical species and kinetic precipitation/ dissolution reactions with application to reactive flow in single phase hydrothermal systems. American Journal of Science 294 (5), 529–592. Taylor, S.W., Milly, P.D., Jaffe, P.R., 1990. Biofilm growth and the related changes in the physical properties of a porous medium 2. Permeability. Water Resources Research 26 (9), 2161–2169. Yamamoto, S., Sugiyama, S., Matsuoka, O., Kohmura, K., Honda, T., Banno, Y., Nozoye, H., 1996. Dissolution of zeolite in acidic and alkaline aqueous solutions as revealed by AFM imaging. Journal of Physical Chemistry 100 (47), 18474–18482.
Computers & Geosciences 31 (2005) 968–977 www.elsevier.com/locate/cageo
Implementation of biofilm permeability models for mineral reactions in saturated porous media$ Vicky L. Freedman, K. Prasad Saripalli, Diana H. Bacon, Philip D. Meyer Pacific Northwest National Laboratory, P.O. Box 999, MSIN K9-36, Richland, Washington 99352, USA Received 8 July 2004; received in revised form 14 February 2005; accepted 22 February 2005
Abstract An approach based on continuous biofilm models is proposed for modeling permeability changes due to mineral precipitation and dissolution in saturated porous media. In contrast to the biofilm approach, implementation of the film depositional models within a reactive transport code requires a time-dependent calculation of the mineral films in the pore space. Two different methods for this calculation are investigated. The first method assumes a direct relationship between changes in mineral radii (i.e., surface area) and changes in the pore space. In the second method, an effective change in pore radii is calculated based on the relationship between permeability and grain size. Porous media permeability is determined by coupling the film permeability models (Mualem and Childs and Collis–George) to a volumetric model that incorporates both mineral density and reactive surface area. Results from single mineral dissolution and single mineral precipitation simulations provide reasonable estimates of permeability, though they predict smaller permeability changes relative to the Kozeny and Carmen model. However, a comparison of experimental and simulated data show that the Mualem film model is the only one that can replicate the oscillations in permeability that occur as a result of simultaneous dissolution and precipitation reactions occurring within the porous media. r 2005 Published by Elsevier Ltd. Keywords: Permeability; Geochemical models; Hydrogeology; Numerical models
1. Introduction There is an increasing interest in reactive transport models that integrate porosity and permeability changes that occur as a result of mineral precipitation and dissolution in a porous medium. Such models are needed for accurately predicting contaminant fate and transport through porous media undergoing hydraulic property $ Code on server at http://www.iamg.org/CGEditor/index.htm Corresponding author. Tel.: +1 509 372 4067; fax: +1 509 372 6089. E-mail address: [email protected] (V.L. Freedman).
0098-3004/$ - see front matter r 2005 Published by Elsevier Ltd. doi:10.1016/j.cageo.2005.02.003
changes. Several models exist in the literature. In the reactive transport literature, for example, a discrete approach is usually taken to describe changes in permeability. Typically, the surface area of a mineral is used as input to the model, or calculated from the given volume fraction and initial sphere radius for the mineral (Ortoleva et al., 1987; Sanford and Konikow, 1989; Carnahan, 1990; Steefel and Lasaga, 1994; LeGallo et al., 1998). When mineral-fluid reactions occur, either the mineral radius or its surface area is recalculated accordingly. In the biological reactive solute transport literature, a continuous biofilm approach has been developed for describing hydraulic property changes due to biomass
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
accumulation in porous media. In this method, it is assumed that the solid particles that make up the porous media are uniformly covered with a biofilm. This approach assumes that the porosity and permeability reduction is thereby uniform within the pore space (Taylor et al., 1990; Cunningham et al., 1991). This uniform change in pore space may also occur for mineral reactions. For example, uniform dissolution and precipitation is consistent with surface complexation models that formulate reaction rates based on the surface concentration of specific chemical groups (Yamamoto et al., 1996; Duckworth and Martin, 2003). Observations of glass corrosion experiments for low-level vitrified waste have also shown that as the glass dissolves over time, the secondary mineral precipitates form as a uniform coating over the glass grains (McGrail et al., 1999). The goal of this work is to develop a film depositional model applicable to mineral precipitation and dissolution reactions. The model is developed to introduce an alternative permeability model that can be coupled to a reactive transport simulator for cases where uniform precipitation and dissolution may occur. Although Taylor et al. (1990) also presented a film-based porosity model, early testing demonstrated that unless the film model accounted for mineral density, it could not be adapted for use in a reactive transport model. Hence, this paper focuses on the permeability model development, and describes the film model theory and the main program functions. The model is implemented utilizing two separate methods for calculating the relative change in pore radii. An example problem involving simultaneous dissolution and precipitation is demonstrated. Model results are compared to laboratory experiments to ensure that some degree of generality can be attached to model conclusions and interpretations.
2. Film model theory The cut-and-random-rejoin approach for describing permeability assumes that porous media contains pores of different radii that are randomly distributed in space. The hydraulic conductivity of porous media depends on the number of interconnected pores and their geometric configuration. Permeability can be expressed in terms of the pore-size distribution function, which Taylor et al. (1990) define as a power function given as b r l1 f ðrÞ ¼ , (1) rmax rmax where f ðrÞ ¼ the volume of full pores of radii between r and r þ dr per unit volume of medium, r ¼ the pore radius ½L, rmax ¼ maximum pore radius ½L within the pore-size distribution ½L, and b ¼ a dimensionless
969
constant. The parameter l represents an index of the pore-size distribution, and whose range is from 0olp1, with larger values being associated with more homogeneous pore sizes. In the section that follows, Eq. (1) is used in the development of the permeability models, which closely follows the presentation given in Taylor et al. (1990). As outlined below, the development of a mineral film model differs from the biofilm model of Taylor et al. (1990) inasmuch as it accounts for a increase in the radii defined by the pore-size distribution.
3. Film model development The porosity ðnÞ of porous media can be expressed as a function of the pore-size distribution ðf ðrÞÞ and is given as Z rmax n¼ f ðrÞoðrÞ dr, (2) rmin
where r ¼ the pore radius ½L and o ¼ the weighting factor for radius r, which is assumed to be unity. Maximum and minimum pore radii are defined by rmax ½L and rmin ½L, respectively. Thus, utilizing the pore-size distribution function defined in (1), the initial diffusive porosity of a clean porous medium ðnÞ is " #1 b rmin l n¼ 1 . (3) l rmax In integral form, the Mualem (1976) and Childs and Collis–George (1950) permeabilities (km and kccg , respectively) are Z rmax 2 tn0:5 km ¼ rf ðrÞ dr , (4) 8 rmin kccg ¼
1 4
Z
rmax
r2 ½n yðrÞf ðrÞ dr,
(5)
rmin
where t ¼ tortuosity ½L, and yðrÞ ¼ the unsaturated moisture content. Substituting (1) into both permeability formulations yields ( " #)2 tb2 r2max n0:5 1 rmin 1þl km ¼ , (6) 1 8 rmax 1þl
kccg
" ( # b2 r2max 1 rmin 2þl 1 ¼ 1 4ln rmax 2þl 2 þ 2l " #) rmin 2þ2l . 1 rmax
ð7Þ
Eqs. (3), (6), and (7) describe a porous medium in the absence of a mineral film. When mineral precipitation occurs, it is assumed that a film of uniform thickness, Lf
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
970
½L, uniformly coats the internal surfaces in the porous media. Mineral film growth completely seals pores having a radius less than Lf , and reduces the radius of those pores larger than Lf . For example, pores with a new radius r, previously had a radius r þ Lf , and their cylindrical volume (assuming long tubular pores) is changed by the ratio of their volumes ½ðphÞr2 =½ðphÞðr þ Lf Þ2 , where h represents the length of the pore tube. If f p ðrÞ describes the pore-size distribution after mineral precipitation and dissolution occurs, then f ðrÞ and f p ðrÞ are related by f p ðrÞ ¼
r2 f ðr þ Lf Þ, ðr þ Lf Þ2
rmin Lf orormax Lf .
ð8Þ
By substituting f p ðrÞ for f ðrÞ in Eq. (2), the porosity resulting from mineral precipitation is Z rmax Lf r2 np ¼ f ðr þ Lf Þ dr, (9) ðr þ Lf Þ2 rmin where rmin ¼ max½rmin Lf ; 0. The permeabilities for mineral precipitation are obtained in the same manner. Substituting f p ðrÞ for f ðrÞ in (4), and (5) yields "Z #2 rmax Lf tn0:5 r3 p f ðr þ L Þ dr , (10) kmp ¼ f 8 ðr þ Lf Þ2 rmin
kccgp ¼
1 4
Z
rmax Lf
½np yr ðrÞ rmin
r4 f ðr þ Lf Þ dr. ðr þ Lf Þ2 (11)
For mineral dissolution, it is assumed that the internal surfaces of the porous media decrease by a uniform thickness. Thus, pores with a new radius r, previously had a radius r Lf , and their volume is increased by the ratio of volumes r2 =ðr Lf Þ2 . If f d ðrÞ describes the voids affected by a decrease in a uniform mineral film thickness, then f d ðrÞ and f ðrÞ are related by f d ðrÞ ¼
r2 f ðr Lf Þ, ðr Lf Þ2
rmin þ Lf orormax þ Lf .
ð12Þ
By substituting f d ðrÞ for f ðrÞ in Eq. (2), the permeabilities resulting from mineral dissolution are "Z #2 rmax þLf tn0:5 r3 d kmd ¼ f ðr Lf Þ dr , (13) 8 ðr Lf Þ2 rmin þLf
kccgd ¼
1 4
Z
rmax þLf
½nd yr ðrÞ rmin þLf
r4 f ðr Lf Þ dr. ðr Lf Þ2 (14)
The film-based permeability equations for Mualem [Eqs. (10) and (13)], and Childs and Collis–George [Eqs. (11) and (14)] are used for calculating permeability in the reactive transport simulator.
4. Program implementation The mineral film permeability model presented in this paper was first developed as a stand-alone program. Required input included porosity (n), tortuosity (t), initial pore radius (r), as well as the minimum and maximum radii of the pore-size distribution (rmin and rmax , respectively). In this work, it is assumed that tortuosity is constant. Also required as input to the program were the incremental changes in pore radii (Lf ). A critical component in the implementation of the film models in a reactive transport simulator was in determining a time-dependent method for the incremental changes in the pore space. For the stand-alone program, and in Taylor et al. (1990), it was assumed that the film thickness was static. In fact, the stand-alone program is capable of calculating both porosity and permeability based on film models that utilize a static film thickness. However, the stand-alone model does not consider mineral surface areas and densities, which is an important aspect for implementation with a reactive transport simulator. Hence, porous media porosity was determined using a volumetric model that considered both density and surface area. For permeability calculations, the mineral film thickness needed to dynamically change as a function of time-dependent reaction processes. To this end, mineral volume changes due to mineral precipitation were interpreted as a change in the film thickness that uniformly altered the radii of all pores in the pore-size distribution. For dissolution, the change in volume also occurred uniformly on the solid-phase surface. Thus, a negative change in the mineral film thickness was translated into an increase in the pore radii throughout the entire pore-size distribution. To calculate changes in hydraulic properties within the reactive transport simulator, porous media permeability was determined by coupling the film permeability models to a volumetric model that incorporated both mineral density and reactive surface area. Once new mineral volumes were determined, the film permeability models required the change in pore radii (Lf ) to be dynamically calculated from the precipitation and dissolution reactions. Two separate methods for determining the change in pore radii were used. In the first method, it was assumed that the cumulative pore radii changes were equal to the cumulative changes in mineral radii: X X t1 Ltf ¼ Rtm Rt1 (15) m þ Lf ,
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
where t ¼ the current timestep, t 1 ¼ the previous timestep, and Rm ¼ mineral radius ½L, which is not related to the pore radius (r) initially described in Eq. (1). Alternatively, an effective Lf was determined based on a relationship between permeability and grain size. Like the film models of Taylor et al. (1990), the Kozeny and Carmen permeability model is based on conduit flow, which assumes steady, laminar flow through bundles of capillary tubes. Because it is difficult in natural systems to determine the tube radius, the Kozeny and Carmen model was modified to relate permeability to the mineral grain size. An earlier version of the Kozeny and Carmen equation [in Pape et al., 2000] attempts to incorporate the effects of effective hydraulic pore radius and the length of parallel tubes k¼
nr2eff
. (16) 8t Although the effective hydraulic radius cannot be directly measured, it can be determined from the following relationship (Pape et al., 2000) reff ¼
2Rm tn . 3ð1 nÞ
(17)
By including an expression for the effective porosity, Eq. (17) accounts for the fact that only a fraction of the pores are interconnected and participate in the flow process. Thus, in the second method, the cumulative Lf ðeff Þ was determined as 2Rm tn t 2Rm tn t1 t Lf ðeff Þ ¼ þ Lft1 (18) ðeff Þ . 3ð1 nÞ 3ð1 nÞ Eq. (18) is also consistent with the Mualem permeability model in as much as it accounts for the tortuosity of the flow path. In Taylor et al. (1990), analytical solutions were developed to calculate permeability. In this work, numerical solutions for the integrals in Eqs. (10), (11), (13) and (14) were generated using the mathematical software Maple (by Maplesoft).
5. Reactive transport simulator The film permeability models were coupled to the reactive mass transport and flow simulator, Subsurface Transport Over Reactive Multiple phases (STORM) (Bacon et al., 2004). In STORM, flow of fluid phases was computed from Darcy’s law. The resulting flow fields were then used to solve conservation equations for reactive aqueous phase transport. Transport of phase components was computed from Fick’s law. Once flow and transport calculations were complete, porosity and permeability calculations were performed. For brevity’s sake, the governing equations describing flow, transport
971
and chemical reactions are not presented in this document, but are described in detail in Bacon et al. (2004). However, a brief description of the porosity and permeability models used in STORM is provided below. The two governing equations describing saturated flow and advective-dispersive transport are discretized into algebraic equations using the control volume method (Patankar, 1980). The nonlinear equation sets are converted into linear form utilizing iterative Newton–Raphson method. Details of the numerical solution methods can be found in Bacon et al. (2004). A primary assumption associated with the solute mass conservation equation is that solute concentrations are infinitely dilute. This assumption implies that fluid properties are independent of solute concentrations, which allows the solute conservation equations to be solved independently of the flow equation. Thus, once the flow equation is solved, the chemical reaction equations and the solute transport equations are solved simultaneously. Reactions that only involve aqueous species are assumed to be at equilibrium. Mineral precipitation and dissolution reactions are treated kinetically. Reaction rates in STORM are based on transition state theory and consider the forward rate and the equilibrium constant 2 3 , jnij j jnij j Y Y W j ¼ AmðjÞ kmðjÞ 4 ðgi ci Þ f j ðgi ci Þ f j ðK rðjÞ Þ1=f 5, nij o0
nij 40
(19) where W j ¼ the reaction rate for reaction j [mol/MT], Am ðjÞ ¼ the effective mineral reaction area ½L2 ], kj ¼ the rate constant [1/t], K r ðjÞ ¼ the equilibrium constant, gi ¼ the activity coefficient of species i, ci ¼ the concentration [moles/M], and f i ¼ the factor of the reaction order. The stoichiometric coefficient of species i in reaction j is represented by nij . The bars j j represent the absolute value and the sign of the reaction is positive for precipitation and negative for dissolution. Once the flow, transport and chemical reaction equations are solved, porosity and permeability updates are then calculated by the STORM texture solver, and are used as input values for flow in the time stepping routine. 5.1. Volumetric equations STORM calculates volume fractions by assuming that mineral grains are spherical in shape, and occupy a given volume based on the number of grains and the mineral radius. Assuming that the number of grains remains constant, changes in porosity are calculated based on a change in the mineral radius. STORM describes the relationship between the mineral volume fraction ðV m ½L3 =L3 Þ and the porosity
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
972
Table 1 Parameter values for numerical simulations
of the porous medium as n ¼ 1:0
Nm X
V m.
(20)
Property
Value
Initial porosity (n) Initial permeability (k) Maximum pore radius (rmax )
0.50 2:92 1010 m2 1:125 103 m (Mualem) 1:441 104 m (Childs and Collis–George) 0.0 m 3.0 1:86 102 m 5:84 1010 m
m¼1
where N m is the total number of minerals. Assuming a spherical geometry, the number of grains (nm ) present for each mineral is defined as nm ¼
3V m 4pR3m
(21)
and the effective reactive surface area (Am ) is Am ¼ 4pR2m nm .
(22)
By assuming that a dissolution/precipitation reaction is a surface-limited kinetic reaction, the change rate of the volume fraction is qV m qRm ¼ 4pR2m nm . qt qt
(23)
Minimum pore radius (rmin ) Pore-size distribution (l) Tortuosity (t) Kozeny and Carmen fitting coefficient ðC o Þ Calcite initial surface area (SA) Quartz initial surface area (SA) Darcy flux ðqÞ Domain length ðLÞ Grid discretization ðDzÞ Total simulation time (T)
1:2084 102 m2 =m3 2:5000 105 m2 =m3 2:5 105 m=s 2.0 m 0.1 m 5.0 days
5.2. Permeability equations Permeability in the STORM simulator is calculated by a relationship defined by Kozeny and Carmen. This permeability equation is empirical and requires the input of fitting parameters. For example, the Kozeny and Carmen ðk ½L2 Þ equation is given as (Bear, 1972) k ¼ Co
n3 , ð1 nÞ2 A2s
(24)
where As is the specific surface area and C o ¼ an empirical fitting coefficient.
6. Model testing Evaluation of the film depositional models was performed by comparing general trends in hydraulic property changes with the Kozeny and Carmen model [Eq. (24)] already implemented in STORM. By adjusting the fitting parameters, the initial permeabilities were approximately equal between the STORM and film models at the beginning of the simulations. To calculate changes in permeability, the models required n, l, rmin , and rmax as input. Because hydraulic conductance largely depends on the connectivity of larger pore sizes, rmax essentially determined the magnitude of permeability. Since the permeability was insensitive to the lower limit of integration, the minimum radius was assumed to have a value close to zero rmin ¼ 1 1015 . These calculations were independent of the initial mineral radius [Eq. (22)] required for input in STORM.
Although several example problems were run to test the depositional models, results from only a few representative simulations are presented in the sections that follow.
6.1. Precipitation A one-dimensional, saturated flow simulation involving the precipitation of calcite ðCaCO3 Þ at 25 1C was used to test the film models. Initially, the column was slightly saturated with respect to calcite. Input parameters appear in Table 1. For flow, a specified flux boundary condition was assigned at the top of the column, and a free gradient boundary at the bottom. An initial conditions boundary type was assigned at both the top and bottom boundaries for solute concentrations. This boundary type fixed the concentration to the initial concentration value at the node adjacent to the boundary surface, and was invariant with time (Bacon et al., 2004). At time greater than zero, a water flux entered the column from the top, which maintained super saturation with respect to calcite. As the calcite precipitated, incremental changes in the pore radii occurred, causing a porosity and permeability reduction. The precipitation of calcite in the column was simulated for 5 days. Results from this simulation are presented below. In Figs. 1 and 2, results for surface area changes and findings from the Mualem film permeability model are plotted with STORM’s implementation of the Kozeny
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
973
Distance from Outlet (m)
2.0 Kozeny and Carmen C and C-G (Lf) C and C-G (Lf-eef)
1.5
1.0
0.5
0.0 0.6
0.8
1.0
1.2
K/Ko
Fig. 1. Relative changes in reactive surface area for calcite precipitation simulation.
Fig. 3. Permeability reduction profiles for the Kozeny and Carmen (grain packing) and Childs and Collis–George film models.
Distance from Outlet (m)
2.0 Kozeny and Carmen Mualem (Lf) Mualem (Lf-eef)
1.6 1.2 0.8 0.4 0.0 0.6
0.8
1.0
1.2
K/Ko
Fig. 2. Permeability reduction profiles for the Kozeny and Carmen (grain packing) and Mualem film models.
and Carmen model, respectively. Results of the Mualem permeability calculations showed that even though surface area changes were the same for all three permeability models, the magnitude of the Kozeny and Carmen permeability reduction is greater than both of the implementations of the Mualem film model. At the top of the column, the relative change in surface area was only 13%, yet the permeability reduction was 31% for the Kozeny and Carmen model, and only 26% and 11% for the Mualem models ðLf and Lf ðeff Þ , respectively). The shape of the permeability reduction profiles was fairly similar for the Kozeny and Carmen and Mualem ðLf Þ model, but the effective ðLf ðeff Þ Þ implementation was less steep as the permeability reduction predicted over the length of the column was less than the other two models. Shown in Fig. 3 are predictions from the implementation of the Childs and Collis–George film permeability model with the Kozeny and Carmen equation. Similar to the Mualem permeability model, the Childs and Collis–George equation experienced a larger permeabil-
ity reduction due to calcite precipitation when Lf was calculated as in Eq. (15). However, this difference was quite small. At the top of the column, the Lf implementation of the Childs and Collis–George model predicted a 5% permeability reduction, whereas the Lf ðeff Þ implementation predicted a 4% reduction. Relative to the Kozeny and Carmen model, the Childs and Collis–George film models estimate a smaller permeability reduction. As is evident from the results shown in Figs. 2 and 3, predictions made by the Kozeny and Carmen model showed a more dramatic decrease in permeability throughout the column relative to all of the film model implementations. Similar behavior was demonstrated in the permeability changes for mineral dissolution presented in the following section. 6.2. Dissolution Inorganic reactions may cause an increase in the porosity and permeability of the porous media by increasing the pore radii due to solid phase dissolution. To test the film growth models in STORM, quartz dissolution simulations for a one-dimensional column were conducted at 100 1C. The column was undersaturated with respect to quartz, but otherwise had the same physical setup as the calcite precipitation test case (Table 1). Shown in Fig. 4 is the relative changes in surface area profile, and in Figs. 5 and 6 the permeability profiles for the Kozeny and Carmen model as well as the Mualem and Childs and Collis–George film models, respectively. Fig. 5 shows that the two lines representing the permeability enhancement profiles for each implementation of the Mualem model coincide, as the predictions were nearly identical. Although the maximum change in surface area was only 5% at the top of the column, both implementations predict only a 4% increase in the
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
974
Fig. 4. Relative changes in reactive surface area for quartz dissolution simulation.
6.3. Column experiment with simultaneous precipitation and dissolution
Distance from Outlet (m)
2.0 1.6 Kozeny and Carmen Mualem (Lf) Mualem (Lf-eef)
1.2 0.8 0.4 0.0 0.8
1.0
1.2
1.4
1.6
K/Ko
Fig. 5. Permeability enhancement profiles for the Kozeny and Carmen (grain packing) and Mualem film models.
Distance from Outlet (m)
2.0 1.6 1.2 Kozeny and Carmen C and C-G (Lf)
0.8
C and C-G (Lf-eff) 0.4 0.0 0.8
two lines representing the Childs and Collis–George permeability profiles coincide in Fig. 6, but relative to the Kozeny and Carmen model, the two implementations predict a smaller permeability enhancement (7%). Although all of the permeability profiles presented for both mineral precipitation and dissolution reactions yield plausible results, the film models estimated smaller permeability changes to the Kozeny and Carmen model. However, comparisons of film-based permeability estimates to the Kozeny and Carmen relation were only used for reference, not as an absolute scale of performance. To this end, permeability predictions for a column experiment are examined in the following section.
1.0
1.2 K/Ko
1.4
1.6
Fig. 6. Permeability enhancement profiles for the Kozeny and Carmen (grain packing) and Childs and Collis–George film models.
permeability at the top of the column compared to a 15% increase predicted by Kozeny and Carmen. A similar result occurred for the predictions made by the Childs and Collis–George models. Not only do the
Matrix acidizing is a process which uses acids to dissolve minerals near a wellbore to increase permeability and enhance oil recovery. However, some of the dissolution products can precipitate as the acid reacts with geologic materials. Commonly, iron hydroxides are observed to plug pore throats and decrease permeability. The simulations presented in this section are based on the experiments of Rege and Fogler (1989). In their investigation, Rege and Fogler (1989) acidized a limestone core using 0.75 M iron chloride with a pH of one. Dissolution of the calcite ðCaCO3 ) by the acid caused the solution pH to increase and ferric hydroxide ðFeðOHÞ3 ) to precipitate. Variations in the permeability occurred as a result of the simultaneous precipitation and dissolution reactions occurring within the core. Both the film models and the Kozeny and Carmen equation are compared to the experimental permeability ratio reported by Rege and Fogler (1989). The onedimensional column measures 0.089 m, and is discretized into 20 equally spaced nodes. Since the initial porosity and permeability of the experimental limestone core are not reported, the porosity was assumed to be equivalent to the value used in the previous test simulations reported in Table 1. The initial permeability of the medium was assumed to be 5:8177 1010 m2 . The kinetic reaction rate for FeðOHÞ3 ðsÞ was an order of magnitude lower ð1 107 mol=m2 s) than that of the CaCO3 ð1 106 mol=m2 s). The flow rate for the simulations is 1 108 m3 =s, and is equivalent to the one used in the acidizing matrix experiments of Rege and Fogler (1989). A reaction front in the simulations is created by injecting a 0.075 M solution of ferric chloride ðFeCl3 Þ with a pH of 1.7 into a one-dimensional column containing calcite. The solution in the column is initially undersaturated with respect to calcite and has a pH of 7.9. This setup differs from the pH and strength of concentration (1.0 and 0.75 M, respectively) of the
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
injected solution in the experiments. However, concentration limits imposed by both the gauss-seidel iterative solver in STORM, and the b-dot ion activity model (Helgeson, 1969), which is valid for ionic strengths p0:1 M, prevented an exact replication of the Rege and Fogler (1989) experiments. Hence, the simulations presented here only attempt to replicate the permeability trends reported by Rege and Fogler (1989), and not produce an exact replication of the chemistry of their system. 6.3.1. Comparison of results The changes in reactive surface area for both calcite and FeðOHÞ3 ðamÞ for the top of the column are found in Fig. 7a, and the calculated permeability responses for this work, as well as the one found by Rege and Fogler (1989) are shown in Fig. 7b. This includes the simulated permeability ratios for the Kozeny and Carmen, Mualem and Childs and Collis–George models at 0.0134 m from the top of the column. In all cases, acid is consumed through the dissolution of calcite, which is shown to decrease its surface area throughout the simulation. As the pH increases, the iron in solution precipitates as iron hydroxide. Because the rate of iron hydroxide precipitation exceeds the rate of calcite dissolution, the permeability initially decreased in the column. For the Mualem predictions, both the Lf and Lf ðeff Þ methods were used to calculate changes in the pore space. Although an attempt was made to use the Lf method for the Childs and Collis–George model, the permeability decreased to zero and convergence was not obtained. Hence, for the Childs and Collis–George permeability predictions shown in Fig. 7, only the Lf ðeff Þ method is shown. In the Rege and Fogler (1989) experiments, fluctuations in the permeability response are observed that are absent in all but the Mualem Lf simulation. Rege and Fogler (1989) explain that these oscillations were due to the deposition of the precipitate ahead of the dissolution front. Initially, the precipitate formed in flow paths with higher flow rates. As deposition of the precipitate occurred, there was an increase in the flow resistance, which caused flow to be diverted to flow paths with lower rates of flow. As the smaller flow paths lengthened, flow resistance decreased. This caused the sudden increase in permeability at 8 pore volumes. In addition to the diverted flow paths, the oscillations and excursions in the permeability profiles of the (Rege and Fogler, 1989) experiments and in the Mualem permeability profile, could also be a result of the cyclical nature of precipitation–dissolution reaction events in the medium. The rate equation for a typical precipitation–dissolution reaction shows that the rate of reaction is proportional to ð1 Q=KÞ, where Q is the activity product of aqueous species involved in reaction, and K is
975
the equilibrium constant for reaction. Although precipitate formation is thermodynamically favorable at or above saturation values (i.e., when Q=K41Þ, precipitation will not generally occur until a significant degree of super-saturation is reached. Deposition of the precipitated mass in the porous medium is a faster, or more sudden event compared to the slow dissolution process. Once such precipitation occurs, the pore fluids are depleted of the solute mass and the system becomes thermodynamically more favorable for dissolution, which starts a fresh precipitation–dissolution cycle. Since the permeability of the porous medium is sensitive to the thickness of the precipitate film [as shown in Eqs. (4) and (5)], it follows that the cyclic or oscillatory precipitation–dissolution behavior would likely lead to a similar oscillatory change in permeability as seen in Fig. 7. Using the Lf method, the Mualem model was able to capture this behavior, due to the sensitivity of the method to mineral volume changes. Not only was the Mualem model able to replicate the oscillations in permeability, but it also predicted a general increase in permeability at 7 pore volumes, similar to the increase in permeability observed at 8 pore volumes. Because of the differences in the concentration and pH of the injected fluid, the simulated permeability decrease occurred less rapidly than the one observed experimentally. With the Lf ðeff Þ method, however, the translation of the mineral volume changes to changes in pore radii were too small to capture this behavior. In fact, a monotonic permeability reduction was predicted by both the Mualem ðLf ðeff Þ ) and Childs and Collis–George ðLf ðeff Þ ) models because of the insensitivity of the method to the increases in the pore radii. A less dramatic, monotonic permeability reduction was predicted by the Kozeny and Carmen model, but the Kozeny and Carmen model also failed to capture the oscillatory behavior of the permeability in the system. In contrast to the Mualem and Kozeny and Carmen models, the Childs and Collis–George film model predicted dramatic changes in permeability. Between 4 and 9 pore volumes, the permeability decreased by nearly 100%. This decrease in permeability was not due to the lengthening of smaller flow paths, rather fluctuations in film thicknesses. Because the Childs and Collis–George film model does not include an explicit relationship to the porosity, changes in permeability are entirely dependent on film thickness. As the ferric hydroxide precipitated and removed iron from the system, the dissolution of calcite caused the net cumulative change in film thickness to initially increase by a small measure (1% between 0 and 2 pore volumes), but then monotonically decrease thereafter. Although a quantitative comparison could not be made, the comparison of experimental and simulated data demonstrated that the Mualem film model was
ARTICLE IN PRESS
Calcite Surface Area (m2/m3)
1.23E+06
4.00E+00
1.20E+06
3.00E+00
1.17E+06
2.00E+00
1.14E+06
1.00E+00
1.11E+06
Fe(OH)3 Surface Area (m2/m3)
V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
976
0.00E+00 0
2
4
8
6
Pore Volumes
(a)
Calcite
Fe(OH)3
1.0
K/Ko
0.8 0.6 0.4 0.2 0.0 0
2
4
6
8
Pore Volumes
(b)
Mualem (Lf)
C & C-G (Lf-eff)
Kozeny and Carmen
Rege & Fogler
Mualem (Lf-eff)
Fig. 7. (a) Changes in surface area for calcite and FeðOHÞ3 (am); and (b) Permeability response for ferric chloride/carbonate system for the experimental results of Rege and Fogler (1989) and the simulated results using the Kozeny and Carmen (grain packing), Mualem (film), and Childs and Collis–George (film) permeability models.
capable of predicting similar trends in permeabilities comparable to the one demonstrated by the experiment. By contrast, the lack of dependence on porosity exhibited by the Childs and Collis–George film model generated a large decrease in the permeability.
7. Summary and conclusions An approach based on the biofilm models of (Taylor et al., 1990) has been developed for modeling changes in permeability of saturated porous media due to inorganic
chemical reactions. Results of the film model implementation in the reactive transport code, STORM, demonstrate that the Mualem film permeability model can be successfully applied to mineral precipitation and dissolution reactions as demonstrated by the comparison of experimental and simulated data. The Mualem model is preferred over the model of Childs and Collis–George, because it directly relates porosity to estimates of permeability. The Mualem model, utilizing the Lf method for determining changes in pore radii, was successful in capturing the oscillatory behavior of permeability
ARTICLE IN PRESS V.L. Freedman et al. / Computers & Geosciences 31 (2005) 968–977
observed experimentally. However, similar to the Childs and Collis–George model, the Mualem model predicted a monotonic decrease in permeability when using the Lðeff Þ method. This demonstrated that the manner in which film thickness is calculated is critical to the successful implementation of the film model within a reactive transport code. This result emphasizes an important difference in the work presented in this manuscript from the work of Taylor et al. (1990). In Taylor et al. (1990), the authors linked film thickness to permeability change as an indirect measure of the biofilm thickness. In this work, mineral film thickness was calculated from precipitation and dissolution reactions and linked to changes in permeability. Other methods of translating the mineral volume changes to estimates of film thickness may improve the film permeability model performance. Other performance enhancements may also include eliminating the assumption of a uniform film thickness, or even a different conceptual model for the shape of the pore space.
References Bacon, D.H., White, M.D., McGrail, B.P., 2004. Subsurface Transport Over Reactive Multiphases (STORM): a general, coupled, nonisothermal multiphase flow, reactive transport, and porous medium alteration simulator, version 3, user’s guide. Technical report PNNL-14783. Pacific Northwest National Laboratory, Richland, WA, p. 218. Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York, p. 764. Carnahan, C.L., 1990. Coupling of precipitation–dissolution reactions to mass diffusion via porosity changes. In: Melchior, D.C., Bassett, R.L. (Eds.), Chemical Modeling of Aqueous Systems II. American Chemical Society, Washington, DC, pp. 235–242 (Chapter 18). Childs, E.C., Collis–George, N., 1950. The permeability of porous materials. Proceedings of the Royal Society London A 201, 392–405. Cunningham, A.B., Characklis, W.G., Abedeen, F., Crawford, D., 1991. Influence of biofilm accumulation on porous media hydrodynamics. Environmental Science and Technology 25 (7), 1305–1311. Duckworth, O.W., Martin, S.T., 2003. Connections between surface complexation and geometric models of mineral
977
dissolution investigated for rhodochrosite. Geochimica et Cosmochimica Acta 67 (10), 1787–1801. Helgeson, H.C., 1969. Thermodynamics of hydrothermal systems at elevated temperatures and pressures. American Journal of Science 267, 729–804. LeGallo, Y., Bilstein, O., Brosse, E., 1998. Coupled reactionflow modeling of diagenetic changes in reservoir permeability, porosity and mineral composit. Journal of Hydrology 209 (1–4), 366–388. McGrail, B.P., Icenhower, J.P., Bacon, D.H., Vienna, J.D., Jiricka, A., Ebert, W.L., Martin, P.F., Schaef, H.T., O’Hara, M.J., Rodriguez, E.A., 1999. Waste form release data package for the 2001 immobilized low-activity waste performance assessment. Technical Report PNNL-13043, Rev. 1. Pacific Northwest National Laboratory, Richland, WA, p. 178. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research 12 (3), 513–522. Ortoleva, P., Chadam, J., Merino, E., Sen, A., 1987. Geochemical self-organization II: the reactive-infiltration instability. American Journal of Science 287 (10), 1008–1040. Pape, H., Clauser, C., Iffland, J., 2000. Variation of permeability with porosity in sandstone diagenesis interpreted with a fractal pore space model. Pure and Applied Geophysics 157 (4), 603–619. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing, Washington, DC, p. 197 Rege, S.D., Fogler, H.S., 1989. Competition among flow, dissolution, and precipitation in porous media. American Institute of Chemical Engineers Journal 35 (7), 1177–1185. Sanford, W.E., Konikow, L.F., 1989. Simulation of calcite dissolution and porosity changes in saltwater mixing zones in coastal aquifers. Water Resources Research 25 (4), 655–667. Steefel, C.I., Lasaga, A.C., 1994. A coupled model for transport of multiple chemical species and kinetic precipitation/ dissolution reactions with application to reactive flow in single phase hydrothermal systems. American Journal of Science 294 (5), 529–592. Taylor, S.W., Milly, P.D., Jaffe, P.R., 1990. Biofilm growth and the related changes in the physical properties of a porous medium 2. Permeability. Water Resources Research 26 (9), 2161–2169. Yamamoto, S., Sugiyama, S., Matsuoka, O., Kohmura, K., Honda, T., Banno, Y., Nozoye, H., 1996. Dissolution of zeolite in acidic and alkaline aqueous solutions as revealed by AFM imaging. Journal of Physical Chemistry 100 (47), 18474–18482.