Experimental and modeling investigation of multicomponent reactive transport in porous media

Journal of Contaminant Hydrology 120–121 (2011) 27–44

Contents lists available at ScienceDirect

Journal of Contaminant Hydrology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j c o n h yd

Experimental and modeling investigation of multicomponent reactive transport in porous media Guy E. Katz a, Brian Berkowitz a,⁎, Alberto Guadagnini b, Maarten W. Saaltink c a b c

Department of Environmental Sciences and Energy Research, Weizmann Institute of Science, Rehovot 76100, Israel Dipartimento di Ingegneria Idraulica, Ambientale, Infrastrutture Viarie, Rilevamento (DIIAR), Politecnico di Milano, Piazza L. Da Vinci, 32, 20133 Milano, Italy Department of Geotechnical Engineering and Geosciences, Technical University of Catalonia, C/Jordi Girona 1-3, 08034 Barcelona, Spain

a r t i c l e

i n f o

Article history: Received 4 June 2009 Received in revised form 15 October 2009 Accepted 13 November 2009 Available online 3 December 2009 Keywords: Precipitation Calcium carbonate Advection Diffusion Dispersion

a b s t r a c t We present an experimental and modeling study of solute transport in porous media in the presence of mixing-induced precipitation of a solid phase. Conservative and reactive transport experiments were performed in a quasi-two-dimensional laboratory flow cell, filled with homogeneous and heterogeneous porous media. Conservative experiments were performed by injecting solutions containing sodium chloride and calcium chloride into the domain. In reactive transport experiments, inlet solutions of calcium chloride and sodium carbonate were injected in parallel, resulting in calcium carbonate precipitation where the solutions mix. Experimental results were used as a benchmark to examine the performance of a reactive transport numerical model. Good agreement between model predictions and experimental results was obtained for the conservative transport experiments. The reactive transport experiments featured the formation of a calcium carbonate mineral phase within the mixing zone between the two solutions, which controlled the spatial evolution of calcium carbonate in the domain. Numerical simulations performed on high resolution grids for both the homogeneous and heterogeneous porous systems underestimated clogging of the system. Although qualitative agreement between model results and experimental observations was obtained, accurate model predictions of the spatial evolution of calcium concentrations at sample points within the flow cell could not be achieved. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Transport of reactive species, like most geochemical phenomena of interest in Earth science, is characterized by a combination of physical and chemical processes. The key role of mixing in reactive transport phenomena is welldocumented (Steefel and MacQuarrie, 1996; De Simoni et al., 2005), and is recognized for its control over the chemical nature of subsurface systems through aqueous speciation and its forcing of multiphase reactions.

⁎ Corresponding author. Tel.: +972 8 9342098. E-mail address: [email protected] (B. Berkowitz). 0169-7722/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2009.11.002

Special attention has been given to the geochemistry of carbonate systems at the interface between fresh water and sea water (e.g., Sanford and Konikow, 1989; Corbella et al., 2003; Singurindy et al., 2004; Rezaei et al., 2005; Romanov and Dreybrodt, 2006), where mixing of waters that are initially in equilibrium with a solid phase drives precipitation and dissolution of minerals. Diagenesis of carbonate rocks is affected by mixing of waters in many field sites. For example, carbonate dissolution was observed in coastal aquifers in Yucatan, Mexico (Back et al., 1979), the coastal springs of Greece (Higgins, 1980), and Andores Island in the Bahamas (Smart et al., 1988). Carbonate coastal aquifer systems constitute a striking example of the intrinsic complexity of subsurface mixingdriven reactions, which are associated with different governing processes. First, solubility of constituents in the water depends

28

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

in a nonlinear fashion on temperature, pH, salinity and other factors (Singurindy and Berkowitz, 2003). Moreover, the importance of precipitation and dissolution reactions is not limited to their effect on groundwater chemistry; in addition these reactions can modify significantly the physical properties of porous media (Lasaga, 1984; Daccord et al., 1993; Dijk and Berkowitz, 1998; Singurindy and Berkowitz, 2003; Singurindy and Berkowitz, 2004; Emmanuel and Berkowitz, 2005; Romanov and Dreybrodt, 2006; Tartakovsky et al., 2008). Under some conditions, large quantities of mass can be transferred between liquid and solid mineral phases, altering the porosity and permeability of the rock formation (Domenico and Schwartz, 1998). The change in physical properties may in turn change the flow pattern, thus causing a feedback on the reaction rates through solute transport variations. Several experimental and theoretical studies have analyzed the influence of dissolution and precipitation on the chemical and physical properties of soluble porous matrices (Fogler and Rege, 1986; Hoefner and Fogler, 1988; Kelemen et al., 1995). Similarly, some understanding of the evolution of hydraulic conductivity caused by precipitation has been gained by experimental and theoretical study (Novak et al., 1989; Novak, 1993; Bolton et al., 1996, 1997; Dijk and Berkowitz, 1998; Cochepin et al., 2008; Tartakovsky et al., 2008). Quantitative analysis and prediction of reactive transport processes in natural groundwater systems is typically accomplished by modeling tools. The latter rely on a variety of mathematical formulations to obtain estimates of spatial distributions of chemical species concentrations and reaction rates when considering transport of multiple species with different types of reactions. Many codes have been developed over the last decades to address these problems (e.g., Cederberg et al., 1985; Liu and Narasimhan, 1989; Kinzelbach et al., 1991; Mangold and Tsang, 1991; Yeh and Tripathi, 1991; Steefel and Lasaga, 1994; Walter et al., 1994; Wunderly et al., 1996; Saaltink et al., 2004). The complexity of reactive transport problems presents challenges for code development, accounting for all types of chemical reactions in a concise mathematical formulation. Available codes typically formulate the reactive transport problem by using components (sometimes referred to as master species). These can be defined as linear combinations of the concentrations of chemical species. In some cases the transport equations of these components become conservative (i.e., independent of reactions) and decoupled from each other (Molins et al., 2004; De Simoni et al., 2005). Identification and quantification of reactive transport processes are usually based on controlled laboratory experiments performed in flow cells (e.g., Tartakovsky et al., 2008; Guadagnini et al., 2009). Typical cell configurations involve one- or two-dimensional domains, usually not longer or wider than tens of centimeters, filled with a desired porous material. Such experiments have several advantages with respect to field-scale aquifer analyses: (i) initial and boundary conditions of both chemical and flow related quantities can be established with relative ease and monitored throughout the experiment; (ii) the initial distribution of physical and geochemical parameters can be identified; and (iii) experiments are relatively fast. On the other hand, direct measurement of the space–time distribution of solute concentrations within the flow cell is

not easy to accomplish, although it is this type of information that is needed to assess the predictive capability of a (conceptual and numerical) model. This problem can be addressed, for instance, by employing sampling ports at selected points in the flow cell, although only limited amounts of liquid can be withdrawn without disturbing the flow field and the reaction dynamics. While ports could allow measurement of the desired concentrations of several species, there are no direct means to measure the space–time evolution of the reaction rates taking place within the cell. Reactions can be measured in only a few cases, e.g., where the products of the reactions lead to changes in color or luminescence and can be measured indirectly (Gramling et al., 2002; Cirpka et al., 2006). In practice reaction rates are usually calculated from a mass balance of reactants in the inflowing and outflowing waters. However, this can only provide an integrated measure of reaction rates and not space–time patterns of local rates. Thus, importantly, this information is insufficient to identify the position and magnitude of reactions that may significantly impact the evolution of rock hydraulic properties (Daccord et al., 1993; Dijk and Berkowitz, 1998; Singurindy and Berkowitz, 2004; Emmanuel and Berkowitz, 2005; Guadagnini et al., 2009). Additional concerns are associated with the ability to upscale the behavior observed in flow cells to field-scale systems. Modeling porosity evolution within mixing regions during reactive transport has important applications in several environmental and engineering problems. Assessing the reliability of numerical codes to this end is not an easy task, due to the lack of benchmark analytical solutions and paucity of exhaustive experimental datasets under flow through conditions. A commonly adopted strategy relies on setting a benchmark problem and comparing the results provided by different numerical models to identify their critical similarities and differences (e.g., Cochepin et al., 2008). A recent example of this is offered by the efforts coordinated in the context of the MoMaS (Modeling, Mathematics and numerical Simulations related to nuclear waste management problems) modeling initiative (http://www.gdrmomas.org/ benchmarks-en.html). Tartakovsky et al. (2008) performed a laboratory experiment and a numerical study of a mixing-induced reaction of sodium carbonate and calcium chloride solutions that featured changes in domain porosity and associated hydraulic properties of an initially homogeneous flow cell. Their experiments evidenced the formation of a narrow calcium carbonate precipitate layer between the two injected solutions. The width of the precipitate layer was approximated from visual inspection of the flow cell; concentrations of solutes were not measured in the cell or at the outlet. To quantify such processes, an interpretive model based on the commonly adopted localized format of the (otherwise nonlocal) space-averaged advection–dispersion–reaction equation (ADRE) is often considered. Tartakovsky et al. (2008) note, however, that several scale separation conditions must be maintained (e.g., the characteristic length associated with the pore space geometry must be much smaller than the characteristic length of the averaging volume) to allow use of the ADRE. These conditions are often violated in transport processes that involve mixing, which suggests that ADRE-based models are of limited applicability

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

to this type of problem. On the other hand, Tartakovsky et al. (2008) showed that solving an ADRE-based model on a high resolution grid can describe qualitatively some features of such a process. In the current study, we use a similar experimental setup and perform a set of detailed measurements within the flow cell, to enable semi-quantitative comparison between experiments and the results of a numerical model. Solutions of Na2CO3 and CaCl2 were injected in parallel into two halves of a flow cell, yielding formation of a solid CaCO3 phase in a narrow zone where the two solutions mixed. Several experiments were conducted with homogeneous and heterogeneous porous medium packing patterns. The main objectives of this research were: (i) to acquire high resolution measurements of solute concentrations in the flow domain, which can also be used as a future benchmark test for numerical models; (ii) to provide a semi-quantitative analysis of the calcite evolution patterns within the flow cell; and (iii) to test the ability of a well-documented numerical model (Retraso Code Bright; Saaltink et al., 2004) to reproduce the evolution of solute concentrations measured within the system during the precipitation experiments. 2. Theoretical basis Modeling reactive transport involves describing mass balances of species and reactions among species. The basic equations, an ADRE formulation, and the strategy adopted to solve them, are described below.

matrix, Nr and Ns being the number of reactions and of chemical species, respectively). Equilibrium is described by the mass action law, which can be written as Se log a = log K, where K is the vector of chemical equilibrium constants and a is the vector of species activities. Some species (usually called constant activity species), such as water and minerals, can be assumed to have unit activity. The mass action law can be written such that the activities of Nr secondary species can be calculated from the activities of Ns −Nr primary species (Steefel and MacQuarrie, 1996; Saaltink et al., 1998; Molins et al., 2004; De Simoni et al., 2005, 2007). Following De Simoni et al. (2007), we choose as primary species the Nc constant activity species plus Ns −Nr −Nc aqueous species. All secondary species are aqueous, as mineral species are considered to have constant activities. Vector a is split as a = (aca′aa″a)T, where ac contains the activities of the Nc constant activity species, a′a contains the Nu (=Ns − Nr −Nc) activities of aqueous primary species, and a″a is composed by the Nr activities of secondary species. Likewise, we subdivide Se into three parts. It is always possible (e.g., Saaltink et al., 1998) to redefine the chemical system so that Se = (Sec|S′ea| − I), where Sec, S′ea and −I correspond, respectively, to the stoichiometric coefficients of constant activity, primary and secondary species. This allows explicit calculation of secondary species activities from mass action laws. In term of concentrations, c = (ccca)T = (ccc′ac″a)T, one can write the mass action law as ′ logc′a −logK⁎ ; logc″a = Sea

ð2Þ

where K⁎ is a vector of equivalent equilibrium constants defined as

2.1. Species transport equations A species mass balance can be written as ∂ðmÞ = MLðcÞ + f : ∂t

29

⁎ ′ j−IÞlogγa; logK = logK−ðSea ð1Þ

Here, vector m contains the moles of species per unit volume of porous medium and vector c contains species concentrations in mol/mass of liquid (m =ϕρc for mobile species in a system of porosity ϕ and liquid density ρ). Matrix M is diagonal and its diagonal terms are unity when a species is mobile and zero otherwise; f is a source/sink term, which we use to represent chemical reactions. Here we define the linear operator L(c) in Eq. (1) as L(c) = −∇ · (ϕρvc) +∇ · (ϕρD∇c), where D is the diffusion/dispersion tensor and v is the fluid velocity, v = q /ϕ, q being Darcy's flux vector (in this application we assume that all mobile species are characterized by the same D and v). 2.2. Equilibrium reactions We start from the assumption that the time scales involved in the various chemical reactions occurring in the system are small with respect to typical diffusive and advective time scales. Under these conditions reactions can be considered at equilibrium. Moreover, we only consider aqueous reactions and precipitation–dissolution of minerals. We define f = STer, where r is the vector of reaction rates (expressed per unit volume of medium) and Se is the stoichiometric matrix of the chemical system (i.e., Se is an Nr × Ns

ð3Þ

with γa being the vector of activity coefficients. The latter are calculated in terms of the ionic strength, I, through the extended Debye–Hückel equation (Helgerson and Kirkham, 1974) logγi = −

pffiffi ∘ Az2i I ∘ pffiffi + b I: B + ai I

ð4Þ

Here, zi and åi are, respectively, the valence and the ∘ ionic radius of aqueous species i; b, A and B are temperaturedependent parameters, and I = 0.5∑iz2i ci. A reactive transport process is fully described once concentrations of the Ns species are obtained, together with the Nr reactions rates (Ns + Nr unknowns). In our case, this requires solution of the Ns coupled mass balance Eq. (1) and the Nr equilibrium Eq. (2). It is important to note that minerals need not be present, and precipitation/dissolution mechanisms need not occur, in the entire domain and at all times. Therefore, the definition of the chemical system (characterized by S′ea, c′a, c″a and K⁎) may change both in space and in time. 2.3. Components Solution of a reactive transport problem can be simplified upon defining components by introducing an auxiliary component matrix, U. Following De Simoni et al. (2007), we

30

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

define U so that USTe = (Sec 0)T. Because Se = (Sec|S′ea| − I), an expression for U is  U=

Uc Ua

Nc

 =

Nu

Nr

$ $ $! I

0

0

I

0 T S ′ea

l Nc : l Nu

ð5Þ

∂ðmc Þ T = Mc Lt ðmc Þ + Sec r; ∂t

ð6Þ

∂ðϕuÞ = Lt ðuÞ ∂t

ð7Þ

where vector mc contains the mass of constant activity species and u is a vector of Nu (aqueous) components defined as ′ Þ c″ u = Ua c = c′a + ðSea a: T

ð8Þ

Eq. (7) represents Nu transport equations. Vector u contains only aqueous species so that we can leave out matrix M in Eq. (7), while Mc in Eq. (6) is the part of M referring to constant activity species. As S′ea may vary in time and space, Nc, Nu and matrix U, may vary as well. If the definition of the chemical system does not change, the easiest way to solve Eq. (7) is to first solve u (with appropriate boundary and initial conditions). Then, (aqueous) solute concentrations are obtained from the expression of u in terms of concentrations of species (Eq. (8)) and the mass action law (Eq. (2)). If the chemical system changes in space and time, then one can solve the problem by substituting Eqs. (8) and (2) into (7) to obtain c′a from the resulting nonlinear equations. In both cases, equilibrium reaction rates are then calculated, e.g., from transport Eq. (1) of secondary species as ∂ðϕc″a Þ : ∂t

ð9Þ

The mass evolution of the constant activity species is finally obtained upon substituting r into Eq. (6). 2.4. Modification of medium properties Dissolution/precipitation processes may lead to changes in medium properties, specifically porosity and hydraulic conductivity, which have a major influence on groundwater velocity. To assess the importance of these effects, we include porosity variability in the numerical simulations of the experiments. Permeability (k) updating is performed for each time step on the basis of the Kozeny relationship (Bear, 1972) k = ko

ϕ3 ð1−ϕo Þ2 ð1−ϕÞ2 ϕ3o

ð10Þ

where ko and ϕo denote initial permeability and initial porosity, respectively. The porosity is calculated from the concentrations of minerals as ∂ϕ ∂mm = − ∑ vm ∂t ∂t m

3. Materials and methods 3.1. Solutions

Multiplying Eq. (1) by U leads to

r = Lðc″a Þ−

where subscript m refers to a mineral and v is the molar volume.

ð11Þ

Two potentially reactive solutions, calcium chloride and sodium carbonate, and a nonreactive solution, sodium chloride, at concentrations of 5 g/kg-water were used in the experiments. After weighing the salts and mixing with double deionized water, solutions were equilibrated with ambient air. Even mass concentrations were chosen to ensure equal densities between the solutions; this prevents stratification at the mixing zone of the two solutions, and allows creation of a quasi-two-dimensional flow and transport pattern within the investigated domain. The high concentrations, which are not typical of natural systems, were employed to increase the rate of reaction and the evolution of mechanical properties. 3.2. Flow cell The experiments were conducted in a rectangular laboratory flow cell filled with porous media, saturated with an initial aqueous solution of sodium chloride. The structure of the flow cell and the setup for the homogeneous packing experiments are illustrated in Fig. 1. The inner dimensions of the flow cell were 25 cm × 10 cm × 0.8 cm in the x-, y- and z-directions, respectively. The inlet face was split crosswise to create two separate inlets, side by side, each with a length of 5 cm. In the reactive transport experiments, calcium chloride and sodium carbonate solutions entered the cell through inlets 1 and 2, respectively. In the conservative transport experiments, sodium carbonate was replaced by the initial sodium chloride solution in inlet 2. The inlets were designed to allow solutions to enter evenly along the entire length of the inlet. The outlet face includes 12 regularly spaced 4.4 mm diameter holes. Ten additional holes were bored in the top face of the cell and sealed with septa, creating sampling ports. BD micro-fine syringes (0.3 mL) were used to sample water from the inside of the flow cell through the septa. During the experiments, minimal volume samples of 50 μL were withdrawn. Calcium concentration was then measured in the samples taken from the cell, as detailed in Section 3.5. The inner space of the flow cell was viewed in the experiment as a two-dimensional domain with an origin in the middle of the inlet face (between inlet 1 and inlet 2, as shown in Fig. 1). Legends and locations of the sampling ports are given in Table 1. It is important to note that sampling from the port interferes with the flow field and therefore also with the spatial/temporal distribution of concentration throughout the cell. To minimize this interference we reduced as much as possible the sampling volume and set large time gaps between samplings, which allowed for more than one pore volume of solution to flow through the cell. The total duration of the experiments was equivalent to 9.6 and 10.8 pore volumes for the situations corresponding to the packing of the cell with homogeneous or heterogeneous materials, respectively.

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

31

Fig. 1. Experimental setup for experiments with homogeneous packing.

3.3. Homogeneous packing experiments — setup and procedure Two columns provided separated inlet reservoirs for the solution used in the experiment. The two reservoirs were connected via the flow cell, thus creating equal pressure conditions at the two inlets. As the water levels in the two reservoirs descended simultaneously, the two column system ensured inlets with even volumetric flow rates. At the outlet boundary, water was pulled from the cell through the outlet holes with the use of a KDS 230 multi-syringe pump. The pump was modified to carry twelve 60 mL Teromo syringes. The setup was designed to create pseudo-one-dimensional flow, parallel to the x-axis. Although the pressures at the boundaries changed during the experiments, the global hydraulic conductivity of the cell remained essentially constant (even for the reactive case, because the relative volume of precipitate was very small and the precipitation strip was generally oriented parallel to the flow. Therefore flow in the system can be modeled by keeping fixed heads at the inlets, and fixed uniform flow at the outlet. The flow cell was packed with 1 mm diameter glass beads with a measured porosity of 0.375. A flow rate of 12 mL/h was set at the outlet syringe pump creating an average water velocity of ∼97 cm/day inside the flow cell. The duration of the reactive transport experiments was 60 h. During this period, sampling from the ports located inside the flow cell was performed seven times. The duration of conservative transport experiments was limited to 30 h, because there was no evolution of mechanical properties (i.e. porosity and

permeability) in the medium throughout an experiment, and the system reached a steady state relatively quickly. Sampling was performed four times during each conservative transport experiment. The small number of sampling was employed to ensure minimum interference of the flow field. We note that the goal of this experiment was not to capture the replacement of one solution by another but rather to focus at much longer times where processes like mineral precipitation influence the spatial distribution of concentrations. Prior to the experiments the establishment of the desired flow field was verified with color pulse experiments. Red and blue dye solutions were injected at the inlets and the evolution of color speciation was observed under the same flow conditions used in the conservative and reactive transport experiments. Reactive and conservative transport experiments were repeated three times and twice, respectively, to ensure repeatability of measurements. Between reactive experiments the cell was cleaned and repacked. Between conservative experiments the cell was washed with several pore volumes of fresh water and initial condition (sodium chloride) solution. 3.4. Heterogeneous packing experiments — setup and procedure Fig. 2 illustrates the setup for the heterogeneous domain experiments. As the establishment of a pseudo-one-dimensional flow was not relevant in this case, a new setup was employed that featured simpler laboratory operation and focused on the effect of nonuniform flow conditions. Peristaltic pumps with

Table 1 Legends and coordinates of the sampling ports. Sampling ports

A1

A2

A3

A4

C1

C2

C3

D1

D2

D3

x coordinate y coordinate

11.4 − 1.5

11.4 − 0.5

11.4 0.5

11.4 1.5

17.4 −1

17.4 0

17.4 1

21.4 −1

21.4 0

21.4 1

Coordinates are given in centimeters, with the origin centered at the inlet face between inlets 1 and 2 (see Fig. 1).

32

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

Fig. 2. Experimental setup for experiments with heterogeneous packing.

flow dampeners were used to set the two inlet flow rates at 0.6 mL/min. All outlet holes were sealed, except for the one located 0.55 cm along the y-axis. This hole was left open to atmospheric pressure. Sand with 30/40 mesh size and average grain size of 0.532 mm (Levy and Berkowitz, 2003) was used as a porous medium, in addition to the 1 mm glass beads that were used in the homogeneous transport experiments. The sand was crushed, sieved and cleaned by UNIMIN Corporation, USA. It is composed of well-rounded quartz particles with minimal surface coating (99.8% pure SiO2 as reported by UNIMIN). We measured the porosity to be 0.32. The sand filled a trapezoidalshaped area inside the cell with corners at the x and y coordinates (in cm) (14, 0), (21.25, 0), (14, 10) and (15.25, 10). The remainder of the domain was packed with the same glass beads used in the homogeneous experiment. Uniformity of the packing along the z-direction was maintained to allow two-dimensional analysis of the domain. An illustration of the packing is given in Fig. 2. Conservative and reactive transport experiments in the heterogeneous packing arrangement were repeated three times and twice, respectively. Between conservative transport experiments, the cell was washed with several pore volumes of fresh water and initial condition solution. Between reactive experiments, the cell was washed with diluted acid followed by fresh water and initial condition solution. 3.5. Analytical method Calcium concentrations were measured using a colorimetric method based on a method described by Pierce Chemical Company (1973). A basic pH solution with 0.1 mM of methylthymol blue metal-binding dye was used as a reagent. Absorbance was read at 611 nm wave length using a Cary 100 UV–Vis spectrophotometer. A calibration graph based on spectrophotometric measurements of standard solutions was used to determine actual calcium concentrations. 3.6. Numerical modeling We modeled the transport experiments using the code Retraso Code Bright (RCB) (Saaltink et al., 2004), which was designed to model multiphase reactive transport processes and accounts for effects of mineral precipitation and dissolution.

RCB calculates flow and reactive transport separately. Therefore, porosity and permeability required for the flow equation are updated from the concentrations of minerals that are calculated by reactive transport at the previous time step. Concentrations of minerals are calculated for each node, whereas the permeability is required for each element. Appendix A provides details of the algorithm used to transform the nodal concentrations of minerals to porosities within each grid element. Permeabilities are then calculated from porosities using Eq. (10). Calcium carbonate was allowed to precipitate only in the form of calcite and all reactions were modeled assuming local chemical equilibrium conditions. A constant temperature of 25 °C was imposed during the simulations. The chemical compositions of the initial solutions used in the experiments were calculated using the code PHREEQC (Parkhurst and Appelo, 1999). Calculations were performed under the assumption of equilibrium with partial CO2(g) pressure of 3.16 × 10− 4 bar. PHREEQC also provides the relevant aqueous complexation reactions. The complete chemical system is formed by 19 chemical species and 12 reactions. All reactions and their equilibrium constants are displayed in Table 2. Table 3 presents the values of ionic radius used for the calculation of the activity coefficients according to Eq. (4). The homogeneous medium experiment was modeled using hydraulic boundary conditions of constant atmospheric pressure along the inlet face and constant volumetric flux

Table 2 Chemical reactions with their log equilibrium constants, log K, used in the reactive transport simulations. Unit activity is assumed for water. A uniform temperature, T = 25 °C, is adopted. Reaction R1 R2 R3 R4 R5 R6 R7 R8 R9 R10 R11 R12

+

log K 2+

−

CaCl = Ca + Cl CaCl2(aq) = Ca2+ + 2Cl− − CaCO3(aq) + H2O = Ca2+ + HCO− 3 + OH 2+ CaHCO+ + HCO− 3 = Ca 3 CaOH+ = Ca2+ + OH− CO2(aq) + OH− = HCO− 3 − − CO2− 3 + H2O = HCO3 + OH + − H + OH = H2O + − − NaCO− 3 + H2O = Na + HCO3 + OH NaHCO3(aq) = Na+ + HCO− 3 NaOH(aq) = Na+ + OH− 2+ CaCO3(s) + H2O = HCO− + OH− 3 + Ca

0.6956 0.6436 − 6.9934 − 1.0467 − 1.1451 7.6504 − 3.6663 13.9951 − 4.1807 − 0.1541 0.1849 1.8487

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

33

Table 3 Ionic radius of the system species. Ion

Ca2+

Cl−

OH−

HCO− 3

Na+

NaCO− 3

H+

CO2− 3

CaOH+

CaHCO+ 3

CaCl+

Ionic radius (Å)

6.0

3.0

3.0

4.0

4.0

4.0

9.0

5.0

4.0

4.0

4.0

along the outlet face. The heterogeneous case was simulated using hydraulic boundary conditions of constant volumetric flux along the inlet face, and constant atmospheric pressure at the narrow outlet. Transport for both homogeneous and heterogeneous cases was solved using an inlet boundary condition of prescribed mass flux, and a free flow outlet boundary. Initial values of permeability (k0) were estimated by the Kozeny–Carman model

k0 = 8:3 × 10

−3 2 d10

ϕ30 ð1−ϕ0 Þ2

ð12Þ

where d10 and ϕ0, respectively, are the 10% cumulative passing particle diameter and medium initial porosity. We obtained k0 = 1.12 × 10− 9 m2 and k0 = 1.47 × 10− 10 m2, respectively for glass beads and sand; these values are close to those estimated from laboratory measurements. It is important to note that for the homogeneous experiments, the initial value of permeability does not affect calculation of the flow field in the numerical simulation. On the other hand, different estimates of permeabilities of glass beads and sand, using variations of Eq. (12), may affect the calculated flow field in the heterogeneous domain. For chemical transport, we assumed that all chemical species are characterized by the same molecular diffusion and dispersivity coefficients. The former was set to a value D0 = 10− 9 m2/s; the dispersivity was determined by a best fit of the numerical simulations against concentrations measured at the ports within the cell during conservative tracer experiments. RCB uses a finite difference scheme to discretize the time derivatives, and a finite element scheme for the spatial derivatives. To achieve high resolution in the mixing zone, and for compatibility with code limitations on the maximum number of grid nodes, we performed a preliminary sensitivity analysis on the grid size and influence of boundary conditions. Final simulations of the homogeneous packing system were carried out using a smaller domain, i.e., only the central 4 cm in the y-direction and the first 22 cm in the x-direction were modeled. Homogeneous packing simulations, for the conservative transport experiments, were based on a symmetric grid of 220 × 61 elements with uniform 1 mm length (along the x-direction). Elements width (aligned with the ydirection) decreased from 2 mm for the grid cells located at the largest distances from the system centerline to 0.1 mm in the middle of the domain, where precipitation was observed. In the heterogeneous packing simulations, a 200 × 80 grid was used with elements of 1.25 mm length, and width between 0.4 mm (in the region where, according to preliminary simulations, precipitation was seen to occur) and 5 mm (near the impermeable boundaries). Further discussion of the calibration of dispersivity is given in Section 4.

4. Results and discussion 4.1. Homogeneous packing 4.1.1. Conservative transport and calibration of dispersivity Calibration of the dispersivity was performed on the basis of the conservative transport experiments. Preliminary simulations confirmed that the longitudinal dispersion (DL) influences the results significantly only at early times (i.e., less than 10 h), during displacement of the fluid which is initially in the flow cell. Transverse dispersion (DT) has an important role because it governs mixing processes in the system. We define transverse dispersion as DT = αTv + D0, where αT, v and D0 denote transverse dispersivity, average pore velocity and diffusion coefficient, respectively. We assigned a typical value of D0 = 10− 9 m2/s, as described in Section 3.6, and simulated the homogeneous packing, conservative experiment, for several αT values. We compared the simulated calcium component concentrations in domain locations that match the locations of the sampling ports for the measured concentrations, and found the overall best fit value αT = 0.2 mm, regardless of the value of αL. We observed that there was little or no change in simulation results over a large range of longitudinal dispersivity values. Fig. 3 depicts the comparison between model results (solid lines) and the measured calcium concentrations (symbols) at the sampling ports within the flow cell (Fig. 1). The model results are associated with the nodes at the exact location of the sampling port (where possible), or with a linear interpolation between the concentration values calculated at the two closest nodes. Maximum concentration refers to the calcium concentration at the calcium chloride inlet boundary. It is difficult to estimate measurement uncertainty and error. Because of the difficulty in attempting to estimate an error range for the measured concentration values, we present all of the measurements in the figures. The error in concentration measurements by UV–vis spectrophotometry is negligible, while meaningful estimation of uncertainty requires numerous repetitions of the experiment. However, at least in qualitative terms, we consider the agreement to be of acceptable quality and do not pursue this point further. Fig. 3a, c, e, and g corresponds to sampling ports located in the “calcium side” of the cell, while Fig. 3b, d, f, and h represents ports in the “carbonate side”, located symmetrically with respect to the centerline of the cell (see Fig. 1). It can be seen from Fig. 3 that the model predicts reasonably accurate concentrations for ports A1 (Fig. 3a) and A2 (Fig. 3c). For port A3, the model predicts concentrations higher than those measured, whereas the opposite behavior is observed at port A4. Both of these ports are in the carbonate side of the system. Therefore, a change in the dispersivity value would improve the fit of one plot but degrade that of the other. The poorer fits to the A port data may be due to the proximity to the

34

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

Fig. 3. Calibration of dispersivity. Total aqueous calcium concentrations in homogeneous packing, conservative transport experiments, at locations of the sampling ports. Model predictions are indicated by a solid line. Results from measurements of the first, second and third experiments (1, 2 and 3) are denoted by the ×, triangle, and circle symbols, respectively. The dashed line represents the total calcium concentration in the calcium chloride boundary solution.

inlets, where the sampled volume is located in a region of relatively high transverse concentration gradients; the presence of local heterogeneities may also contribute to this effect. Fig. 3e–h, representing ports C and D, show good agreement between measurements and simulations on both sides of the cell. One of the concerns that arises when sampling from ports is that the small volume of water extracted from the domain might not accurately represent the water chemistry at the

location of the port. In this case, the differences between model and experimental results will increase in regions in which the concentration gradient is large (e.g., the calcium gradient near the inlets along the y-axis). This difficulty is reflected in the quality of the agreement between model and experimental results in ports C and D, compared to A. The evolution of the modeled concentrations can be described as having two phases: (i) chemical breakthrough phase, characterized by a

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

sigmoid shaped curve (this phase ended at all ports after about 8 h); and (ii) a steady state phase characterized by a plateau. We note that measured calcium concentrations at ports A2, C1, and D1 did not reach the value of the input solution because of mixing between the two waters injected in the cell. Overall, a model based on the ADRE appears to capture the spatial distribution of total calcium concentrations for the homogeneous, conservative case.

4.1.2. Reactive transport Calcium carbonate precipitated along the mixing zone between the two solutions containing CaCl2 and Na2CO3, which flowed in parallel during the reactive transport experiments. The formation of a white strip of calcium carbonate precipitate along the cell centerline can be seen clearly in Fig. 4, which shows the domain 30 h after the start of a reactive transport experiment. The precipitation can already be seen along the entire centerline of the cell, although less calcium carbonate precipitated near the outlet, where the strip of precipitation is seen to be weaker. Another feature of the precipitation region that can be seen in Fig. 4 is the mild shift to the “carbonate side” of the cell (see detailed explanation below accompanying Fig. 6). The precipitation of calcium carbonate in the cell also causes an evolution in the mechanical properties of the medium, as pores in the mixing zone are filled by the precipitating mineral. Fig. 5 compares the experimental and simulated calcium concentrations at the sampling port locations for the reactive scenario. We begin our discussion by focusing on the sampling ports located in the calcium side of the domain (Fig. 5a, c, e, and g). The experimental results display very good repeatability and can be described as follows. During a first phase one can distinguish an increase in concentrations. It is difficult to determine from the limited measurements whether this rise is more properly described by a linear increase than by a sigmoid. The rising phase ended when the calcium concentrations at the sampling ports reached the values of the calcium chloride boundary solution. Concentrations remained at this level until the end of the experiment. The duration of the first phase was about 28 h for the A1 port and 38 h for the A2, C1 and D1 sampling ports. For comparison, note that one pore volume corresponds to about 6.2 h.

Fig. 4. Calcium carbonate precipitate along the centerline of the flow cell, 30 h after initiation of a preliminary experiment. The cell is packed with homogeneous media.

35

Comparing the measured concentrations after 13 h at ports A1, A2, C1, and D1 (Fig. 5a, c, e, and g) to the corresponding values detected for the conservative scenario (Fig. 3a, c, e, and g), one can note that concentrations were significantly lower in the reactive case than in the conservative scenario at all sampling ports. This is associated with the fact that calcium is subtracted from the solution in the mixing region as a direct result of calcium carbonate precipitation. The rise of sampled concentrations to those of the boundary water at ports A2, C1, and D1 indicates a reduction in the reaction rate in the region between the inlet face and the ports, and also a major reduction in calcium migration towards the carbonate side of the cell due to dispersion and diffusion. An explanation for this behavior is clogging of pores in the mixing zone. The evolution of the precipitation region from the inlets toward the outlet, as seen in Fig. 4, is consistent with the observed gradual arrival of concentration to the maximal value in the different ports (first A, and then C and D). The time scale of the evolution of the precipitation line, as observed visually in the experiment (tens of hours), is consistent also with the time required for the measured concentrations (in Fig. 5) to reach their final value. Similarly to the conservative transport case, the simulation results for the reactive scenario observed at ports A1, A2, C1, and D1 can be described by two phases. The first phase corresponds to a breakthrough of concentration, which does not last longer than 8 h. Two aspects can be noticed about the calculated concentrations reached at the end of this first phase, in all of these sampling ports: (i) they are lower than the corresponding concentrations simulated in the conservative transport scenario because of the occurrence of precipitation; and (ii) they are higher (almost double) than the measured calcium concentrations. The latter might indicate that the rate of calcium precipitation rendered by the model is lower than the real rate of precipitation. The second part of the curve displays a steady, moderate rise of concentration that persists until the end of the experiment. The positive slope, in contrast to the steady state phase observed in the conservative transport simulation, implies a continuous, albeit slow, reduction in the degree of mixing and the rate of precipitation. This could be due to clogging effects. However, comparing the slope of the modeled curve to the steeper slope of the experimentally measured concentrations, it is clear that the behaviors predicted by the model occur much slower than what was observed in the laboratory. Indeed, even after 60 h, the concentrations in ports A1, A2, C1 and D1 were significantly lower than the concentration of the calcium chloride inlet boundary solution. Fig. 5b, d, f, and h (ports A3, A4, C3, D3 which are located in the half of the domain where sodium carbonate is injected), provides a view of the system dynamics that is somewhat noisier, but still in agreement, with the measurements obtained from the calcium side ports. Focusing on the C3 and D3 sampling ports, it is seen that the concentrations measured after 13 h in the experiment were much lower than those measured after 6 h. This sharp decrease in concentration may be attributed to the influence of the initial sodium chloride solution (as discussed below), in addition to a clogging effect. At later times the clogging effect was observed as the concentrations continued to decrease until asymptotically reaching a minimal value of about 2–3 × 10− 5 (mol/kgwater) after about 40 h. For some experiments this value is

36

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

Fig. 5. Total aqueous calcium concentrations in homogeneous packing, reactive transport experiments, at locations of the sampling ports. Model predictions are indicated by a solid line. Results from measurements of the first, second and third experiments (1, 2 and 3) are denoted by the ×, triangle, and circle symbols, respectively. The dashed line represents the total calcium concentration in the calcium chloride boundary solution.

not significantly different from zero concentration. Tartakovsky et al. (2008) claimed that if the precipitation separates the two solutions, then they may become undersaturated (with respect to calcite) and subsequently begin to dissolve the calcite. It is possible that such a mechanism occurred in our experiment, and was responsible for the low concentrations of calcium observed at the C3 and D3 ports at late times. However, the level of accuracy of our measurements cannot serve as a basis for this discussion.

The concentration measurements from the A4 port (Fig. 5b) showed different trends between repetitions of the reactive experiment. Nevertheless, all of the measurements were within a small range of concentrations, of about 10− 4 (mol/kg-water). The measured concentrations from the A3 ports did not show a very good repeatability. This could have been the result of measuring at a location with high proximity to the precipitation line. In spite of that, it appears that the A3 concentration values at about 6 h after the beginning of the

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

experiment were higher than those observed during the rest of the experiment, and that the minimum measured concentration did not reach zero. The numerical results displayed in Fig. 5b, d, f, and h suggest that the relatively high concentrations measured after about 6 h from the beginning of the experiment were not experimental errors. A spike is observed in the simulation curves in all of the carbonate side plots. This spike simply indicates that the highest concentration of calcium at these points (ports A4, A3, C3, and D3) appeared when the two boundary solutions were still very much diluted by the initial sodium chloride solution. Once the initial solution was flushed from the flow cell, more calcium precipitated in the region between the domain inlets and sampling ports and less reached the sampling ports. Shortly after the breakthrough of the solutions, the concentrations in the locations corresponding to A4, A3, C3, and D3 decreased and stabilized at a minimal value that is determined mostly by the rate of the precipitation reaction. Although the concentrations predicted by the model at 60 h agrees relatively well with the measurements, the processes that drive these results are different (as explained above). To complete the picture, Fig. 6a and b displays, respectively, the calculated calcite accumulation (in logarithmic scale) and the calculated porosity distribution within the domain after 60 h from the beginning of the experiment. Concentrations near the inlets reached values of the order of 104 mol/m3. These large values dropped by orders of magnitude within a few centimeters along the direction of flow; it is seen that calcite precipitated all along the cell. We observe that a large reduction in porosity is predicted near the inlets, with values reaching a minimum of the order of 10− 3. No change in porosity was predicted from 10 cm along the cell up to the outlet face. This is in contrast to the observations

37

from Fig. 4. Although no quantitative results on porosity and calcite distributions can be inferred from Fig. 4, it is reasonable to assume that if precipitation is detected at some locations, there is also some local reduction in porosity. The precipitation strip in the simulation did show some deviation towards the carbonate side of the flow cell, possibly due to the fact that calcium carbonate solubility increases with acidity. The pH in the sodium carbonate solution was higher than that of the calcium chloride (9.8 as opposed to 5.6, respectively), so that the best conditions for precipitation occur when the mixing ratio is to a small extent in favor of the sodium carbonate. In summary, our computational results suggest that the model recognizes some of the key processes that determine the evolution of calcium concentrations in the flow domain, but underestimates their magnitude. In other words, the model reproduces the observed changes in concentrations, but not their rates. One reason for this may be the transformation of nodal concentrations of minerals to elemental porosities. The algorithm described in Appendix A smears the strip of reduced porosity and thus reduces the clogging effect. Moreover, the concentrations of minerals are not bounded by the volume available for precipitation. Therefore calcite concentrations could (in principle) reach values from which a straightforward application of Eq. (11) would yield a negative (nodal) porosity. 4.2. Heterogeneous media 4.2.1. Conservative transport We devote this Section to investigate the ability of the ADRE-based transport model (as implemented in the algorithms of RCB) to simulate the observed dynamics of the flow and precipitation processes in a heterogeneous medium. A

Fig. 6. Simulation results of reactive transport in homogeneous packing after 60 h (equal to the entire duration of the experiment). (a) Calcite accumulation displayed as the logarithm of the concentrations; (b) spatial distribution of porosity.

38

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

Fig. 7. Preliminary color pulse experiment in the heterogeneous packing. Dye solutions are concentrated near the outlet on the right.

trapezoidal-shaped block of sand was located in the middle of the glass bead domain (see Fig. 2). The outlet boundary was modified to permit flow only through a small outlet. This allowed establishment of a nonuniform flow pattern in the flow cell and had a significant effect on the shape of the mixing zone and the distribution of chemical elements throughout the domain. A tracer experiment was performed by using two different dye solutions, to obtain a preliminary assessment of the flow pattern within the cell. As an example, Fig. 7 depicts the spatial distribution of dyes after 1 h. We note that the relatively sharp mixing zone between the blue and red solutions formed a curved interface that ended near the narrow outlet where the two solutions converged. The shape of this interface was a direct outcome of the velocity field. The latter was, in turn, determined by the boundary conditions and the hydraulic conductivity field. These observations were complemented by the depiction of the calculated streamlines, as reported in Fig. 8 (the gray area marks the sand inclusion). The pattern of streamlines is consistent with the observed shape of the mixing region between the two solutions. Fig. 9 describes the comparison of modeled and experimentally measured calcium concentrations under conservative transport conditions. Because dispersivity is often said to be correlated to grain size (Dullien, 1979), a ratio of 2:1 was maintained for the dispersivities of glass beads and sand, respectively. Using the calibrated dispersivity of the system, (i.e., αT = 0.2 mm), we attempted to simulate the heterogeneous packing experiments with transverse dispersivities of 0.2 mm and 0.1 mm for glass beads and sand, respectively. However, the code did not converge to a solution, even with variation of the time step and with further reduction in the

grid size. We believe that the lack of convergence is caused by the increase in the Peclet grid number (Pl) associated with the y-direction. Discussion regarding the relation between Pl and the stability of the numerical scheme is beyond the scope of this work. The interested reader is referred to the work of Kinzelbach (1986). In practice the values of αT were increased by small increments, manually, until a solution was achieved; convergence was obtained for αT = 0.8 mm and αT = 0.4 mm for glass beads and sand respectively. The modeled calcium concentration are slightly lower than most measured values at the sampling ports that are relatively deep in the calcium side of the solutions interface, i.e., at ports A1, C1, C2, and D1 (Fig. 9a, e, f, and h). After the breakthrough phase, concentrations at A1, C1, and D1 stabilized at a concentration similar to the calcium chloride boundary solution. On the other hand, it can be seen that calculated concentrations at ports A4 and D3, which are located on the carbonate side of the cell, were much higher than their measured counterparts. It is possible that the large difference between measurements and predictions appearing at port A4 is due not only to the use of values of αT that are too large, but also to small deviations between the computed and the actual flow fields, which can impact the calculated concentrations. There is a good match between model and measurements for ports C3 and D2 (Fig. 9g and i), which are close to the interface between the two solutions; on the other hand, at port A3 (Fig. 9d) which is also very close to the interface of the solutions, the modeled and measured concentrations display large relative and absolute differences. 4.2.2. Reactive transport Fig. 10 depicts the measurements and model results for the reactive transport experiment. The figure shows that concentration at the A1, C1, and D1 (Fig. 10a, e, and h) ports stabilized at the maximal concentration at the end of the breakthrough phase, similarly to the observations for conservative transport. However, the model renders lower concentrations at these ports, which can be attributed to the choice of a large transverse dispersivity value which is too large (as well as the possibility that the ADRE does not fully account for the dynamics of the system). Changes in the concentrations caused by an evolution of the media properties were shown by the moderate rise in measured concentrations after the breakthrough phase at ports C2, C3 and D2 (Fig. 10f, g, and i). In comparison to C2, the subplots associated with ports D2 and C3 show a good match between

Fig. 8. Streamlines in the heterogeneous packing, calculated from fluid velocities in a simulation of conservative transport.

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

39

Fig. 9. Total aqueous calcium concentrations in heterogeneous packing, conservative transport experiments, at locations of the sampling ports. Model predictions are indicated by a solid line. Results from measurements of the first, second and third experiments (1, 2 and 3) are denoted by the ×, triangle, and circle symbols, respectively. The dashed line represents the total calcium concentration in the calcium chloride boundary solution.

model and measurements, in terms of the values of concentration at which the breakthrough phase ends. We remark that ports D2 and C3 are closest to the precipitation region from the side of the inlet calcium solution. The fact that the measured concentrations at D2 and C3 did not reach the maximum boundary concentration indicates that the pores in the precipitation region were not completely clogged.

Similarly to the homogeneous case, the temporal evolution of concentrations at the sampling port locations, which is due to clogging, is underestimated by the model. This can be deduced from the moderate slope in the second phase of the modeled curves at ports C3, D2, and A2, in comparison to the measured values. The concentrations measured at port A4 (Fig. 10b) were low but showed a slow decreasing trend that

40

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

Fig. 10. Total aqueous calcium concentrations in heterogeneous packing, reactive transport experiments, at locations of the sampling ports. Model predictions are indicated by a solid line. Results from measurements of the first, second and third experiments (1, 2 and 3) are denoted by the ×, triangle, and circle symbols, respectively. The dashed line represents the total calcium concentration in the calcium chloride boundary solution.

is most likely caused by the clogging effect. The concentrations predicted by the model are lower than the experimental measurements by a factor of five for early time measurements. Note, however, that because the concentrations are so low, any constant measurement error (which is difficult to estimate) could change this factor. For port D3 (Fig. 10j), the model predicted concentrations that were one order of magnitude higher than the measured values. However, a decrease in concentrations that is related to the clogging is seen clearly. A similar trend was seen in the experiments,

although this is not visible from the figure because of the vertical scale. Fig. 11 shows the flow cell after 11 h of the reactive experiment. A white strip of calcium carbonate precipitate appeared approximately where the mixing region between the dye solutions appeared in the preliminary experiment (Fig. 7). Fig. 12a shows the calculated cumulative calcium carbonate distribution at 11 h from the beginning of the experiment. A white area appears in the figure where no calcium carbonate precipitation was predicted by the model.

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

41

Fig. 11. Calcium carbonate precipitate at the end of a reactive transport experiment with the heterogeneous packing. Shown here is a mirror image of the bottom face of the cell (sampling ports are not visible).

Fig. 13b and c shows the simulated distribution of porosity in the domain after 11 h and 22 h, respectively. Qualitatively, the adopted model predicts the occurrence of a precipitation strip which resembles the one observed in Fig. 11, in terms of

both shape and location. Relatively high accumulation of calcium carbonate is predicted near the inlets and the outlet. High accumulation is indeed expected near the inlets, where mixing occurred from the beginning of the experiment.

Fig. 12. Simulation of reactive transport in heterogeneous packing. (a) Distribution of the logarithm of CaCO3 concentrations (expressed in mole per cubic meter) after 11 h (equal to the duration of the experiment); (b) Spatial distribution of porosity after 11 h; (c) Spatial distribution of porosity after 22 h.

42

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

Moreover, one should also take into account that the solutions at the inlet side of the flow cell were not yet depleted from the precipitating elements. On the other hand, the convergence of streamlines near the outlet created enhanced mixing that caused high accumulation of calcium carbonate. Visual inspection of Fig. 12 indicates that the model predicts considerably less precipitation in the sand area than in the glass bead region. This is most likely due to our choice of dispersivity values, i.e., we adopted a lower dispersivity for the sand than for the glass beads (although this effect could also result from the different pore sizes in the two types of porous media). Fig. 11 suggests that the precipitation strip was less wide inside the sand than inside the glass beads area. However, no quantitative analysis was performed to corroborate the visual inspection. The porosity decrease predicted by the model after 11 h and even 22 h is lower than what was observed and suggested from Figs. 11 and 12. For example, inside the sandy region, the model hardly predicts any visible change in porosity even though a solid precipitation region is clearly visible in the experiments throughout the entire section of sand.

reactions, and/or to the need to transform the mineral concentrations which are typically calculated at grid nodes to porosities at elements. The latter is needed because permeability depends on porosity and, while species concentrations are calculated at nodes, permeability values are assigned to elements. This transformation is a common problem for all Eulerian numerical methods (e.g., finite elements, finite differences), as it tends to smooth the zones with reduced porosity; as a result, the transformation yields regions of reduced porosity that are wider, but with smaller reductions in local porosity. Only use of very fine grids may solve this problem. The measurements provided in this study are unique, and can provide a benchmark for further analysis of modeling tools, considering both ADRE-based and other modeling approaches. In future experimental investigations, additional information could be extracted by employing an analytical method that can measure more than one chemical element in a sample. For example, measurement of sodium in the reactive transport experiments and comparison to the conservative case can provide a view of the spatial evolution of concentration that is affected only by temporal changes in the hydraulic properties.

5. Concluding remarks Acknowledgements We performed a set of experiments to obtain direct, detailed measurements of reactive transport and precipitation in a flow cell containing both homogeneous and heterogeneous porous media. These data were used as a benchmark to enable semi-quantitative comparison to the results of a reactive transport numerical model. Our work leads to several main conclusions. Multiple experiments on all four systems, conservative and reactive transport in both homogeneous and heterogeneous domains, showed reasonable reproducibility. A feature of the adopted setup is that it enabled spatial measurements of concentration, which could be coupled to a numerical model to estimate transverse dispersivity and capture some of the relevant space–time dynamics of a reactive transport scenario. We found very good agreement between the ADRE-based model and the experiments in the case of conservative transport in homogeneous media. Due to limitations of the numerical scheme, the simulations of the heterogeneous media scenarios employed a larger value of transverse dispersivity relative to that calibrated for the conservative homogeneous case. This might have contributed to reducing the accuracy of the predicted distributions of spatial concentrations. In spite of this limitation, the model coped well with the introduction of complexity in the form of a heterogeneous medium. In the two reactive transport scenarios, and notwithstanding the limitations of applying an ADRE-based model (e.g., the inherent assumptions of scale separation and domain homogeneity), the model mimicked the precipitation process and the trends in evolution of calcium concentrations that occurred in the domain. However, the model severely underestimated the clogging phenomenon and its influence on concentration distributions throughout the reactive experiments. This may be due to the use of an ADRE-based continuum model in the presence of pore-scale

The financial support of the Israel Science Foundation (grant 575/08) is gratefully acknowledged (GK and BB). Appendix A. Algorithm to transform concentration of nodes to porosities within elements Porosity at elements is calculated from concentrations of minerals at nodes. Each element j is subdivided into Nn (elementj) patches of equal volume, V(patch)ij, Nn(elementj) being the number of nodes of element j (see Fig. A1). Each node i is then associated with a volume (or cell) of influence, V(cell)i, given by

Vðcelli Þ =

Ne ðnodei Þ

∑

j=1

Vðpatchij Þ

ðA1Þ

where Ne(nodei) is the number of elements sharing node i. Updating the porosity of an element from the nodal mineral concentrations is accomplished upon distributing the

Fig. A1. Sketch depicting the volume (or cell) of influence, V(cell)i, associated with node i, shared by different elements in the finite element mesh.

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

porosity calculated at nodes (from mineral concentrations) to the elements shared by the node according to Nn ðelementj Þ

∑

i=1

Δϕðelementj Þ =

 f ðpatchij ÞΔVp ðnodei Þ Vðelement j Þ

ðA2Þ

where Δϕ(elementj) is the change of porosity in element j, ΔVp(nodei) is the change in pore volume that is calculated from precipitation and dissolution as Nm

ΔVp ðnodei Þ = − ∑ Vðcelli ÞΔmm ðcelli Þvm m=1

ðA3Þ

where subscript m refers to the mineral, Δmm is the change in concentration of mineral m (in mol per volume of porous medium), calculated through the reactive transport equations for each node i, Nm is the number of minerals and vm is the molar volume of the mineral. The term f(patchij) in Eq. (A2) is a distribution factor for each patch ij. The latter is calculated taking into account the volume of pores at the previous time step, Vkp,  f ðpatchij Þ =

n Vpk ðpatchij Þ

Ne ðnodei Þ 

∑ j

n Vpk ðpatchij Þ :

ðA4Þ

Here n is an exponent that assigns more or less weight to the porosity of the previous time step. When n tends to infinity, all pore reduction occurs in the element with the highest porosity, whereas a value of zero does not take into account at all the porosity value calculated during the previous time step. The volume of pores Vkp is calculated as Vpk ðpatc hij Þ =

Vðelement j Þ Nn ðelement j Þ

k

ϕ ðelementj Þ i = 1; …; Nn ðelementj Þ

ðA5Þ where ϕk(elementj) is the porosity of element j at the previous time step. Note that the sum of distribution factors for each cell equals unity Ne ðnodei Þ

∑ j

f ðpatchij Þ = 1:

ðA6Þ

References Back, W., Hanshaw, B.B., Pyle, T.E., Plummer, L.N., Weidie, A.E., 1979. Geochemical significance of groundwater discharge and carbonate solution to the formation of Caleta Xel Ha, Quintana Roo, Mexico. Water Resources Research 15, 1521–1525. Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York. 764 pp. Bolton, E.W., Lasaga, A.C., Rye, D.M., 1996. A model for the kinetic control of quartz dissolution and precipitation in porous media flow with spatially variable permeability: formulation and examples of thermal convection. Journal of Geophysical Research-Solid Earth 101 (B10), 22157–22187. Bolton, E.W., Lasaga, A.C., Rye, D.M., 1997. Dissolution and precipitation via forced-flux injection in a porous medium with spatially variable permeability: kinetic control in two dimensions. Journal of Geophysical Research-Solid Earth 102 (B6), 12159–12171.

43

Cederberg, G.A., Street, R.L., Leckie, J.O., 1985. A groundwater mass-transport and equilibrium chemistry model for multicomponent systems. Water Resources Research 21 (8), 1095–1104. Cirpka, O.A., Olsson, A., Ju, Q.S., Rahman, M.A., Grathwohl, P., 2006. Determination of transverse dispersion coefficients from reactive plume lengths. Ground Water 44 (2), 212–221. Cochepin, B., Trotignon, L., Bildstein, O., Steefel, C.I., Lagneau, V., Van der lee, J., 2008. Approaches to modelling coupled flow and reaction in a 2D cementation experiment. Advances in Water Resources 31, 1540–1551. Colorimetric Determination of Calcium in Biologic Fluids. Pierce Chemical Company (Rockford, IL), United States. Corbella, M., Ayora, C., Cardellach, E., 2003. Dissolution of deep carbonate rocks by fluid mixing: a discussion based on reactive transport modeling. Journal of Geochemical Exploration 78–9, 211–214. Daccord, G., Lenormand, R., Lietard, O., 1993. Chemical dissolution of a porousmedium by a reactive fluid 1. Model for the wormholing phenomenon. Chemical Engineering 48 (1), 169–178. De Simoni, M., Carrera, J., Sanchez-Vila, X., Guadagnini, A., 2005. A procedure for the solution of multi-component reactive transport problems. Water Resources Research 41, W11410. doi:10.1029/2005WR004056. De Simoni, M., Sanchez-Vila, X., Carrera, J., Saaltink, M.W., 2007. A mixing ratios-based formulation for multicomponent reactive transport. Water Resources Research 43, W07419. doi:10.1029/2006WR005256. Dijk, P., Berkowitz, B., 1998. Precipitation and dissolution of reactive solutes in fractures. Water Resources Research 34 (3), 457–470. Domenico, P.A., Schwartz, F.W., 1998. Physical and Chemical HydrogeologySecond edition. John Wiley, New York. Dullien, F.A.L., 1979. Porous Media: Fluid Transport and Pore Structure. Academic Press, New York. Emmanuel, S., Berkowitz, B., 2005. Mixing-induced precipitation and porosity evolution in porous media. Advances in Water Resources 28 (4), 337–344. Fogler, H.S., Rege, S.D., 1986. Porous dissolution reactors. Chemical Engineering Communications 42 (4–6), 291–313. Gramling, C.M., Harvey, C.F., Meigs, L.C., 2002. Reactive transport in porous media: a comparison of model prediction with laboratory visualization. Environmental Science and Technology 36 (11), 2508–2514. Guadagnini, A., Sanchez-Vila, X., Saaltink, M.W., Bussini, M., Berkowitz, B., 2009. Application of a mixing-ratios based formulation to model mixingdriven dissolution experiments. Advances in Water Resources 32 (5), 756–766. doi:10.1016/j.advwatres.2008.07.005. Helgerson, H.C., Kirkham, D.H., 1974. Theoretical prediction of the thermodynamic behaviour of aqueous electrolytes at high pressure and temperature: II Debye–Hückel parameters for activity coefficients and relative partial molal properties. American Journal of Science 274, 1199–1261. Higgins, C.G., 1980. Nips, notches, and the solution of coastal limestone — overview of the problem with examples from Greece. Estuarine Coastal Marine Science 10 (1), 15–30. Hoefner, M.L., Fogler, H.S., 1988. Pore Evolution and channel formation during flow and reaction in porous-media. AIChE American Institute of Chemical Engineers Journal 34 (1), 45–54. Kelemen, P.B., Whitehead, J.A., Aharonov, E., Jordahl, K.A., 1995. Experiments on flow focusing in soluble porous-media, with applications to melt extraction from the mantle. Journal of Geophysical Research 100 (B1), 475–496. Kinzelbach, W., 1986. Groundwater modelling: an introduction with sample programs in BASIC. Elsevier, New York. 343 pp. Kinzelbach, W., Schafer, W., Herzer, J., 1991. Numerical modeling of natural and enhanced denitrification processes in aquifers. Water Resources Research 27 (6), 1123–1135. Lasaga, A.C., 1984. Chemical-kinetics of water–rock interactions. Journal of Geophysical Research 89 (Nb6), 4009–4025. Levy, M., Berkowitz, B., 2003. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. Jounal of Contaminant Hydrology 64 (3–4), 203–226. Liu, C.W., Narasimhan, T.N., 1989. Redox-controlled multiple-species reactive chemical-transport.1. Model development. Water Resources Research 25 (5), 869–882. Mangold, D.C., Tsang, C.F., 1991. A summary of subsurface hydrological and hydrochemical models. Reviews of Geophysics 29 (1), 51–79. Molins, S., Carrera, J., Ayora, C., Saaltink, M.W., 2004. A formulation for decoupling components in reactive transport problems. Water Resources Research 40, W10301. doi:10.1029/2003WR002970. MoMaS (Modeling, Mathematics and numerical Simulations related to nuclear waste management problems) modeling initiative (http:// www.gdrmomas.org/benchmarks-en.html), 2008. Novak, C.F., 1993. Modeling mineral dissolution and precipitation in dual-porosity fracture matrix systems. Journal of Contaminant Hydrology 13 (1–4), 91–115. Novak, C.F., Schechter, R.S., Lake, L.W., 1989. Diffusion and solid dissolution– precipitation in permeable media. American Institute of Chemical Engineering Journal 35 (7), 1057–1072.

44

G.E. Katz et al. / Journal of Contaminant Hydrology 120–121 (2011) 27–44

Parkhurst, D.L., Appelo, C.A.J., 1999. User's guide to PHREEQC (version 2) — a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. U.S. Geological Survey Water-Resources Investigations Report 99-4259, p. 310. Rezaei, M., Sanz, E., Raeisi, E., Ayora, C., Vázquez-Suñé, E., Carrera, J., 2005. Reactive transport modeling of calcite dissolution in the fresh-salt water mixing zone. Journal of Hydrology 311 (1–4), 282–298. Romanov, D., Dreybrodt, W., 2006. Evolution of porosity in the saltwater– freshwater mixing zone of coastal carbonate aquifers: an alternative modelling approach. Journal of Hydrology 329 (3–4), 661–673. Saaltink, M.W., Ayora, C., Carrera, J., 1998. A mathematical formulation for reactive transport that eliminates mineral concentrations. Water Resources Research 34 (7), 1649–1656. Saaltink, M.W., Batlle, F., Ayora, C., Carrera, J., Olivella, S., 2004. RETRASO, a code for modeling reactive transport in saturated and unsaturated porous media. Geologica Acta 2, 235–251. 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. Singurindy, O., Berkowitz, B., 2003. Evolution of hydraulic conductivity by precipitation and dissolution in carbonate rock. Water Resources Research 39 (1), 1016. doi:10.1029/2001WR001055, 2003. Singurindy, O., Berkowitz, B., 2004. Dedolomitization and flow in fractures. Geophysical Research Letters 31, L24501. Singurindy, O., Berkowitz, B., Lowell, R.P., 2004. Carbonate dissolution and precipitation in coastal environments: laboratory analysis and theoretical consideration. Water Resources Research 40, W04401. doi:10.1029/ 2003WR002651.

Smart, P.L., Dawans, J.M., Whitaker, F., 1988. Carbonate dissolution in a modern mixing zone. Nature 335 (6193), 811–813. Steefel, C.I., MacQuarrie, K.T.B., 1996. Approaches to modeling of reactive transport in porous media. In: Lichtner, P.C., Steefel, C.I., Oelkers, E.H. (Eds.), Reactive Transport in Porous Media: Reviews in Mineralogy, vol. 34. Mineralogical Society of America, Washington DC, pp. 83–125. 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. Tartakovsky, A.M., Redden, G., Lichtner, P.C., Scheibe, T.D., Meakin, P., 2008. Mixing-induced precipitation: experimental study and multiscale numerical analysis. Water Resources Research 44, W06S04. doi:10.1029/ 2006WR005725. Walter, A.L., Frind, E.O., Blowes, D.W., Ptacek, C.J., Molson, J.W., 1994. Modeling of multicomponent reactive transport in groundwater 1. Model development and evaluation. Water Resources Research. 30 (11), 3137–3148. Wunderly, M.D., Blowes, D.W., Frind, E.O., Ptacek, C.J., 1996. Sulfide mineral oxidation and subsequent reactive transport of oxidation products in mine tailings impoundments: a numerical model. Water Resources Research 32 (10), 3173–3187. Yeh, G.T., Tripathi, V.S., 1991. Model for simulating transport of reactive multispecies components — model development and demonstration. Water Resources Research 27 (12), 3075–3094.