Coupled thermohaline groundwater flow and single-species reactive solute transport in fractured porous media

Advances in Water Resources 30 (2007) 742–771 www.elsevier.com/locate/advwatres

Coupled thermohaline groundwater flow and single-species reactive solute transport in fractured porous media Thomas Graf *, Rene´ Therrien De´partement de Ge´ologie et Ge´nie Ge´ologique, Universite´ Laval, Ste-Foy, Que., Canada G1K 7P4 Received 15 December 2005; received in revised form 26 June 2006; accepted 6 July 2006 Available online 17 August 2006

Abstract A 3D numerical model has been developed to solve coupled fluid flow, heat and single-species reactive mass transport with variable fluid density and viscosity. We focus on a single reaction between quartz and its aqueous form silica. The fluid density and viscosity and the dissolution rate constant, equilibrium constant and activity coefficient are calculated as a function of the concentrations of major ions and temperature. Reaction and flow parameters, such as mineral surface area and permeability, are updated at the end of each time step with explicitly calculated reaction rates. Adaptive time stepping is used to increase or decrease the time step size according to the rate of temporal variation of the solution to prevent physically unrealistic results. The time step size depends on maximum changes in matrix porosity and/or fracture aperture. The model is verified against existing analytical solutions of heat transfer and reactive transport in fractured porous media. The complexity of the model formulation allows studying chemical reactions and variable-density flow in a more realistic way than done previously. The newly developed model has been used to simulate illustrative examples of coupled thermohaline flow and reactive transport in fractured porous media. Simulations indicate that thermohaline (double-diffusive) transport impacts both buoyancy-driven flow and chemical reactions. Hot zones correspond to upwelling and to quartz dissolution while in cooler zones, the plume sinks and silica precipitates. The silica concentration is inversely proportional to salinity in high-salinity regions and proportional to temperature in lowsalinity regions. Density contrasts are generally small and fractures do not act like preferential pathways but contribute to transverse dispersion of the plume. Results of a long-term (100 years) simulation indicate that the coexistence of dissolution and precipitation leads to self-sealing of fractures. Salt mass fluxes through fractures decrease significantly due to major fracture aperture reduction in the precipitation zone. The system is the most sensitive to temperature because it impacts both the dissolution kinetics (Arrhenius equation) and the quartz solubility. The system is least sensitive to quartz surface area in the fracture because the volumetric fraction of a fracture is small compared to the volumetric fraction of the porous matrix.  2006 Elsevier Ltd. All rights reserved. Keywords: Numerical modeling; Nuclear waste; Fracture; Quartz; Reactive transport; Density; Heat; Thermohaline

1. Introduction 1.1. Problem definition Many countries generate electrical power using nuclear fuel. Altogether, 436 nuclear power plants around the world operate in 31 countries. Nuclear energy production *

Corresponding author. Tel.: +1 418 6562131; fax: +1 418 6567339. E-mail address: [email protected] (T. Graf).

0309-1708/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2006.07.001

eventually creates waste in the form of spent nuclear fuel. Deep burial of radioactive waste in low-permeability geological formations is a disposal option considered worldwide [11]. Studies show that, with increasing depth, groundwater and ambient rock temperature increase according to geothermal gradient and water composition can reach that of hot saline Na–Ca–Cl brines [17,71]. Because of these temperature and composition changes, deep-fluid properties such as viscosity and density can no longer be assumed constant. For deep waste disposal,

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

743

Nomenclature The use of symbols for main variables is consistent throughout the entire text All symbols represent scalar variables denoted in normal italic letters Latin letters (2b) [L] fracture aperture Aqz [M1 L2] specific surface area in the matrix 1 Afr L2  specific surface area in the fracture qz ½M 2 As [L ] active surface area B [MOL1 L3] coefficient in the Jones–Dole equation ~c ½L2 T2 #1  specific heat C [M L3] solute concentration, expressed as volumetric mass D [–] Marshall–Chen coefficient Dd [L2 T1] free-solution diffusion coefficient Dij [L2 T1] hydrodynamic matrix dispersion tensor 2 1 Dfr ij ½L T  hydrodynamic fracture dispersion coefficient Ea [MOL1 M L2 T2] activation energy g [L T2] acceleration due to gravity h0 [L] equivalent freshwater head I+, I [–] fracture–matrix interface k [M L T3 #1] thermal conductivity k 0þ ½MOL L2 T1  dissolution reaction constant in deionized water k corr ½MOL L2 T1  dissolution reaction constant in þ saltwater Kad [MOL1 M] equilibrium adsorption coefficient Kd [M1 L3] equilibrium distribution coefficient K fr d ½L fracture–surface distribution coefficient Keq [MOL M1] equilibrium constant K 0ij ½LT1  coefficients of hydraulic conductivity tensor of freshwater 1 K fr ½LT  hydraulic freshwater conductivity of the 0 fracture ‘v [L] geometry of the model domain; v = x, y, z Lv [L] geometry of a block element; v = x, y, z m [MOL M1] molal concentration M [MOL L3] molar concentration Mw [M] mass of water P [M L1 T2] dynamic pressure of the fluid PPM [–] mass parts per million qi [L T1] Darcy flux r [MOL M1 T1] precipitation (backward) reaction rate r+ [MOL M1 T1] dissolution (forward) reaction rate rnet [MOL M1 T1] net molal production rate rM [MOL L3 T1] net molar production rate R [–] retardation factor R* [MOL1 M L2 T2 #1] universal gas constant R = 8.3144 mol1 kg m2 s2 K1 fr R [–] fracture retardation factor SS [L1] specific storage of the porous matrix 1 S fr S ½L  specific storage of an open fracture

t [T] time T [#] absolute temperature in Kelvin TC [#] relative temperature in centigrade vi [L T1] linear flow velocity Vqz [MOL1 L3] molar volume of quartz Greek letters afl [M1 L T2] coefficient of the compressibility of the fluid due to fluid pressure or hydraulic head variations al [L] matrix longitudinal dispersivity afr l ½L longitudinal fracture dispersivity am [M1 L T2] coefficient of the compressibility of the porous medium due to fluid pressure or hydraulic head variations at [L] matrix transverse dispersivity afr t ½L transverse fracture dispersivity asalt [M1 L3] solutal expansion coefficient cr [–] activity coefficient of species r Cm variable mass sources and sinks dij [–] Kronecker delta function gj [–] indicator for flow direction hr [–] fraction of sites occupied by cation r jij [L2] coefficients of the intrinsic permeability tensor k [T1] decay constant K [M T3] convective–dispersive–conductive loss or gain of heat l [M L1 T1] dynamic viscosity of the fluid q [M L3] density qr [–] relative fluid density q~c ½M L1 T2 #1  heat capacity s [–] factor of tortuosity / [–] porosity of the rock matrix /qz [–] quartz volume fraction u [1] fracture incline x [–] fracture roughness coefficient X [M M1 T1] advective–dispersive–diffusive loss or gain of solute mass Sub- and superscripts 0 [–] reference fluid b [–] bulk fr [–] fracture i, j [–] spatial indices init [–] initial time level l [–] liquid phase L [–] time level n [–] normal direction s [–] solid phase r [–] species Special symbols o [–] partial differential operator n [–] difference

744

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

critical safety questions arise due to the presence of fractures in low-permeability rock formations. Fractures have a great impact on mass transport because they represent preferential pathways where accidentally released contaminants migrate at velocities that are several orders of magnitude larger than within the rock matrix itself. Significant increases in temperature can promote rockfluid interactions such as mineral dissolution and precipitation. Physical properties, such as matrix permeability and fracture aperture are modified if chemical rock-fluid interactions occur. These modifications in physical properties can be significant because the Cubic Law states that, for example, an increase of the fracture aperture by 26% doubles the discharge through this fracture. The high number of conceivable feedback scenarios between variable-density flow and reactive solute transport demonstrates that the two processes are strongly coupled. This is especially the case in fractured media where high groundwater flow velocities enable rapid transport of reactive species to the location of the chemical reaction and away from it. Clearly, the ability to predict the transport behavior of hazardous chemicals leaked to the geosphere is essential. Processes that affect the fluid properties on one hand and induce reactions on the other hand are elevated temperatures and high salt contents. Numerical models are very helpful tools for studying the behavior and predicting long-term effects in such complex systems. Models that couple variable-density flow with reactive transport are relatively new and ‘‘the development of these codes has only just begun’’ [50]. Already available models vary greatly in their coupling method and in model sophistication [50]. 1.2. Prior studies Reactive transport models are typically limited to the chemical system being investigated. Several modelling studies considered multiple mobile species undergoing reactions [25,56,85,65,81,88,58,6,68,67,59,26,51,28,55,27,21,41,22]. Some models consider reactions with aqueous silica (H4SiO4) as the only mobile reactive component [34,66,83]. In this case, the reactive transport equation remains linear and iterative solvers for nonlinear equations, such as those described by Steefel and MacQuarrie [68], do not need to be applied. Several authors examined reactive solute transport in porous media assuming constant water density (e.g. [68,69]). Bolton et al. [6] and Freedman and Ibaraki [21] investigated the impact of density-driven flow on chemical reactions. Bolton et al. [6] studied coupled thermal convection and quartz dissolution/precipitation for large spatial and temporal scales. They found that long-term changes of porosity and permeability can either increase the flow velocities and the degree of subsaturation (in regions of dissolution) or decrease flow rates and the degree of supersaturation (in regions of precipitation). However, Bolton et al. [6] did not account for the salinity dependency of the kinetic rate law nor for the salinity effect on water density

and viscosity. Freedman and Ibaraki [21] numerically simulated the horizontal migration of a dense plume in an unfractured porous medium where density varies with salinity but not with temperature. The results were compared with simulations where chemical reactions are ignored. The most important outcome of the Freedman and Ibaraki [21] study was that chemical reactions do not significantly impact density-driven flow in porous media. However, this finding may not be always valid because the type of chemical system studied will influence the outcome. In Freedman and Ibaraki [21], the solubility of the solid phase (calcite) is fairly low and coupling between flow and reactive transport is therefore weak. In addition, [21] focused on a small temporal scale and did not study long-term effects. They also ignored the influence of temperature and salinity on both solute solubility and reaction kinetics. Simulations of reactive transport in fractured systems have previously been carried out by a number of authors [65,67,28,27]. Not all of the studies addressed the question of how dissolution/precipitation reactions will alter fracture aperture and matrix permeability and, thus, impact the flow field. Modifications of flow parameters were either not considered [28] or only applied to the permeability of the porous matrix [27]. However, other investigations showed that chemical reactions within open fractures trigger complex reaction-flow feedback scenarios [65,67] and that fracture aperture cannot be assumed constant. Steefel and Lichtner [67] studied the infiltration of a hyper-alkaline fluid along a discrete fracture. They observed that within tens of years, the permeability feedback between reaction and transport is significant. They also found that fluid flow through the fracture is likely to be restricted, or even stopped, due to self-sealing if reaction rates in the fracture are only one order of magnitude larger than in the adjacent matrix. On the other hand, if the rates are of the same order of magnitude, the porous matrix will be cemented first. While most studies ignored density variations when investigating reactive transport in fractured media [67,28,27,65,49] fully accounted for thermal density-driven flow. According to Steefel and Lasaga [65], geothermal convection cells in reactive fractured media are never stable because upwelling fluids cool and the resulting precipitation of minerals significantly reduces permeability leading to highly dispersive plumes. On the other hand, if fluids move downward to a zone of higher temperature, dissolution reactions locally increase permeability leading to channelling of flow. However, Steefel and Lasaga [65] did not account for the impact of salinity on both reaction kinetics and fluid properties. Being designed for simulating flow in salt domes, the TOUGH2 model accounts for thermohaline impact on fluid density and viscosity [51]. Oldenburg and Pruess [49] used this model to simulate thermohaline convection in porous media. In addition, TOUGH2 can simulate flow and transport in fractured media using the dual permeabil-

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

ity approach. TOUGH2 also simulates reactive transport of salt (NaCl), including temperature-dependent solubility and dissolution/precipitation of salt. However, Pruess et al. [51] did not account for the impact of salinity on reaction kinetics. Nevertheless, other studies suggested that the rate-enhancing effect of salt is significant [13] and that fluid salinity also impacts the reactive species solubility [40,36]. Table 1 shows a selection of previous studies on coupled variable-density flow and reactive transport. It highlights subtle differences between model assumptions made by various authors. The studies are not listed chronologically but according to increasing model complexity. This paper presents the development and application of a numerical model that simulates dense plume migration in a chemically reactive fractured geologic material, for nonisothermal conditions. This new model is based on the existing FRAC3DVS model, which solves variable-saturated and multi-component transport in discretely fractured porous media [76]. The model can be used to simulate the long-term behavior of coupled thermohaline flow and reactive solute transport in fractured media. The model developed in the present study continues the series of increasing model complexity and provides simulation capabilities previously lacking (Table 1).

745

2. Physicochemical system 2.1. Chemical system We have chosen the quartz-water chemical system because silicate minerals are the most abundant minerals in the earth’s crust, accounting for 90% of the total mass [35]. The focus will be on a-quartz, the most common SiO2 polymorph in the upper crust. Several studies underlined the importance to examine the reactive nature of quartz-rich rock. Fournier [18,20] demonstrated that quartz precipitation may significantly decrease permeability and thus lead to self-sealing of the rock matrix. The essential role of quartz dissolution has been studied since 1884 [9] and the high number of more modern studies [7,14,75,15,83,13] indicates that quartz dissolution is also a recent subject of intensive research. It is assumed that aqueous silica is the only species that undergoes chemical interactions (dissolution/precipitation) with the rock matrix and it is thus termed ‘‘reactive species’’. Other species that undergo sorption reactions but not dissolution/precipitation will not be termed ‘‘reactive’’. Therefore, the model presented here is a single-species reactive transport model.

Table 1 A selection of previous studies on reactive solute transport in porous and fractured porous media Authors

Steefel and MacQuarrie [68] Johnson et al. [34] Steefel and Yabusaki [69] Freedman and Ibaraki [21] White and Mroczek [83] Bolton et al. [6] Ghogomu and Therrien [28] Geiger et al. [27] Steefel and Lichtner [67] Steefel and Lasaga [65] Pruess et al. [51] Present study

Simulated processes Reactive transport

Heat transfer

Density from

in PMa pe p /,A p /,A p /,j,A p /,j,A p /,j,A pe pj p /,A p /,A p /,j p /,j,A

in FMb

in PM

in FM

Cc

– – – – – – pe pe p (2b),A p (2b),A p /,j,DP p (2b),A

– – – – – p

– – – – – – – – – p p p

– – – p

– – – p p p

Viscosity from

Reaction kinetics from

Reactive species solubility from

Td

C

T

C

T

– – – p

– – – – – p

– – – – p

– p p

– p

– – – – – p

C pf

– – – – p p

– – – p p p

– – – – – – p p

– – – p p p

– – – – – – p

– p p – – – p p p

–qz pf pf pqz pf pf pf pf pf pf pqz

T – p p – p p – – – p p p

The models are listed in order of increasing complexity. If density is a function of salinity and/or temperature, the model couples reactions with variabledensity flow. DP Dual permeability approach used. a Porous media. b Fractured media. c Salinity. d Temperature. e No change of simulation parameters considered. f Multi-species reactive transport. qz Single-species reactive transport (silica) and nonreactive transport (electrolytes). (2b) Change of fracture aperture considered. j Change of matrix permeability considered. / Change of matrix porosity considered. A Change of specific mineral surface area considered.

746

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

Rimstidt and Barnes [53] experimentally studied the quartz-water system, described by the reaction rþ

SiO2ðsÞ þ 2H2 OðaqÞ  H4 SiO4ðaqÞ r

ð1Þ

where r+ and r [both MOL M1 T1] are the dissolution and precipitation rates of quartz (SiO2), respectively, and where H4SiO4 is aqueous silica. Upon applying the law of mass action to reaction (1), the net rate of silica production in the porous matrix, rnet, can be written as [53,15,13]   cH4 SiO4 corr rnet ¼ /qz k þ Aqz 1  mH4 SiO4 ð2Þ K eq where /qz [–] is the volume fraction of quartz, k corr ½MOL L2 T1  is the dissolution rate constant, corþ rected for salt water, Aqz [M1 L2] is the specific quartz surface area, cH4 SiO4 ½– is the activity coefficient of silica, Keq [MOL M1] is the quartz solubility or equilibrium constant of reaction (1) and mH4 SiO4 ½MOL M1  is the molal concentration of silica. It is assumed here that pure solids and pure liquids, for example SiO2 and H2O, have activities equal to unity [35]. For a discrete fracture, the net rate of silica production is given by   cH4 SiO4 fr corr fr fr rnet ¼ /qz k þ Aqz 1  mH4 SiO4 ð3Þ K eq where superscript ‘‘fr’’ refers to fracture and where other variables are defined similarly to those used for the porous medium. This chemical model is based on transition state theory and simulates a zeroth/first order chemical kinetic reaction. The model is similar to that used by Rimstidt and Barnes [53] and White and Mroczek [83], except that the rate law in (2) and (3) is also a function of the quartz volume fraction, /qz, as proposed by Johnson et al. [34]. The following paragraphs present in more detail parameters k corr þ , cH4 SiO4 and Keq. 2.1.1. Corrected dissolution rate constant, k corr þ The corrected dissolution rate constant, k corr þ , is obtained from the dissolution rate constant in deionized water, k 0þ , which is given by the Arrhenius equation as (e.g. [37,69])    Ea 1 1  k 0þ ¼ k 025 exp ð4Þ R T 298:15 k 025

2

assumed here as proposed by Rimstidt and Barnes [53] and used in Steefel and Lasaga [65] as well as in the software packages OS3D and GIMRT [69]. Eq. (4) holds for deionized water but [14] have shown that the presence of electrolytes in the fluid can increase the reaction rate by 1.5 orders of magnitude. The increase is caused by adsorbed cations that change the structure of the mineral surface and make the surface more vulnerable to water dipole attacks. According to Dove and Nix [15], the concentration of the bivalent (IIA) cations (Mg2+, Ca2+) have the greatest impact on dissolution rate while the effect of monovalent (IA) cations (Na+, K+) is minor due to their less effective adsorption. Dove [13] demonstrated that the fraction of adsorption sites occupied by species Na+, Mg2+ and Ca2+ can be expressed by a Langmuir model for equilibrium adsorption [54,5,14] as þ

þ K Na ad mNa

hNa ¼ þ

2þ

ð5Þ

2þ

K Mg mMg2þ ad

hMg2þ ¼

2þ

þ

2þ

Mg þ 1 þ K Na mMg2þ þ K Ca ad mNa þ K ad ad mCa2þ

ð6Þ

2þ

K Ca ad mCa2þ

hCa2þ ¼

2þ

þ

2þ

Mg þ 1 þ K Na mMg2þ þ K Ca ad mNa þ K ad ad mCa2þ

ð7Þ

where hr [–] is the fraction of sites occupied by cation r, mr [MOL M1] and K rad ½MOL1 M are the molal concentration and the equilibrium adsorption coefficient of cation r, respectively. Because sorption constants, K rad , for cations on quartz at hydrothermal temperatures are unknown [15], we used K rad values at 20 C. These values are equal to 101.78, 103.7 and 103.35 mol1 kg for sodium, magnesium and calcium, respectively [13]. The constant k rþ is computed from a fit to experimental data published by Dove [13] who measured the dependence of quartz dissolution rates on different electrolyte concentrations of sodium chloride, magnesium chloride and calcium chloride at 200 C. The logarithm of the dissolution rate constant for Na+, Mg2+ and Ca2+ can be written as þ

log k Na 200 ¼  2þ

log k Mg 200

2:8  104  6:35 mNaþ

¼

1

½MOL L T  is the known dissolution rate where constant in deionized water at 25 C, Ea [MOL1 M L2 T2] is the activation energy necessary to overcome the potential energy maximum of the transition state and T [#] is absolute temperature. Values of the universal gas constant, R* [MOL1 M L2 T2 #1], and constant k 025 are given in Rimstidt and Barnes [53] as 8.3144 mol1 kg m2 s2 K1 and 4.3 · 1014 mol m2 s1, respectively. Values of Ea for quartz dissolution have been reported in the range between 36 and 96 kJ mol1, and a value 75.0 kJ mol1 is

2þ

þ

Mg þ 1 þ K Na mMg2þ þ K Ca ad mNa þ K ad ad mCa2þ

2þ

log k Ca 200 ¼ 

2:2  104  6:80 mMg2þ

1:3  106 ðmCa2þ Þ

2

 6:35

ð8Þ ð9Þ ð10Þ

where log denotes the decadic logarithm log10. With the help of the Arrhenius equation (4) and with Eqs. (8)– (10), the dissolution rate constant of species r at any concentration and temperature can be formulated as    Ea 1 1  k rþ ¼ k r200 exp ð11Þ R T 473:15

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

The corrected dissolution rate constant, k corr þ , accounts for the presence of Na+, Mg2+ and Ca2+ in the solution and is defined here as þ

2þ

2þ

Na Mg þ k corr hMg2þ þ k Ca hCa2þ þ ¼ k þ hNa þ k þ þ h  i þ k 0þ 1  hNaþ þ hMg2þ þ hCa2þ

ð12Þ

where Dove’s [13] idea of competitive adsorption is adapted in order to account for protons ‘‘adsorbed’’ on the remaining sites, expressed by the last term in Eq. (12). This last term does not occur in Dove’s [13] original formulation of the dissolution rate constant in a mixed electrolyte solution but it is necessary to obtain a correct rate constant in water of very low salinity. In that case, the values of mr and hr approach zero and k corr becomes þ approximately equal to k 0þ . Note that competitive adsorption of cations is only used in a qualitative sense to assess the effect of salinity on dissolution rates. This is because the model does not explicitly simulate competitive adsorption but sorption of cations is simulated with a linear isotherm (Eqs. (39) and (42) in Section 3.2). Electrolyte concentration and fluid temperature are the main factors that affect quartz dissolution rates. Between pH 8 and pH 12, quartz dissolution is also a function of water acidity [7]. However, the expected range of pH values of water in the quartz-rich rock considered here is below 7. This pH is controlled by silica dissolution and subsequent buffering by the silicic acid buffer [47]. For pH values below 7, changes of dissolution rates are ‘‘small [. . .] and difficult to interpret’’ [3]. Therefore, the pH dependency of quartz dissolution rates is neglected in this study. 2.1.2. Activity coefficient, cH4 SiO4 In an electrolyte solution, the solubility of a neutral species, such as H4SiO4, is a function of the amount of dissolved salt and temperature. Marshall and Chen [40] proposed a modified form of the Setchenow equation to calculate the activity coefficient of H4SiO4 in a mixed electrolyte solution at any given temperature: X Dr mr ð13Þ log cH4 SiO4 ¼ r

where Dr is the dimensionless, temperature dependent Marshall–Chen coefficient of ion r and mr is the molal concentration [MOL M1] of r. Marshall and Chen [40] provide temperature-dependent values of Dr for species Na+, Mg2+, Cl and SO2 4 in the range 25–300 C. Due to the physicochemical similarity of Mg2+ and Ca2+, their Marshall–Chen coefficients are assumed to be identical, such that DCa2þ ¼ DMg2þ . The further assumption is made that Dr can be extrapolated beyond the 25–300 C temperature range down to 0 C. 2.1.3. Equilibrium constant, Keq In this model, the equilibrium constant, Keq, is expressed as a function of the absolute temperature, T, over the range 0–300 C [52] such that

log K eq ¼ 

1107:12  0:0254 T

747

ð14Þ

Other models use an apparent solubility, which also accounts for the impact of salinity [19,80,62,46]. However, the present model takes ion activity into consideration by calculating nonzero Marshall–Chen coefficients in Eq. (13), thus using silica activity coefficients that are greater than or equal to one. In the calculation of the quartz solubility, it is further assumed that below pH 9, there is no influence of pH on quartz solubility [35,72] and that the pressure effect on quartz solubility is not significant within the temperature range considered in the scope of this study [83]. 2.2. Physical system In thermohaline convective systems, the fluid density and viscosity depend on temperature and salinity, while the effect of pressure on fluid properties can be ignored [6]. In the model, both fluid quantities  are first calculated as a function of temperature alone q0T ; l0T and then corrected according to salinity. 2.2.1. Fluid density, q Under isobaric conditions, the fluid density is calculated as a function of temperature for different temperature ranges [32]: 8 2 > 1000  ð1  ð½T C  3:98 =503; 570Þ > > > > > > > ð½T C þ 283=½T C þ 67:26ÞÞ > > > for 0  C < T C 6 20  C > > > > 4 > > < 996:9  ð1  3:17  10 ½T  298:15 q0T ¼ 2:56  106 ½T  298:152 Þ > > > > for 20  C < T C 6 175  C > > > > > 1758:4 þ 1000  T ð4:8434  103 > > > > > þT ð1:0907  105  T  9:8467  109 ÞÞ > > > : for 175  C < T C 6 300  C ð15Þ where TC [#] and T [#] are the temperatures in centigrade and Kelvin, respectively. In a second step, the fluid density at any given salinity and temperature is evaluated using the freshwater density at temperature, q0T , and from the sum of all species concentrations using the following empirical relation: X q ¼ q0T þ asalt  Cr ð16Þ r

salt

1

[M L3] is the solutal expansion coefficient. where a Due to the low solubility of quartz, the impact of dissolved silica on fluid density is not significant [19,46] and, therefore, ignored in Eq. (16). The model calculates density from the concentration of eight major ions found in natural 2 waters: Na+, K+, Ca2+, Mg2+, Cl, SO2 and 4 , CO3  HCO3 . The Pitzer ion interaction model is used [44,45],

748

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

where the fluid density is derived from the partial electrolyte volumes. The Monnin model is used to derive an empirical expression for asalt as a function of the groundwater chemistry in the Canadian Shield given by Farvolden et al. [17] in the form ! X salt a ¼ 0:0829  ln C r þ 1:1415 ð17Þ r

where the unit of Cr [M L3] is mg l1. 2.2.2. Fluid viscosity, l Although fluid viscosity is assumed constant in some thermohaline transport studies [70,78,79,16,8,86], it is recommended to relate viscosity to both temperature [6] and salinity [21] because it can increase by a factor of two between pure water and a dense brine [49]. Relations to calculate fluid viscosity for different temperature ranges are 8 1:787  103  expðð0:03288 þ 1:962  104  T C Þ  T C Þ > > > > > for 0  C < T C 6 40  C > > > 3 < 10  ð1 þ 0:015512  ½T C  20Þ1:572 l0T ¼ > for 40  C < T C 6 100  C > > > > ^ > > 0:2414  10 ð247:8=½T C þ 133:15Þ > : for 100  C < T C 6 300  C ð18Þ

The relation for temperatures below 40 C has been used by Molson et al. [43] and the other two relations are given by Holzbecher [32]. Viscosity can be expressed as a function of salinity and temperature by substituting the temperature-dependent freshwater viscosity, l0T , in the Jones–Dole equation: ! X 0 l ¼ lT  1 þ Br M r ð19Þ r

where Mr is the molar concentration of species r. Marcus [39] gives values of the B-coefficients [L3 MOL1] for all species considered here. 2.3. Solid phase properties Chemical reactions have a significant impact on a number of physical flow and transport properties. In the model, the individual quartz volume fraction, /qz [–], is recalculated using [69,65]: o/qz ¼ V qz rLþ1 M ot

ð20Þ

where Vqz [MOL1 L3] is the molar volume of quartz and 3 rLþ1 T1  is the molar reaction rate at time level M ½MOL L L+1. In finite difference form, this equation becomes L Lþ1 /Lþ1 qz ¼ /qz  DtV qz r M

ð21Þ

where Dt = tL+1  tL [T] is the time step size. The molar volume of a mineral is the ratio of its molecular weight

to its density [36]. The model updates porosity from the sum of all mineral volume fractions: X Lþ1 /r ð22Þ /Lþ1 ¼ 1  r

where it is assumed that quartz is the only reactive solid species and that /r is constant for any given mineral r except quartz. The specific surface area in the porous matrix is updated by means of the two-thirds power relation given by Steefel and Yabusaki [69] as 8h  i2=3 > Lþ1 init Lþ1 init > ð/ =/ Þ  / =/ > qz qz > > < Lþ1 init dissolution of quartz ð23Þ Aqz ¼ Aqz > Lþ1 init 2=3 > ð/ =/ Þ > > > : precipitation of quartz 2 1 where Ainit qz ½L M  is the initial specific surface area in the matrix and where /init [–] and /init qz ½– are the initial matrix porosity and quartz fraction, respectively. The matrix permeability, jij [L2], is calculated from porosity for the special case of dissolution and precipitation of quartz [82]: 8 2 !1:58 30:46 9 < = Lþ1 init 41  /  /c 5 ð24Þ ¼ j  1  jLþ1 ij ij : ; /init  /c 2 where jinit ij ½L  is the initial permeability and /c [–] is the critical porosity at which jij = 0. This relation is obtained from theoretical considerations of deposition and dissolution of quartz grains, arranged in a rhombohedral array of uniform spheres [82]. Similar to the matrix porosity, fracture apertures are recalculated from [66,67]: Lþ1 L  ð2bÞ ð25Þ ¼ ð2bÞ  1 þ DtV qz rLþ1 M

Finally, the specific surface area in the fracture is updated using [66]: ! Lþ1 ð2bÞ fr;Lþ1 fr;init Aqz ¼ Aqz  ð26Þ init ð2bÞ where Afr;init ½L2 M1  is the initial specific surface area in qz the fracture and where (2b)init [L] is the initial fracture aperture. The initial surface area in a 2D rectangular fracture element is the ratio between active surface area, given by As = 2xLxLz, and the mass of water stored in the 2D element, Mw = q Æ (2b)LxLz, where x [–] is the fracture roughness coefficient (Fig. 1). Thus, a fluid moving through a large fracture will encounter less mineral surface area per unit fluid mass than a fluid moving through a narrow fracture. That relationship is expressed by the following equation for the initial specific surface area: Afr;init ¼ qz

x qb

ð27Þ

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

749

2

K fr0 ¼

Fig. 1. Fracture roughness coefficient for rough-walled (left) and smooth fractures.

where (2b) [L] is the fracture aperture. The application of Darcy’s law in fractures (31) requires that the Reynolds number be smaller than 1 [2]. The equation that governs variable-density, variable-viscosity flow in porous media has the following 3D form [24,30]:    o oh0 0 l0 oh0 K þ q r gj ¼ SS ; i; j ¼ 1; 2; 3 ð33Þ oxi ij l oxj ot The specific storage, SS [L1], accounts for both matrix and fluid compressibility and is defined as [1] ð34Þ 1

3.1. Variable-density, variable-viscosity flow The model uses the equivalent freshwater head h0 [L], defined by Frind [24] as P þz q0 g

ð32Þ

S S ¼ q0 gðam þ /afl Þ

3. Governing equations

h0 ¼

ð2bÞ q0 g 12l0

ð28Þ

where P [M L1 T2] is the dynamic fluid pressure, q0 [M L3] is the reference fluid density, g [L T2] is the gravitational acceleration and z [L] is the elevation above datum. The 3D variable-density, variable-viscosity Darcy flux in porous media can be completely expressed in terms of freshwater properties [24,30]:   0 l0 oh0 qi ¼ K ij þ qr gj ; i; j ¼ 1; 2; 3 ð29Þ l oxj where the assumption of a horizontal datum (i.e., oz/ oz = 1) is made and where gj [–] represents the direction of flow with gj = 0 in the horizontal directions and gj = 1 in the vertical direction [24]. The relative density, qr [–] is given by qr = (q/ q0)  1. Using the ratio between the reference fluid viscosity, l0, and the fluid viscosity, l [both M L1 T1], the hydraulic conductivity of the porous medium for freshwater, K 0ij ½L T1 , is [1] jij q0 g K 0ij ¼ ð30Þ l0 The 2D Darcy flux in differently oriented 2D fracture faces is calculated using a form of the Darcy equation similar to that presented by Graf and Therrien [30]:  fr  oh0 fr fr l0 fr qi ¼ K 0 fr þ qr gj cos u ; i; j ¼ 1; 2 ð31Þ l oxj where gj is 0 in the horizontal direction and 1 along the fracture incline. The incline of a fracture face is given by u with u = 0 for a vertical face and u = 90 for a horizontal face. In the case of flow within fractures, a local 2D Cartesian coordinate system is assumed. The freshwater hydraulic conductivity of the fracture, K fr0 ½L T1 , is derived from the parallel plate model as

2

where am and afl [both M L T ] are the matrix and fluid compressibility, respectively. Flow in an open discrete fracture is assumed 2D here. Therefore, the corresponding governing equation is defined in a local 2D coordinate system. The governing flow equation in fractured media is similar to that presented by several authors [4,73,64,77,30]:  fr 

 fr  o oh0 fr oh0 fr l0 fr ð2bÞ K þ qr gj cos u  S S þ qnjI þ oxi 0 lfr oxj ot  qnjI  ¼ 0;

i; j ¼ 1; 2

ð35Þ

where the last two terms represent normal components of fluid flux across the boundary interfaces (I+ and I) that separate the fracture and the porous matrix. In the conceptual model, fractures are idealized as 2D parallel plates. Therefore, both the total head, hfr0 , and the relative density, qfrr , are uniform across the fracture width. The specific storage in an open fracture, S frS ½L1 , can be derived from (34) by assuming that the fracture is essentially incompressible, such that am = 0, and by setting its porosity to 1: S frS ¼ q0 gafl

ð36Þ

3.2. Reactive solute transport The governing reactive transport equation in porous media has the 3D form [1]   o oC oð/RCÞ ; i; j ¼ 1; 2; 3 /Dij  qi C þ C m ¼ oxi oxj ot ð37Þ where C [M L3] is solute concentration. In this form of the transport equation, the assumptions of fluid incompressibility and constant fluid density are made. The coefficients of the hydrodynamic dispersion tensor, Dij [L2 T1], are given by Bear [1] as /Dij ¼ ðal  at Þ

qi qj þ at jqjdij þ /sDd dij ; jqj

i; j ¼ 1; 2; 3 ð38Þ

750

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

where al [L] and at [L] are the longitudinal and transverse dispersivity, respectively, s [–] is matrix tortuosity, Dd [L2 T1] is the free-solution diffusion coefficient and dij [–] is the Kronecker delta function. The dimensionless retardation factor, R, is given by Freeze and Cherry [23] as q R ¼ 1 þ b Kd ð39Þ / where qb [M L3] is the bulk density of the porous medium and Kd [M1 L3] is the equilibrium distribution coefficient describing a linear Freundlich isotherm. The source/sink term, Cm [M L3 T1], is /kRC for radioactive components with decay constant k [T1]. For chemically reactive species such as silica, the governing transport equation is obtained from (37) by replacing the concentration, C, by the silica molality, m, and by setting the source/sink term, Cm [now MOL M1 T1], equal to the net reaction rate, rnet, given by (2). Yeh and Tripathi [85] argue that precipitation/dissolution reactions and sorption cannot be simulated simultaneously if the aqueous component is the primary dependent species. Thus, the distribution coefficient of silica must be set to zero in Eq. (39). Therrien and Sudicky [76] give the equation that describes 2D solute transport in a discrete fracture as 

 fr fr  o fr fr oC fr fr fr oC ð2bÞ Dij  qi C þ Cm  R þ XnjI þ oxi oxj ot  XnjI  ¼ 0;

i; j ¼ 1; 2

ð40Þ

where Dfrij ½L2 T1  is the hydrodynamic dispersion coefficient of the fracture, calculated as qfri qfrj Dfrij ¼ ðafrl  afrt Þ fr þ afrt jqfr jdij þ Dd dij ; jq j

i; j ¼ 1; 2 ð41Þ

where afrl and afrt [both L] are the longitudinal and transverse fracture dispersivity, respectively. The dimensionless fracture retardation factor, Rfr, is given by [23] Rfr ¼ 1 þ

2K frd ð2bÞ

ð42Þ

where K frd ½L is the fracture–surface distribution coefficient. The source/sink term in Eq. (40) is Cfrm ¼ kRfr C fr for radioactive chemicals and Cfrm ¼ rfrnet for the silica species. The last two terms in (40) represent advective–dispersive–diffusive loss or gain of solute mass across the fracture–matrix interfaces I+ and I [73]. Sorption reactions of silica must be neglected, thus K frd ¼ 0 for silica [85]. Note that Eqs. (39) and (42) simulate a linear sorption isotherm but calculating corrected dissolution rate constants (Eq. (12)) assumes competitive adsorption of cations. Although assumed in Section 2.1, competitive adsorption is not actually simulated but only used in a qualitative sense to quantify the rate-enhancing effect of salt. The reactive source/sink term always consists of a first order reaction term representing precipitation and a con-

stant term of zeroth order describing dissolution. Thus, both solute transport equations in porous and fractured media are linear, allowing a one-step solution. Therefore, neither an iterative operator splitting, two-step scheme nor a computationally demanding fully-coupled, one-step approach are required. Performing a Newton Iteration and formulating Jacobian matrix entries would have highly complicated the model development. 3.3. Heat transfer The convective–dispersive–conductive heat transfer equation in porous media can be written in a form similar to that given by Molson et al. [43] as   oT o  oT ; i; j ¼ 1; 2; 3 k b þ /Dij ql~cl  qi ql~cl T ¼ qb~cb oxi oxj ot ð43Þ where kb [M L T3#1] is the bulk thermal conductivity, q [M L3] is density and ~c ½L2 T2 #1  is specific heat. The absolute temperature, T [#], is the average temperature between the solid and the liquid phase [12]. The subscripts ‘‘l’’ and ‘‘b’’ refer to the liquid and the bulk phase, respectively. It is assumed that the gaseous phase is absent and that external heat sinks and sources due to chemical reactions (dissolution/precipitation) are negligibly small. Bulk properties qb~cb ½M L1 T2 #1  and kb can be quantified considering the volume fractions of the solid and the liquid phases according to Bolton et al. [6] qb~cb ¼ ð1  /Þqs~cs þ /ql~cl k b ¼ ð1  /Þk s þ /k l

ð44Þ ð45Þ

where subscript ‘‘s’’ refers to the solid phase. Heat transport in an open discrete fracture can be described with a 2D equation similar to Eqs. (40) and (43) in the form 

   oT fr o oT fr fr fr fr ð2bÞ k l þ Dij ql~cl  qi ql~cl T  ql~cl oxi oxj ot þ KnjI þ  KnjI  ¼ 0;

i; j ¼ 1; 2

ð46Þ

The last two terms represent convective–dispersive–conductive loss or gain of thermal energy across the fracture–matrix interfaces I+ and I. The temperature is uniform across the fracture width. Furthermore, it is assumed that, along the fracture–matrix interface, the temperature in the fracture and the adjoining matrix are identical. 4. Numerical Formulation 4.1. The FRAC3DVS model FRAC3DVS is a 3D variable-density saturated-unsaturated numerical groundwater flow and multi-component solute transport model. A detailed description of the model can be found elsewhere [76,77,30,31] and is not repeated

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

here. The FRAC3DVS model was modified to simulate the chemical system described in Section 2.1. The extended model couples reactive transport with variable-density, variable-viscosity flow and with changes of solid phase properties. 4.2. Coupling flow and reactive transport The processes of variable-density flow and reactive solute transport are naturally coupled. Density variations cause weak nonlinearities in the flow equation. In the numerical model, they are treated by means of a sequential iterative approach (SIA), also called Picard iteration, which links the two governing equations for flow and transport. This method alternately solves the two governing equations during each time step until convergence is attained. Mineral dissolution/precipitation has a direct impact on a variety of physicochemical and material properties during the simulation. A change of porosity and fracture aperture affects the active surface area, which, in turn, changes the net rate of reaction (1). The change of such parameters is

no

first time step

751

naturally fully coupled with flow, heat transfer and solute transport. However, mineral volume fractions change much more slowly than do the solute concentrations in the fluid [56,65,69,55]. Therefore, in the present model, like in other common geochemical models, fluid properties are updated after each iteration of the Picard loop whereas material properties are updated after each time step rather than after each iteration (Fig. 2). This procedure of recalculating material parameters at the end of each time step is called the quasi-stationary state approximation and has first been introduced by Lichtner [38]. Using the reaction rate at time level L+1 (implicit time weighting scheme) to update all model parameters ensures numerical stability [65]. The decoupled, two step approach to update material properties works well for relatively small time step sizes. However, if nonuniform time step sizes are used to accelerate the simulation, the time increment may become too large. As a consequence, high reaction rates may lead to unrealistically great changes in quartz volume fraction during a single time step. In this case, variable-density, vari-

threshold exceeded ?

update material properties

yes

initial Cσ , T, h 0

radioactive TRANSPORT

next time step

reduce Δt repeat time step

Niterations = 0

reactive TRANSPORT

Picard Iteration Loop no Niterations += 1

convergence AND Niterations > 1

yes

update fluid properties

ρ, μ

FLOW

T

Cσ

HEAT

TRANSPORT

h0

DARCY

q

Fig. 2. Flow chart of the Picard Iteration with chemistry loop to couple variable-density, variable-viscosity flow and solute transport with external chemical reactions and parameter updates.

752

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

able-viscosity flow and reactive transport are not satisfactorily coupled with parameter changes because the time increment is too large. As a solution, the model dynamically reduces the time step size. This adaptive time stepping scheme eliminates unphysically large changes of quartz fraction in order to stabilize the simulation. With adaptive time stepping, the time step size depends on the absolute change of porosity according to ðDtÞ

Lþ1

¼

/ L ðDtÞ max j/Lþ1  /L j

ð47Þ

where /* is the maximum absolute change in porosity allowed during a single time step. Because quartz is considered as the only reactive mineral, the threshold /* applies to absolute changes of both porosity and quartz fraction. The model includes a verification of newly calculated quartz fractions to ensure that they are not negative. If the maximum change in porosity is greater than the allowed threshold, the fraction in Eq. (47) is less than 1 and the updated new time step size is smaller than the previous one. In this case, the old time step is repeated using the new reduced time increment, (Dt)L :¼ (Dt)L+1, without updating the material properties. In fractured systems, the adaptive time stepping can also be based on absolute changes in fracture aperture by using an expression similar to (47). If both time step size controllers (porosity control and aperture control) are used, the new time step size is calculated from the material whose time step multiplier is smaller. Therrien and Sudicky [76] previously used adaptive time stepping in the simulation of variably-saturated flow. However, the time stepping presented above is different from that used by Therrien and Sudicky [76] where time steps are not repeated and where new time step sizes always apply to the following time step. 5. Model verification Simulations are presented here to verify the model formulation for reactive solute transport and heat transfer. The verification problems shown here for a single fracture are identical to the fracture–matrix geometry used by Tang et al. [74]. 5.1. Reactive solute transport This verification problem examines 2D advective–reactive transport in a single fracture, embedded in a porous matrix. It is assumed that solutes migrate due to advection in the fracture and due to molecular diffusion in the porous matrix. Chemical reactions take place in the fracture and in the matrix. Molecular diffusion and mechanical dispersion in the fracture are neglected, allowing an easier formulation of the analytical solution of Tang et al. [74] with no need to numerically integrate. Groundwater in the fracture migrates at a constant velocity. It is assumed that ground-

0 water is free of dissolved electrolytes, thus k corr þ ¼ k þ and cH4 SiO4 ¼ 1. For simplicity, the molal concentration of silica will be written as m. Heat transfer is not considered here and a constant background temperature is imposed. It is further assumed that the material properties (i.e., matrix porosity, hydraulic conductivity, fracture aperture, mineral surface area) are constant in time. Different mineral surface areas in the porous matrix and in the fracture are used, resulting in two different net reaction rates. This assumption does not correspond to Tang et al. [74] where the radioactive decay rates in fracture and matrix are identical. Initially, the entire domain is in thermodynamic equilibrium. Silica-free freshwater enters the fracture at a constant rate during the entire simulation, diluting the silica-saturated fluid in the fracture. All boundaries, except the fracture inlet and outlet, are impermeable for flow and are assigned zero-dispersive transport rates. The resulting drop of silica molality creates a thermodynamic disequilibrium and initiates quartz dissolution. Eventually, the system reaches equilibrium between dilution and dissolution. 0 With m = m 0 + Keq and mfr ¼ mfr þ K eq , the governing equations of this problem using the new variables, m 0 0 and mfr , are given by Steefel and Lichtner [66] in the form: 0

om0 o2 m0 /qz k þ Aqz 0  Dd 2 þ m ¼ 0; ot ox /K eq

b6x61

ð48Þ

and 0 0 fr fr0 fr0 / k A omfr om 0 /D om qz d þ qz þ vfr þ mfr  ot oz K eq b ox ¼ 0;

06z61

x¼b

ð49Þ

for reactive transport in the porous matrix and in the discrete fracture, respectively. The groundwater velocity in the fracture is given by vfr [L T1]. Using the new governing equations (48) and (49), both initial and boundary conditions are identical to those used in Tang et al. [74]. They are formulated mathematically by Steefel and Lichtner [66] who presented the steady state as well as the transient analytical solutions. In the numerical simulation, the finite element domain is similar to the fracture–matrix system used by Tang et al. [74]. It is spatially discretized in the x-direction using a gradually increasing Dx with factor 1.1 from Dx = 0.005 cm near the fracture to Dx = 0.1 cm at the domain boundary. In the flow direction, the Dz increases with factor 1.25 from Dz = 0.001 cm near the source to Dz = 0.1 cm at the domain boundary. All model parameters are summarized in Table 2. Fig. 3 shows the concentration profile versus distance along the fracture for both the analytical and the numerical solution. Because initial thermodynamic equilibrium has been assumed, silica concentrations are initially high and close to the equilibrium constant (Fig. 3a). Inflowing freshwater dilutes the silica-saturated water and decreases silica concentration (Fig. 3b and c). This leads to subsaturated

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771 Table 2 Model parameters used in the verification example for 2D reactive silica transport in fractured porous media

Constant background temperaturea (TC) Matrix porosity (/) Quartz volume fractionb (/qz) Specific surface area in the matrixa (Aqz) Specific surface areain the fracturec Afr qz Free-solution diffusion coefficientd (Dd) Fracture apertured (2b) Groundwater velocity in the fractured (vfr)  Dissolution rate constante k 0þ e Equilibrium constant (Keq) Domain sizea (‘x, ‘z) Location of cross-sections (z1, z2) Output times (t1, t2, t3, t4)

239 C

a b c d e

v

fr

0.35 0.65 54.2 m2 kg1

(b)

6 150

5 1

Si [mmol kg-1]

2

6.15 m kg

1.0 · 1010 m2 s1 200 lm 1.9727 · 105 m s1 1.3298 · 108 mol m2 s1 6.4996 · 103 m kg1 2.0 cm, 3.1 cm 0.1‘z, 0.5‘z 500 s, 1000 s, 2000 s and steady state

fracture

m0 = K

eq

m0fr = K

eq

(b) 4 Analytical

(a)

3

100

Numerical 2

50

1 0

0 0

0.1

0.2

0.3

distance into matrix [cm]

Fig. 4. Concentration profiles of 1D reactive transport of silica in discretely fractured porous media. Shown are the silica molalities in the matrix at steady state at the distances (a) 0.31 cm and (b) 1.55 cm from the solute source.

sion. However, as the simulation proceeds in time, this inconsistency diminishes and eventually vanishes after an infinitely long period of time. Perfect match between the analytical solution and the results from this model are obtained with the molal concentrations in the matrix as shown in Fig. 4.

7

5.2. Heat transfer

6

(a)

(b)

150

(c)

5 4

(d)

100

3 2

50

Analytical

1

Si [PPM]

Si [mmol kg-1]

m0fr = K eq

fracture

location of cross-sections

porous matrix m =0

v

7

Johnson et al. [34]. 1  /. From Eq. (27) with x = 1.0. Steefel and Lichtner [66]. Computed by this model for deionized water at TC = 239 C.

fr 1

m =0

fr

Si [PPM]

Value

m0 = K eq

porous matrix fr 1

(a)

Parameter

753

Numerical

0

0 0

1

2

3

distance along fracture [cm] Fig. 3. Concentration profiles of 1D reactive transport of silica in discretely fractured porous media. Shown are the silica molalities in the fracture at (a) 500, (b) 1000 and (c) 2000 s and at (d) steady state.

This test case verifies 2D heat transfer in a single fracture embedded in a porous matrix. This example is based on analytical results presented by Meyer [42] for advective heat transfer in a fracture coupled with heat conduction in a surrounding porous matrix. Mechanical heat dispersion as well as conduction within the fracture are not considered, making numerical integration unnecessary. The groundwater flow velocity in the fracture is constant. Under these assumptions, the governing equations for heat transport in the matrix and in the discrete fracture simplify from (43) and (46) to qb~cb

oT o2 T  k b 2 ¼ 0; ot ox

b6x61

ð50Þ

and conditions promoting dissolution of quartz. Because the ambient temperature is high (239 C), silica concentration changes within minutes [34]. Eventually, the system reaches a steady-state condition where dilution and dissolution are in equilibrium (Fig. 3d). The discrepancy between the analytical and numerical solution at early times was previously described by Steefel and Lichtner [66], who interpreted this as numerical disper-

fr oT fr k b oT fr fr oT ql~cl þ ql~cl v  ¼ 0; ot oz b ox x¼b

06z61

ð51Þ

The last term in (51) expresses conductive loss of heat from the fracture into the matrix on the fracture–matrix interface. Initially, the entire system has a uniform temperature equal to T0. The fluid entering the fracture has the constant temperature equal to T1. All boundaries, except the fracture

754

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

inlet and outlet, are impermeable for groundwater flow and for heat exchange. According to Meyer [42], the transient solution along the fracture is ! pffiffiffiffiffiffiffiffiffiffiffiffiffiffi z k b qb~cb T fr  T 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ erfc ð52Þ T1  T0 2vfr ql~cl b ðt  z=vfr Þ

ð53Þ The fracture–matrix geometry is identical to that used by Tang et al. [74]. The finite element domain was spatially

o

T1 = 15 C

v (a)

fr

T0 = 10 C

fracture

fr 0

o

T = 10 C

(b)

location of cross-sections 15

14

temperature [°C]

Using the analytical results presented by Tang et al. [74], it can be shown that the transient solution along a crosssection from the fracture into the porous matrix is given by ! pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi qb~cb ðx  bÞ z k b qb~cb T fr  T 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ pffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ erfc T1  T0 2vfr ql~cl b ðt  z=vfr Þ 2 k b ðt  z=vfr Þ

o

porous matrix fr

Analytical

13

Numerical

(a) 12

(b)

11

10

Table 3 Model parameters used in the verification example for 2D heat transfer in a single fracture embedded in a porous matrix Parameter

Value 3

Bulk thermal conductivity (kb) Heat capacity of solid ð~cs Þ Solid density (qs) Heat capacity of water ð~cl Þ Fluid density (ql) Matrix porosity (/) Groundwater flow velocity in the fracture (vfr) Initial temperature (T0) Boundary temperature (T1) Domain size (‘x, ‘z) Location of cross-sections (z1, z2) Output times (t1, t2)

1

3.4 kg m s K 908 m2 s2 K1 2550 kg m3 4192 m2 s2 K1 997 kg m3 0.2 0.05 m s1 10 C 15 C 2 m, 10 m 0.1 m, 0.61 m 5000 s and 10,000 s

All parameters are identical to those used by Meyer [42].

0

0.1

0.2

0.3

0.4

0.5

distance into matrix [m]

Fig. 6. Temperature profiles of 1D heat transfer in discretely fractured porous media. Shown are the temperatures in the matrix at 10,000 s simulation time at the distances (a) 0.1 m and (b) 0.61 m from the heat source.

discretized in the x-direction by gradually increasing Dx with constant factor 1.1 from Dx = 0.01 m near the fracture to D x = 0.1 m at the domain boundary. In the flow direction, D z also increases gradually from Dz = 0.1 m near the elevated temperature to Dz = 0.5 m at the domain boundary. All other parameters are presented in Table 3 and the simulation results are exhibited in Figs. 5 and 6. 6. Illustrative examples

o

porous matrix fr

T1 = 15oC

vfr

T0 = 10 C

To simulate heat transfer and chemical reactions in a fracture network, existing studies of reactive transport in porous media must first be expanded to include heat transfer. Section 6.2 presents simulations for an unfractured porous medium, pm, that focus on (i) reactive transport, (ii) thermohaline variable-density, variable-viscosity flow, and (iii) coupled thermohaline reactive transport. These processes will in turn be considered for a fractured medium, fm, in Section 6.3. In Section 6.4, the long-term effect of chemical reactions on density-driven mass flux through fractured porous media, fm_long, will be demonstrated. The description of problems pm and fm is given in Section 6.1 while problem fm_long will be described in Section 6.4.

T0fr = 10oC

fracture

15

temperature [°C]

14

13

Analytical Numerical

12

(a) (b)

11

6.1. Problem description 10 0

1

2

3

4

5

6

7

8

9

10

distance along fracture [m]

Fig. 5. Temperature profiles of 1D heat transfer in discretely fractured porous media. Shown are the temperatures in the fracture at (a) 5000 and (b) 10,000 s.

The domain for simulations pm and fm is the vertical cross-section shown in Fig. 7. A similar domain has been used for simulations presented by Schincariol et al. [61], Ibaraki [33] and Freedman and Ibaraki [21] and it corresponds to the laboratory tank used in the experimental

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

755

source of constant concentration and temperature

1.0 m 0.085 m 0.03 m 0.25 m 0.135 m

flow

Fig. 7. Model domain and location of the solute source for numerical simulations of reactive silica transport and variable-density thermohaline flow. The domain is used for simulations in unfractured and discretely fractured porous media.

work of Schincariol and Schwartz [60]. This domain was chosen to highlight the coupling between variable-density flow and reactive solute transport as done previously by Freedman and Ibaraki [21]. The 2D simulation domain has dimensions of 1.0 m · 0.25 m with a unit thickness. The domain has been spatially discretized using 63,000 rectangular finite elements, which are smaller at the left boundary (Dx = 0.5 mm, Dz = 2.0 mm), and whose size increase towards the right (Dx = 2.0 mm, Dz = 2.0 mm). Fractures are shown in Fig. 7 but the first simulations presented here are for an unfractured porous medium. Simulations that incorporate fractures are presented later. We consider that the domain has a uniform initial temperature of 239 C [34] and that the only dissolved species initially in the fluid is aqueous silica (H4SiO4). The fluid is initially in thermodynamic equilibrium (rnet = 0), where mH4 SiO4 ¼ K eq =cH4 SiO4 ¼ 6:4996 mmol kg1 is the initial silica molality at Cr = 0.0 mg l1 and TC = 239 C. A horizontal flow field is imposed by assigning constant fluid fluxes (q = 1.045 · 106 m s1) along the left and right boundaries, with top and bottom boundaries being impermeable to fluid flow. Also shown in Fig. 7 is the location of a source of fluid entering the domain at a constant concentration and temperature. It is assumed that this fluid contains the four common ions Na+, Ca2+, Mg2+ and Cl. This choice of fluid composition is based on average concentrations of major ions in groundwater of the Canadian Br

SO4

HCO3

Cl

TDS

Shield at different depths (Fig. 8). In that environment, water below 1000 m is a Ca–Na–Cl brine with dissolved solids exceeding 100,000 mg l1. Stober and Bucher [71] have stated that all deep waters below 500 m in the continental crystalline crust are brines of the same chemical composition. Although the deep water is depleted in Mg2+, this cation was used as the fourth mobile species because Dove [13] has shown that even at low concentrations, Mg2+ can significantly enhance silica dissolution rates. Fig. 8 also illustrates that water at depth 500 m is rich 2 in SO2 4 . According to Marshall and Chen [40], SO4 forms a sulphate–silicic-acid complex at temperatures above 150 C, which increases silica solubility. The significance of sulphate was verified by calculating the silica activity coefficient with typical Na+, Ca2+, Mg2+, Cl and SO2 4 concentrations at 500 m depth and at the temperature maximum (300 C), where the solubility increasing effect of sulphate is highest. This verification indicated that the activity coefficient is approximately equal to 1 and the presence of SO2 4 is therefore assumed not to influence silica dissolution and it is neglected here. In the numerical simulations, first-type boundary conditions are imposed at the source for solute and heat transport, with a constant concentration equal to 1000 mg l1 assigned to the four species Na+, Ca2+, Mg2+ and Cland a constant temperature equal to TC = 247 C, higher than the initial temperature. First-type boundary conditions are assigned to the remainder of the left boundary of the Sr

K

1E+00

1E+01

Mg

Na

Ca

TDS

depth [m]

0

500

1000

1500 1E+00

1E+01

1E+02

1E+03

1E+04

concentration [mg/l]

1E+05

1E+06

1E+02

1E+03

1E+04

concentration [mg/l]

Fig. 8. Average ion concentrations in groundwater of the Canadian Shield [17].

1E+05

1E+06

756

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

domain, outside the source, with zero concentration and a constant temperature equal to 239 C, to ensure that the only solute input in the domain is through the source. The top, bottom, and right boundaries are zero-conductive heat transfer and zero-dispersive solute flux boundaries. Thus, heat and solutes cannot cross the top and bottom boundaries but they are able to cross the right boundary by convection and advection, respectively. As opposed to the four source species, the aqueous silica molality at the source is not constant but recalculated at each time step. All boundaries are zero-dispersive flux boundaries for aqueous silica. Initially, the equilibrium silica molality at the source is mH4 SiO4 ¼ K eq =cH4 SiO4 ¼ 6:9361 mmol kg1 at Cr = 1000.0 mg l1 and TC = 247 C. The choice of initial thermodynamic equilibrium makes it easy to identify every deviation from the silica equilibrium as the result of a chemical reaction (dissolution or precipitation). Because silica is the only species whose concentration will not be influenced by water density, we will give silica concentrations in the density-independent unit molality (mol kg1). The simulations cover a time of 3 days with increasing time step sizes. Thermal deformations of the rock are not considered. The spatial and temporal discretization as well as all simulation parameters are summarized in Table 4. 6.2. Coupled thermohaline and reactive transport in porous media

Table 4 Model parameters used in reactive transport studies Parameter Domain size (‘x, ‘z) Spatial discretizationd (Dx, Dz) Temporal discretizatione (Dt) Longitudinal dispersivitya,b,c (al) Transverse dispersivitya,b,c (at) Tortuosityb,c (s) Average Darcy fluxa,b,c (q) Free-solution diffusion coefficientb,c (Dd) Distribution coefficient (Kd)

Reference fluid densityg (q0) Reference fluid dynamic viscosityg (l0) Fluid compressibilityh (afl) Matrix compressibilityh (am) Initial porosityb,c (/init) Initial hydraulic freshwater  conductivityb K 0;init ij Initial specificsurface  area in the matrixi Ainit qz Solid phase densityj (qs) Specific heat of solidj ð~cs Þ Specific heat of liquidj ð~cl Þ Thermal conductivity of solidj (ks) Thermal conductivity of liquidj (kl) a b c d

The first simulation, entitled pm_reac, is for reactive transport but assumes that the fluid has constant density and viscosity, equal to those of the ambient groundwater. In pm_reac, the time step size changes dynamically, based on porosity changes. The maximum permitted change in porosity per time step, /*, was set to 103 (0.1%). The initial and maximum time step sizes chosen were 1 min and 2 h, respectively. Fig. 9 shows simulation results for pm_reac after 3 days. Thermal energy is predominantly transferred by conduction in both the longitudinal and the transverse direction (Fig. 9a). Because buoyancy forces are not considered in pm_reac, the plume is mainly transported by advection and migrates laterally across the domain. Molecular diffusion and transverse dispersion slightly increase the plume extension in the vertical direction. The chloride concentration (Fig. 9b; no retardation) and the magnesium concentration (Fig. 9c; highest retardation) illustrate this transport behavior. Concentration contours for the nonreactive and nonsorptive chloride indicate the position of the advective front. Fig. 9d shows the molal concentration of silica and reveals an interesting simulation result for pm_reac. Near the source, temperatures are relatively high such that quartz dissolves. However, further away from the source, the temperatures are close to the background temperature, 239 C and solute concentrations are high, which decreases the solubility of silica. Therefore, salinity controls the silica concentration further away from the source. Conversely,

Value a

e f g h i j

1.0 m, 0.25 m 0.5 mm, . . . , 2.0 mm, 2.0 mm 1 min, . . . , 2 h 3.0 · 104 m 0.0 m 0.35 1.045 · 106 m s1 1.6 · 109 m2 s1 [Na+] 3.0 · 106 kg1 m3 [Ca2+] 5.0 · 105 kg1 m3 [Mg2+] 1.0 · 104 kg1 m3 [Cl] 0.0 kg1 m3 [H4SiO4]f 0.0 kg1 m3 815.969 kg m3 1.1184 · 104 kg m1 s1 4.4 · 1010 kg1 m s2 1.0 · 108 kg1 m s2 0.38 5.6 · 104 m s1 54.2 m2 kg1 2650 kg m3 738 J kg1 K1 4186 J kg1 K1 5.0 W m1 K1 0.6 W m1 K1

Freedman and Ibaraki [21]. Schincariol et al. [61]. Ibaraki [33]. To fulfill the Pe´clet criterion, Pe < 2.3, used by b and c. To fulfill the Courant criterion, Cr 6 1.0, used by b. Yeh and Tripathi [85]. Computed by this model for deionized water at TC = 239 C. Shikaze et al. [64]. Johnson et al. [34]. Bolton et al. [6].

the silica concentration follows the isotherms near the source as well as in regions of low salinity above and below the plume. Fig. 10 is a vertical cross-section located at a distance of x = 0.12 m from the source. The figure shows that the silica concentration is proportional to temperature in low-salinity zones and inversely proportional to salinity in high-salinity zones. Clearly, these observations demonstrate the solubility-lowering effect of salt and the solubility-increasing effect of temperature as discussed by Dove [13]. Fig. 9e1 finally shows the distribution of the hydraulic freshwater conductivity. As expected, the area around the source became more conductive because of quartz dissolution. However, the elongated bluish fields located to the right of the source indicate a conductivity value smaller than the initial one. Apparently, dissolved silica is transported by advection to the right. Silica is assumed to be nonsorptive and its transport rate is, therefore, comparable to the chloride transport rate. If silica molecules are 1 For interpretation of colour in this figure, the reader is referred to the web version of this article.

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

757

Fig. 9. Results of reactive transport simulations in an unfractured porous medium with constant fluid density and viscosity (pm_reac). Shown are (a) temperature, (b) chloride and (c) magnesium ion concentration, (d) molal concentration of aqueous silica and (e) freshwater hydraulic conductivity at 3 days.

transported laterally to regions of lower temperature and high salinity, the system becomes locally supersaturated and some of the previously dissolved silica precipitates, resulting in lower hydraulic conductivity. The second simulation, called pm_dens, is for variabledensity flow in an unfractured porous medium but ignores chemical reactions. Water density and viscosity are calculated from temperature and salinity. Unlike the previous simulation, time step sizes are prescribed and gradually increase from 1 min to 2 h. The results show that the magnitude of buoyancy is controlled by water density, which is a function of both temperature (Fig. 11a) and salinity (Fig. 11b and c). Fig. 11b demonstrates that density effects cause vertical fluid movement. The figure shows concentration profiles of Cl at 3

days, highlighting the mixed convective flow character. Forced convection (advection) remains the main lateral transport mechanism whereas buoyancy-induced free convection controls the shape of the plume in the vertical direction. Different diffusivities explain the completely different transport behavior of thermal energy and solutes, which is known as double diffusive transport [48]. In the pm_dens simulation, heat transfer is practically independent of groundwater flow, while water flow dominates solute transport. This difference results in an interesting density distribution (Fig. 11d). In the near field of the source, temperature appears to control water density, while in the far field, the salt concentration has the greatest influence on water density. The inflowing hot saline water has a density of 808.443 kg m3, lower than the reference density

758

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771 0.25

silica chloride 0.2

zelevation [m]

temperature 0.15

0.1

0.05

0 6.48

6.5

6.52

6.54

6.56

-1

silica [mmol kg ]

0

200

400

600

800

1000

-1

chloride [mg l ]

239.7

239.8

239.9

240.0

temperature [°C] Fig. 10. Vertical cross-section at x = 0.12 m from the source for the simulation pm_reac at 3 days.

(815.969 kg m3). In this case, the inflowing water is less dense than the ambient water. As a consequence, the relative density, qr, is negative (9.223 · 103), resulting in a positive buoyancy effect near the source. However, further away from the source, the influence of the solutes on density dominates because advective solute transport is more efficient than conductive heat transfer. Therefore, the water density exceeds its reference value and the density contrast is positive (in the range of qr  103), which results in a noticeable sinking of the plume. The last simulation for an unfractured porous medium, pm_reac_dens, couples the effect of density with chemical reactions. Similar to the pm_reac simulation, hydraulic conductivity and matrix porosity change with time as a result of reactions. Time step sizes adapt to porosity changes with the maximum permitted porosity change, /*, set to 103 (0.1%). This contrasts with pm_dens, where conductivity and porosity remain constant over time. Fig. 12 shows that, as before, the temperature (Fig. 12a) is the important factor in the near field of the source where it controls quartz solubility (Fig. 12d) and water density (Fig. 12e). In the far field, however, the salt content (Fig. 12b and c) dictates variable-density flow and chemical reactions. Although conductivity and permeability vary with time in pm_reac_dens and are constant in pm_reac, the results of the two simulations are not significantly different. This observation is in agreement with findings by Freedman

and Ibaraki [21], who simulated the chemistry of calcite, coupled with density-driven flow. They concluded that time scales of a few days are too short to perceive a major impact of the reactions. In addition, the solid phase of both the calcite system studied by Freedman and Ibaraki [21] and the quartz-water system studied here have a fairly low solubility. Therefore, coupling between flow and reactive transport is weak. Nevertheless, the three simulations presented here illustrate the coupling between variable-density flow, heat transfer and reactive transport in porous media. They also show that adaptive time stepping is a useful tool and certainly competitive compared with the conventional use of predefined time step sizes (Table 5). 6.3. Coupled thermohaline and reactive transport in fractured media A second series of three simulations assumes the presence of fractures (fm) oriented transversely to the ambient flow direction as shown in Fig. 7. Using the random fracture generator developed by Graf and Therrien [31], a total of 60 random fractures are generated following the two main orientations 60 and 120 with the standard deviation of the Gaussian distribution, sigma = 1. All fractures are 0.1 m in length and have a uniform aperture equal to 100 lm. Initial and boundary conditions of the fm simulations are identical to those used in the previous example in

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

759

Fig. 11. Results of density dependent nonreactive transport simulations in porous media (pm_dens) at 3 days. Shown are (a) temperature, (b) chloride and (c) magnesium ion concentration and (d) water density.

porous media. Table 4 presents the simulation parameters while Table 6 shows additional parameters for the simulations including fractures. The first simulation, fm_reac, ignores density effects but simulates reactive transport. The time step sizes adapt to changes in matrix porosity and/or fracture aperture. Maximum permitted changes of porosity and aperture are chosen as /* = 103 (=0.1%) and (2b)* = 0.1 lm, respectively. Fig. 13 shows the results after 3 days. The fractures have a substantial impact on fluid migration because their hydraulic conductivity is more than 100 times greater than that of the porous matrix between the fractures. Because the value of matrix permeability has originally been used for simulations without fractures, the permeability ratio of 100 is not typical for fractured rocks. Other flow studies in fractured rocks [64,30,31] assume a ratio that is more than three orders of magnitude greater than that observed here. In these cases, fractures have a most significant impact on flow. Nevertheless, the fractures do increase the transverse dispersion of the plume. This results in a larger vertical extension of the plume and reduced lateral migration (Fig. 13b and c), compared with results in porous media (Fig. 9b and c). Fig. 13d exemplifies the thermohaline influ-

ence on silica solubility. The plume is now more dispersed and solute concentrations in the far field are low, causing less silica precipitation. Therefore, silica precipitation does not significantly lower hydraulic conductivity (Fig. 13e). The second simulation, fm_dens, is for variable-density flow but ignores chemical reactions. Time step sizes are prescribed and gradually increase from 1 min to 2 h. Results for this simulation are shown in Fig. 14. Because heat conduction is generally the principal heat transfer mechanism, fractures do not affect the temperature distribution (Fig. 14a). The simulated concentration profile of the nonreactive, nonsorptive chloride ion at 3 days (Fig. 14b) indicate that, in the present case, fractures do not act like preferential pathways as in simulations by Shikaze et al. [64] and Graf and Therrien [30,31]. However, because the plume is more dispersed in the fm_dens simulation (Fig. 14b) than in the pm_dens simulation (Fig. 11b), water density in the far field is smaller in fractured media (816 kg m3) than in porous media (817 kg m3). As a consequence, density contrasts in the far field of the fractured medium are generally small and in the range of qr  104. Therefore, the buoyancy effect in the fractures is minor. This contrasts with simulations carried out by

760

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

Fig. 12. Results of density dependent reactive transport simulations in porous media (pm_reac_dens) at 3 days. Shown are (a) temperature, (b) chloride and (c) magnesium ion concentration, (d) molal concentration of aqueous silica and (e) freshwater hydraulic conductivity.

Table 5 Simulations and CPU times in porous media (pm) Simulation pm_reac pm_dens pm_reac_dens a

Chemical reactions p – p

Density variations – p p

Time stepping Adaptive Prescribed Adaptive

CPU timea 32 min 1 h 24 min 1 h 32 min

Computed on a Pentium 4, 2.6 GHz, 500 MB RAM.

Shikaze et al. [63] and Graf and Therrien [30,31] where relative density has values up to 0.2, more than three orders of magnitude greater than those encountered here. The last fm_reac_dens simulation is for variable-density flow and reactive transport using adaptive time stepping as in fm_reac. The dissolved ions are now distributed over a larger cross-sectional area with smaller salt concentrations

Table 6 Additional model parameters used in reactive transport studies in fractured media Parameter

Value

Fracture dispersivitya,b (afr) Initial fracture apertureb,c (2b)init   Initial specific surface area in the fractured Afr;init qz Fracture roughness coefficiente (x)

0.1 m 100 lm 49.021 m2 kg1 1.0

a b c d e

Therrien and Sudicky [76]. Tang et al. [74]. Sudicky and Frind [24]. Eq. (27). Steefel and Lasaga [65].

(Fig. 15b and c) than before (Fig. 12b and c). The silica concentrations in the far field are smaller than in the pm_reac_dens example in unfractured porous media

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

761

Fig. 13. Results of nondensity dependent reactive transport simulations in fractured media (fm_reac) at 3 days. Shown are (a) temperature, (b) chloride and (c) magnesium ion concentration, (d) molal concentration of aqueous silica and (e) freshwater hydraulic conductivity.

(Fig. 15d), which is a result of the silica solubility-lowering effect of dissolved salt. Because silica precipitation is reduced in fractured media, the hydraulic conductivity of the porous matrix does not change significantly (Fig. 15e). Aperture changes are insignificant for the time scale of this simulation and are not shown. The CPU times are typically greater than in the pm simulations because fractures are present (Table 7). It is again shown that adaptive time stepping is a helpful means to accelerate and control the simulation process. 6.4. Long-term simulation of coupled thermohaline flow and reactive transport in fractured porous media Long-term simulations in a fractured medium (fm_long) have been carried out to demonstrate the significant effect

of chemical reactions on thermohaline flow. The simulation domain is similar to that used by Shikaze et al. [64] and shown in Fig. 16. The domain of dimension 10 m · 10 m is spatially discretized into rectangular elements of size Dx = Dz = 0.1 m. Mobile species are Na+, Cl, H4SiO4 and temperature. Initially, the system is free of salt, silica concentration corresponds to thermodynamic equilibrium, ensuring reaction rates equal to zero in the entire domain and temperature is uniform equal to 10 C. Boundary conditions are also shown in Fig. 16. Lateral boundaries are impermeable for flow, while both top and bottom boundaries are assigned constant hydraulic heads equal to zero. The domain top is assigned a constant temperature equal to 10 C and a constant concentration of 0.0 mg l1 for both Na+ and Cl, corresponding to freshwater at air temperature. The bottom of the domain is

762

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

Fig. 14. Results of density dependent nonreactive transport simulations in fractured media (fm_dens) at 3 days. Shown are (a) temperature, (b) chloride and (c) magnesium ion concentration and (d) water density.

assigned a constant temperature of 239 C and a constant concentration of 1000 mg l1 for both Na+ and Cl, corresponding to salty hot groundwater. Unlike salt concentrations and temperature, the silica molality at bottom and top is not imposed but recalculated at each time step. A temperature of 239 C has been chosen because (i) it has been used by Johnson et al. [34] to experimentally study quartz chemistry, and (ii) temperature around a radioactive waste facility is in the range of 239 C [10,84]. With concentrations of 1000 mg l1, water density at the bottom of the domain is 818.137 kg m3, giving a relative density value of 0.19, which is similar to the value of 0.2 used by Shikaze et al. [64]. Flow and transport parameters are identical to those used by Shikaze et al. [64]; Heat transfer and reaction parameters are identical to those used in the previous two sections and given in Table 4. The maximum permitted changes (per time step) of porosity, /*, and fracture aperture, (2b)*, are 0.005 (0.5%) and 1.0 lm, respectively. Simulation time is 100 years. The 2D vertical slice shown in Fig. 16 represents the typical situation of cold freshwater above hot saltwater found in many aquifers. The block is assumed to be heated from below due to nuclear fuel waste heat generation. The spatial dimension of 10 m · 10 m is atypical for a geological

formation above a deep radioactive waste facility. However, this dimension has been chosen in order to create an illustrative scenario of unstable density-driven flow that is similar to that studied by Shikaze et al. [64]. Fig. 17 displays simulated chloride concentrations at different times. Flow in the entire domain is upwards because denser fluid (cold freshwater) overlies less dense fluid (hot saltwater). Thus, chloride is transported by upwards density-driven advection and a steady-state is reached at about 30 years where the chloride concentration in nearly the entire domain is close to 1000 mg l1. Inspection of the velocity field indicated that convection does not occur because isotherms are horizontal, leading to negative relative densities throughout the domain. Fig. 18 shows results of the fm_long simulation after 100 years. The Fig. illustrates thermal and chemical profiles along the fracture at x = 3.5 m. Because heat conduction is the most important heat transfer mechanism, the steady-state temperature distribution is linear between 239 C (bottom) and 10 C (top) and reached after about 2 years. The hot zone near the bottom corresponds to quartz dissolution. Therefore, silica molality in the dissolution zone is high (about 5 mmol kg1), reaction rates are positive ( ) and the fracture widens up to 247 lm. Dis-

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

763

Fig. 15. Results of density dependent reactive transport simulations in fractured media (fm_reac_dens) at 3 days. Shown are (a) temperature, (b) chloride and (c) magnesium ion concentration, (d) molal concentration of aqueous silica and (e) freshwater hydraulic conductivity.

Table 7 Overview of the simulation trials in fractured porous media (fm) Simulation fm_reac fm_dens fm_reac_dens a

Chemical reactions p – p

Density variations

Time stepping

CPU timea

– p p

Adaptive Prescribed Adaptive

36 min 1 h 32 min 1 h 48 min

Computed on a Pentium 4, 2.6 GHz, 500 MB RAM.

solved silica is transported by upwards density-driven advection through the porous matrix and the fractures. At z P 0.5 m, the combination of low temperature and upward migrating silica creates local supersaturation, causing silica precipitation with negative reaction rates ( ) and reduced fracture apertures down to 39.7 lm. The major

aperture reduction is the result of a self-sealing process of fractures, which is due to the coexistence of a dissolution zone (where silica is produced) and a precipitation zone (where imported silica precipitates). As shown for the simulations presented in Sections 6.2 and 6.3, the model domain is divided into a near field of quartz dissolution and a far field of silica precipitation. A second simulation was run where chemical reactions were neglected. The chloride mass flux through a cross-section at z = 5 m (shown in Fig. 16) was compared for the simulations without and with reactions and the comparison is shown in Fig. 19. For the case without reactions, mass fluxes through fractures, porous matrix and total mass flux increase at early times and attain a plateau after 30 years. This corresponds to the time when both temperature and

764

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771 h0 = 0

discrete fractures

Cσ = 0 o

T = 10 C 10 impermeable lateral boundaries

z - distance [m]

8

6 horizontal cross-section

4

2

h0 = 0

0

-1

Cσ = 1000 mg l 0

2

4

6

8

10

o

T = 239 C

x -distance [m]

Fig. 16. Simulation geometry and boundary conditions for the long-term simulation of coupled thermohaline flow and reactive transport (fm_long).

salinity have reached steady-state and when density-driven upwelling is therefore uniform in the entire domain. However, in the case with reactions, mass fluxes are significantly lower after 10 years. Ultimately, at 100 years, mass fluxes through the fractures of the reactive domain are about 25% lower compared to the nonreactive case. The reason is that the minimum fracture aperture of 39.7 lm is the limiting factor that controls mass flux through the fracture. With time, the ongoing dissolution and precipitation will further enhance fracture self-sealing such that fracture mass fluxes decrease even more. Fig. 19 demonstrates that the impact of reactions and flow parameter changes on density-driven flow is a long-term effect that can only be observed after decades. 7. Sensitivity analysis Additional simulations were performed to assess the impact of parameter uncertainties on reactive solute transport. The verification example from Section 5.1 was used as a base case with original input parameter values. Starting with the base case, two simulations were run for every parameter tested, using a lower and a higher parameter value compared to the base case. We conducted two sensitivity analyses: (i) a mathematical analysis and (ii) an analysis for visualization purpose. All simulations of the mathematical sensitivity test were characterized with a dependent variable that appropriately describes the simulation. This quantity is named n and the terms nlow, norg and nhigh denote the results when simulating with a low, unmodified and high value of parameter p: plow, porg and phigh, respectively. A dimensionless mathematical sensitivity of the model parameter p, Xp, is evaluated using an equation presented by Zheng and Bennett [87]:

Xp ¼

on=norg op=porg

ð54Þ

According to Zheng and Bennett [87], the partial derivative of the dependent variable, n, with respect to the input parameter, p, is normalized by the original value of the variable, norg, and the parameter, porg. The steady state silica concentration at the fracture outlet (z = 3.1 cm) was chosen as the dependent variable, n. The simulation of the base case yields the characteristic number norg = 3.5887 · 103 mol kg1. The reactive transport simulation was run with modified values of four parameters: the specific quartz surface area in the fracture and matrix, quartz volume fraction and temperature. The choice of the range over which the input parameter is varied, is subjective. However, if parameter changes (i.e., perturbations) are too small, computer round-off errors may conceal differences of the dependent variable. On the other hand, perturbations which are too large may yield inaccurate sensitivities, especially if the relation between dependent variable and parameter is nonlinear. In the present mathematical sensitivity analysis, a uniform perturbation size of 5% is applied as suggested by Zheng and Bennett [87]. The results of the mathematical sensitivity test are shown first. Table 8 presents the original value of each input parameter, which was lowered and increased uniformly by 5%. Fig. 20 shows the calculated dimensionless sensitivity for each input parameter. In addition to the mathematical sensitivity analysis, further simulations were carried out with much wider ranges of the input parameters (Table 8), in order to visualize parameter sensitivity. The perturbations are not identical for all parameters and are not used in a mathematical

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

765

Both the surface area in the matrix and the quartz volume fraction are directly proportional to the reaction rate (Eq. (2)). Thus, the sensitivity of the two parameters is similar (Fig. 21b and c). However, uncertainties of the quartz volume fraction also impact the fracture reaction rate, which may cause the slightly higher sensitivity (0.34) than that of the matrix surface area (0.32). Temperature variations have the most significant influence on results (Fig. 21d). If, for example, the ambient temperature increases, the dissolution reactions proceed faster (Arrhenius equation) and, in addition, more quartz dissolves because quartz solubility increases. Both geochemical processes considerably increase the net reaction rate, hence the high sensitivity of temperature with a value of 4.28. 8. Summary and conclusion

Fig. 17. Chloride concentration at (a) 0.5 years, (b) 5 years and (c) 30 years. Concentration contours are 200 mg l1 (thin lines) and 600 mg l1 (bold lines).

sense. The result of this visual sensitivity analysis is shown in Fig. 21. Both sets of simulations show that the uncertainties of the fracture surface area have a negligible impact on the results (Fig. 21a), expressed by the low sensitivity of 0.01. However, the fracture surface area is about one order of magnitude smaller than in the matrix. Consequently, the fracture reaction rate is also one order of magnitude slower. Therefore, the fast reaction in the matrix is dominant and suppresses unprecise fracture reaction rates. However, the sensitivity might be greater than 0.01 if more than one fracture were present or if fracture surface area was in the same order of magnitude than matrix surface area.

Underground disposal of nuclear waste has been determined by the international scientific community as the best option for permanently isolating high-level radioactive waste from the biosphere. Countries like Finland, Sweden and Japan regard low-permeability fractured crystalline rock as a potential host rock for spent fuel. In the present study, a new numerical model has been presented that can be used to assess potential sites for underground disposal of nuclear waste in these countries. The model is free software when used for research purpose. Interested readers are welcome to contact one of the authors to get a DOS or LINUX version of the model along with input files of verification problems and illustrative examples. The new model simulates coupled variable-density, variable-viscosity flow and kinetically controlled reactive transport in nonisothermal fractured porous media. We focused on the chemistry of the common quartz-water system with aqueous silica as the only mobile reactive species. The flow equation is linked with the heat transfer and solute transport equations through an iterative Picard approach. After each iteration, the fluid properties density and viscosity are updated from the individual ion concentrations and temperature as primary variables. Chemical reactions are calculated outside the Picard Iteration because the reactive species silica does not significantly impact the fluid properties. According to the quasi-stationary state approximation [38], flow and reactive transport parameters are also updated at the end of a time step. The model calculates water density from temperature using different equations for different temperature ranges. The dissolved solid effect on density is accounted for with an empirical relationship calibrated by means of density calculations with the VOPO model [44,45]. Fluid viscosity is also evaluated differently for various temperature ranges. The developed model applies the Jones–Dole equation to represent the salt impact on viscosity. The new model also accounts for the reaction rate enhancing effect of dissolved solutes and elevated temperature using a modified form of the Arrhenius equation. The equilibrium constant of the

766

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

5

T H4SiO4

flow

z elevation [m]

4 3

precipitation zone

rnetfr (2b)

2 1

dissolution zone

0 temperature [oC] 150

100

200

net rate [10-8 mol kg-1 sec-1]

250

-1.75

-1 -0.5

0

0.5 1

1.5

H4SiO4 [mmol kg-1] 1

2

3

4

aperture [μm]

5

40

100

150

200

250

Cl-1 mass flux [g sec-1]

Fig. 18. Profiles of temperature, silica molality, fracture reaction rate and fracture aperture up to elevation z = 5 m in the fracture located at x = 3.5 m and for t = 100 years. 0.2

total 0.15

without reactions

fractures

0.1

with reactions

0.05

porous matrix 0 0

25

50

75

100

time [yr] Fig. 19. Chloride mass flux at z = 5 m through the porous matrix, fractures and total mass flux versus time for simulations without (bold lines) and with (thin lines) chemical reactions.

Table 8 Model parameter modifications used for visualization only in the sensitivity analysis of reactive solute transport Parameter Specific quartz surface area in the fracture ðAfr qz Þ Specific quartz surface area in the matrix (Aqz) Quartz volume fraction (/qz) Temperature (TC)

Low value 2

1

Original value 2

1

High value 2

1

1.15 m kg

6.15 m kg

11.15 m kg

34.2 m2 kg1

54.2 m2 kg1

74.2 m2 kg1

0.6

1.0

1.0

209 C

239 C

269 C

Specific quartz surface area in 0.01 the fracture Specific quartz surface area in the matrix

0.32

Quartz volume fraction

0.34

4.28

Temperature 0

chemical system is calculated from the local temperature with a common equation presented by Rimstidt [52]. We adapted the method presented by Marshall and Chen [40] to calculate activity coefficients as a function of salinity, accounting for the haline influence on quartz solubility. Adaptive time stepping is used to calculate the time step sizes in the numerical simulations. New time increments depend on maximum changes in matrix porosity and/or fracture aperture. Adaptive time stepping eliminates

1

2

3

4

5

Dimensionless sensitivity

Fig. 20. Mathematical dimensionless sensitivity of model parameters in reactive solute transport simulations in order from least (top) to most (bottom) sensitive.

unphysically great changes of quartz fraction and stabilizes the simulation. The new model was used to simulate coupled thermohaline groundwater flow and kinetically-controlled reactive transport in porous and fractured porous media. Results of short-term simulations (3 days) demonstrate the effects of temperature and salinity on silica concentrations. The

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

(a)

(b)

6

6 fr

5

high Aqz

4

original Aqz

3

fr Aqz

Si [mmol kg -1 ]

Si [mmol kg-1]

767

fr

low

2 1

5

high Aqz

4

original Aqz

3

low Aqz

2 1

0

0 0

1

2

3

0

1

distance along fracture [cm]

2

3

distance along fracture [cm]

(c) 6

6

5

5

(d) high temperature

φqz

original temperature

Si [mmol kg-1]

Si [mmol kg -1 ]

original 4 low

φqz

3 2

4

low temperature

3 2 1

1

0

0 0

1

2

3

distance along fracture [cm]

0

1

2

3

distance along fracture [cm]

Fig. 21. Visual sensitivity input parameters at steady state. Shown is the steady state quartz concentration in the fracture if the following parameters are uncertain: (a) specific quartz surface area in the fracture and (b) in the matrix, (c) quartz volume fraction and (d) temperature.

results indicate that silica concentration is inversely proportional to salinity in high-salinity regions and proportional to temperature in low-salinity regions. However, the results also show that chemical reactions do not significantly impact density-driven flow and that salt mass fluxes remain unchanged. As in Freedman and Ibaraki [21], the solubility of the solid phase (here quartz) is fairly low and, therefore, coupling between flow and reactive transport is weak when simulating over 3 days. However, the model is designed to also simulate on a larger temporal scale and an additional long-term simulation (100 years) in a fractured porous medium was conducted. The results demonstrate that, after 100 years, fracture apertures increase fivefold in the hot dissolution zone and decrease by 20% in the less hot precipitation zone. This self-sealing effect of fractures significantly reduces salt mass fluxes after about 10 years, which cannot be observed on a small temporal scale. Although fracture apertures increase substantially, reduction of fracture mass flux is attributed to the aperture minimum, which represents the bottleneck for upwelling fluid. Ultimately, at 100 years, fracture–sealing diminishes the fracture mass flux by 25% relative to the case where reactions and aperture changes are neglected. In a sensitivity analysis, parameter uncertainties were investigated. We used the simulation of reactive transport

in fractured porous media presented by Steefel and Lichtner [66] as the base case. It was found that the fracture surface area has a negligible impact on the result because the reaction rate in the matrix is about one order of magnitude faster than in the fracture and therefore dominates the reaction. The matrix surface area and the quartz volume fraction are almost equally sensitive. Both are directly proportional to the reaction rate. Variations of temperature have the greatest influence on the results because temperature impacts both the reaction kinetics and the quartz solubility. Thus, field measurements of temperature must be very accurate. The presented results of thermohaline flow and reactive transport simulations are numerically stable and obtained from the developed and fully verified model. However, the model has not been validated because appropriate field data are currently lacking. The used model input is mostly hypothetical although input parameters are representative of the type of geological material considered, and it contains some uncertainties. For example, exact fracture locations are hard to measure in the field but were shown to have a crucial impact on plume migration. Another uncertainty is the mineral surface area in the matrix and especially in the fracture. The fracture roughness is difficult to determine, such that prior studies have commonly assumed perfectly smooth fracture surfaces for simplicity. Moreover, this

768

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

model assumes that matrix permeability is only a function of absolute porosity (Eq. (24)) and ignores the effects of pore structure and connectivity. However, the latter effects may be controlled in various ways by dissolution/precipitation reactions [57]. As a consequence of uncertain model input, the 2D results presented thus far only allow an analysis and interpretation in a conceptual way. The necessary step from obtaining conceptual 2D results to making reliable 3D long-term predictions will involve an iterative cycle of further model development – sensitivity analysis – data gathering – numerical modeling – model development as proposed by Glynn and Plummer [29]. Completing this cycle, however, is highly challenging because ‘‘there are relatively few studies that have used 3D geochemical transport codes’’ [29]. Prior research that would help complete the cycle described above is rare

and it remains greatly demanding to simulate a long-term 3D thermohaline flow – reactive transport feedback scenario in a numerically stable fashion. The complexity of nature and the need to find secure deep repositories for hazardous waste are the reasons why exploring the coupled system of thermohaline flow Table A.1 Water chemistry at different depths in the Canadian Shield; all concentrations are in mg l1 [17] Solute

Freshwater 0m

Brackish water 500 m

Saltwater 1000 m

Dense brine 1500 m

Na+ Mg2+ Ca2+ Cl

9 2 15 30

360 90 630 730

3550 95 7600 24,000

34,000 25 60,000 150,000

3 viscosity [10-3 kgm-1 sec-1]

1200

fluid density [kg m-3]

1100 1000 900 800 Pitzer’smodel

700

2

1

0

600 0

100

200

300

0

200

300

temperature [ C] -5

-1

-6

log k + [mol m-2 sec-1 ]

-1.5

log Keq [mol kg-1]

100

o

o

temperature [ C]

-2 -2.5 -3 -3.5 all fluids

-4 -4.5

-7 -8 -9 -10 -11 -12 -13 -14

-5 0

100

200

300

-15 0.001

0.002

0.003

0.004

-1

temperature [oC]

1/T [K ]

activity coefficient [--]

3 2.5

dense brine

2

saltwater

1.5

brackish water

1

freshwater

0.5 0 0

100

200

300

o temperature [ C]

Fig. A.1. All physicochemical parameters calculated by FRAC3DVS are functions of temperature and salinity.

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

and reactive transport ‘‘will be an area of ongoing research’’ [50]. The model developed here was shown to be a very useful tool that can advance this future research. Acknowledgements We thank the Canadian Water Network (CWN) as well as the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support of this project. We are grateful to P.M. Dove for providing valuable information on quartz dissolution kinetics, C.J. Neville for providing details of the Meyer [42] solution and Christophe Monnin for providing the subroutine for calculating density in this code. We also thank eight anonymous AWR reviewers whose constructive comments have helped improve the manuscript. Appendix A. Parameter dependency on temperature and salinity The model computes every physicochemical system parameter over the low-temperature range 0–300 C and a wide range of salinity, shown in Table A.1. Fig. A.1 illustrates the corresponding model parameters. References [1] Bear J. Dynamics of fluids in porous media. New York: Elsevier; 1988. 764 pp. [2] Bear J, Verruijt A. Modeling groundwater flow and pollution. Dordrecht: Reidel Publishing Company; 1987. 414 pp. [3] Bennett PC. Quartz dissolution in organic-rich aqueous systems. Geochim Cosmochim Acta 1991;55:1781–97. [4] Berkowitz B, Bear J, Braester C. Continuum models for contaminant transport in fractured porous formations. Water Resour Res 1988;24(8):1225–36. [5] Blum A, Lasaga AC. Role of surface speciation in the lowtemperature dissolution of minerals. Nature 1988;331:431–3. [6] Bolton EW, Lasaga AC, Rye DM. 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. J Geophys Res 1996;101(B10):22157–87. [7] Brady PV, Walther JV. Controls on silicate dissolution rates in neutral and basic pH solutions at 25 C. Geochim Cosmochim Acta 1989;53:2823–30. [8] Brandt A, Fernando HJS. Double-diffusive convection. Geophysical monograph, vol. 94. Washington (DC): American Geophysical Union; 1995. 334 pp. [9] Dana JD. On the decay of quartzyte, and the formation of sand, kaolin and crystallized quartz. Am J Sci 1884;28:448–52. [10] Davison CC, Chan T, Brown A. The disposal of Canada’s nuclear fuel waste: site screening and site evaluation technology. Atomic Energy of Canada Limited research (AECL-10719). Pinawa (MB): Whiteshell Laboratories; 1994. 255 pp. [11] Davison CC, Chan T, Brown A. The disposal of Canada’s nuclear fuel waste: the geosphere model for postclosure assessment. Atomic Energy of Canada Limited research (AECL-10713). Pinawa (MB): Whiteshell Laboratories; 1994. 497 pp. [12] Domenico PA, Schwartz FW. Physical and chemical hydrogeology. New York: John Wiley; 1998. 506 pp. [13] Dove PM. The dissolution kinetics of quartz in aqueous mixed cation solutions. Geochim Cosmochim Acta 1999;63(22):3715–27.

769

[14] Dove PM, Crerar DA. Kinetics of quartz dissolution in electrolyte solutions using a hydrothermal mixed flow reactor. Geochim Cosmochim Acta 1990;54:955–69. [15] Dove PM, Nix CJ. The influence of the alkaline earth cations, magnesium, calcium, and barium on the dissolution kinetics of quartz. Geochim Cosmochim Acta 1997;61(16):3329–40. [16] Evans GE, Nunn JA. Free thermohaline convection in sediments surrounding a salt column. J Geophys Res 1989;94:2707–16. [17] Farvolden RN, Pfannkuch O, Pearson R, Fritz P. The Precambrian shield. In: The Geological Society of America, editor. The geology of North America. Hydrogeology 1988;O-2:101–14 [chapter 15]. [18] Fournier RO. Self-sealing and brecciation resulting from quartz deposition within hydrothermal systems. In: Proceedings of the fourth international symposium on water rock interaction, Misasa, Japan, 1983, p. 137–40. [19] Fournier RO. A method of calculating quartz solubilities in aqueous sodium chloride solutions. Geochim Cosmochim Acta 1983;47:579–86. [20] Fournier RO. The behaviour of silica in hydrothermal solutions. In: Gerger BR, Bethke PM, editors. Geology and Geochemistry of Epithermal Systems. Rev Eco Geol 1985;2:45–61. [21] Freedman V, Ibaraki M. Effects of chemical reactions on densitydependent fluid flow: on the numerical formulation and the development of instabilities. Adv Water Resour 2002;25(4):439–53. [22] Freedman V, Ibaraki M. Coupled reactive mass transport and fluid flow: issues in model verification. Adv Water Resour 2003;26(1): 117–27. [23] Freeze RA, Cherry JA. Groundwater. Englewood Cliffs (NJ): Prentice Hall; 1979. 604 pp. [24] Frind EO. Simulation of long-term transient density-dependent transport in groundwater. Adv Water Resour 1982;5(2):73–88. [25] Garven G, Freeze RA. Theoretical analysis of the role of groundwater flow in the genesis of stratabound ore deposits 1. Mathematical and numerical model 2. Quantitative results. Am J Sci 1984;284(12): 1085–174. [26] Garven G, Appold MS, Toptygina VI, Hazlett TJ. Hydrogeologic modeling of the genesis of carbonate-hosted lead–zinc ores. Hydrogeol J 1999;7:108–26. [27] Geiger S, Haggerty R, Dilles JH, Reed MH, Mattha¨i SK. New insights from reactive transport modelling: the formation of the sericitic vein envelopes during early hydrothermal alteration at Butte, Montana. Geofluids 2002;2:185–201. [28] Ghogomu NF, Therrien R. Reactive mass transport modeling in discretely-fractured porous media. In: Bentley LR et al., editors. Computational methods in water resources, vol. XIII. Netherlands: Rotterdam; 2000, ISBN 90-5809-123-6. p. 285–92. [29] Glynn PD, Plummer LN. Geochemistry and the understanding of ground-water systems. Hydrogeol J 2005;13:263–87. [30] Graf T, Therrien R. Variable-density groundwater flow and solute transport in porous media containing nonuniform discrete fractures. Adv Water Resour 2005;28(12):1351–67. [31] Graf T, Therrien R. Variable-density groundwater flow and solute transport in irregular 2D fracture networks. Adv Water Resour 2006. doi:10.1016/j.advwatres.2006.05.003. [32] Holzbecher EO. Modeling density-driven flow in porous media. Berlin: Springer; 1998. 286 pp. [33] Ibaraki M. A robust and efficient numerical model for analyses of density-dependent flow in porous media. J Contaminant Hydrol 1998;34(10):235–46. [34] Johnson JW, Knauss KG, Glassley WE, DeLoach LD, Tompson AFB. Reactive transport modeling of plug-flow reactor experiments: quartz and tuff dissolution at 240 C. J Hydrol 1998;209(10):81–111. [35] Krauskopf KB, Bird DK. Introduction to geochemistry. New York: McGraw-Hill; 1995. 647 pp. [36] Langmuir D. Aqueous environmental geochemistry. Upper Saddle River (NJ): Prentice Hall; 1997. 600 pp. [37] Lasaga AC. Chemical kinetics of water–rock interactions. J Geophys Res 1984;89(B6):4009–25.

770

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771

[38] Lichtner PC. The quasi-stationary state approximation to coupled mass transport and fluid–rock interaction in a porous medium. Geochim Cosmochim Acta 1988;52:143–65. [39] Marcus Y. Ion solvation. Chichester: John Wiley; 1985. 306 pp. [40] Marshall WL, Chen TAC. Amorphous silica solubilities V. Predictions of solubility behavior in aqueous mixed electrolyte solutions to 300 C. Geochim Cosmochim Acta 1982;46(2):289–91. [41] Mayer KU, Frind EO, Blowes DW. Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water Resour Res 2002;38(9):1174. doi:10.1029/2001WR000862. [42] Meyer JR. Development of a heat transport analytical model for a single fracture in a porous media. Unpublished project report of the course Earth 661. In: Sudicky EA, Neville CJ, editors. Analytical solutions in hydrogeology. Waterloo Center for Groundwater Research. 2004. 46 pp. [43] Molson JWH, Frind EO, Palmer C. Thermal energy storage in an unconfined aquifer 2: model development, validation and application. Water Resour Res 1992;28(10):2857–67. [44] Monnin C. An ion interaction model for the volumetric properties of natural waters: density of the solution and partial molal volumes of electrolytes to high concentrations at 25 C. Geochim Cosmochim Acta 1989;53:1177–88. [45] Monnin C. Density calculation and concentration scale conversions for natural waters. Comput Geosci 1994;20(10):1435–45. [46] Mroczek EK, Christenson B. Solubility of quartz in hypersaline brine – implication for fracture permeability at the brittle–ductile transition. In: Proceedings world geothermal congress, Kyushu-Tohoku, Japan. 2000. p. 1459–62. [47] Myers J, Ulmer GC, Grandstaff DE, Brozdowski R, Danielson MJ, Koski OH. Developments in the monitoring and control of Eh and pH conditions in hydrothermal experiments. In: Barney GS, Navratil JD, Schulz WW, editors. Geochemical behavior of disposed radioactive waste. Washington (DC): American Chemical Society; 1984. p. 97–216. [48] Nield DA, Bejan A. Convection in porous media. New York: Springer; 1999. 546 pp. [49] Oldenburg CM, Pruess K. Layered thermohaline convection in hypersaline geothermal systems. Transport Porous Media 1998;33:29–63. [50] Post VEA. Fresh and saline groundwater interaction in coastal aquifers: is our technology ready for the problems ahead? Hydrogeol J 2005;13:120–3. [51] Pruess K, Oldenburg C, Moridis G. TOUGH2 user’s guide version 2.0. Report LBNL-43134. Berkeley (CA), USA: Lawrence Berkeley National Laboratory; 1999. 197 pp. [52] Rimstidt JD. Quartz solubility at low temperatures. Geochim Cosmochim Acta 1997;61(13):2553–8. [53] Rimstidt JD, Barnes HL. The kinetics of silica–water reactions. Geochim Cosmochim Acta 1980;44:1683–99. [54] Rimstidt JD, Dove JD. Mineral/solution reaction rates in a mixed flow reactor: Wollastonite hydrolysis. Geochim Cosmochim Acta 1986;50:2509–16. [55] Saaltink MW, Carrera J, Ayora C. On the behavior of approaches to simulate reactive transport. J Contaminant Hydrol 2001;48(4):213–35. [56] Sanford WE, Konikow LF. Simulation of calcite dissolution and porosity changes in saltwater mixing zones in coastal aquifers. Water Resour Res 1989;25(4):655–67. [57] Saripalli KP, Meyer PD, Bacon DH, Freedman VL. Changes in hydrologic properties of aquifer media due to chemical reactions: a review. Crit Rev Environmental Sci Technol 2001;31(4):311–49. [58] Scha¨fer W, Therrien R. Simulating transport and removal of xylene during remediation of a sandy aquifer. J Contaminant Hydrol 1995;19(9):205–36. [59] Scha¨fer D, Scha¨fer W, Kinzelbach W. Simulating of reactive processes related to biodegradation in aquifers 1. Structure of the three-dimensional reactive transport model. J Contaminant Hydrol 1998;31(5):167–86.

[60] Schincariol RA, Schwartz FW. An experimental investigation of variable-density flow and mixing in homogeneous and heterogeneous media. Water Resour Res 1990;26(10):2317–29. [61] Schincariol RA, Schwartz FW, Mendoza CA. On the generation of instabilities in variable-density flows. Water Resour Res 1994;30(4):913–27. [62] Shibue Y. An empirical equation for quartz solubility in NaCl solution. J Mineral Petrol Eco Geol 1994;89:203–12. [63] Shikaze SG, Sudicky EA, Mendoza CA. Simulations of dense vapour migration in discretely-fractured geologic media. Water Resour Res 1994;30(7):1993–2009. [64] Shikaze SG, Sudicky EA, Schwartz FW. Density-dependent solute transport in discretely-fractured geologic media: is prediction possible? J Contaminant Hydrol 1998;34(10):273–91. [65] Steefel CI, Lasaga AC. A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. Am J Sci 1994;294(5):529–92. [66] Steefel CI, Lichtner PC. Multicomponent reactive transport in discrete fractures: I. Controls on reaction front geometry. J Hydrol 1998;209:186–99. [67] Steefel CI, Lichtner PC. Multicomponent reactive transport in discrete fractures: II. Infiltration of hyperalkaline groundwater at Maqarin, Jordan, a natural analogue site. J Hydrol 1998;209:200–24. [68] Steefel CI, MacQuarrie KTB. Approaches to modeling of reactive transport in porous media. In: Lichtner PC, Steefel CI, Oelkers EH, editors. Reviews in mineralogy, vol. 34. Washington (DC): Mineralogical Society of America; 1996. p. 3–129 [chapter 2]. [69] Steefel CI, Yabusaki SB. OS3D/GIMRT: Software for modeling multicomponent-multidimensional reactive transport; user manual and programmer’s guide. Technical Report PNL-1116. Richland (WA), USA: Pacific Northwest National Laboratory; 1996. 56 pp. [70] Stern ME. The ‘‘salt-fountain’’ and thermohaline convection. Tellus 1960;12(1):172–5. [71] Stober I, Bucher K. Deep fluids: neptune meets pluto. Hydrogeol J 2005;13:112–5. [72] Stumm W, Morgan JJ. Aquatic chemistry. New York: John Wiley; 1996. [73] Sudicky EA, McLaren RG. The Laplace transform Galerkin technique for large-scale simulation of mass transport in discretelyfractured porous formations. Water Resour Res 1992;28(2):499–514. [74] Tang DH, Frind EO, Sudicky EA. Contaminant transport in fractured porous media: analytical solution for a single fracture. Water Resour Res 1981;17(3):555–64. [75] Tester JW, Worley WG, Robinson BA, Grigsbay CO, Feerer JL. Correlating quartz dissolution kinetics in pure water from 25 to 625 C. Geochim Cosmochim Acta 1994;58(11):2407–20. [76] Therrien R, Sudicky EA. Three-dimensional analysis of variably saturated flow and solute transport in discretely-fractured porous media. J Contaminant Hydrol 1996;23(6):1–44. [77] Therrien R, McLaren RG, Sudicky EA, Panday SM. HYDROSPHERE – a three-dimensional numerical model describing fullyintegrated subsurface and surface flow and solute transport. Universite´ Laval, University of Waterloo; 2004. 275 pp. [78] Turner JS. Buoyancy effects in fluids. Cambridge: Cambridge University Press; 1979. 368 pp. [79] Tyvand PA. Thermohaline instability in anisotropic porous media. Water Resour Res 1980;16:325–30. [80] von Damm KL, Bischoff JL, Rosenbauer RJ. Quartz solubility in hydrothermal seawater: an experimental study and equation describing quartz solubility for up to 0.5 m NaCl solutions. Am J Sci 1991;291(12):977–1007. [81] Walter AL, Frind EO, Blowes DW, Ptacek CJ, Molson JWH. Modelling of multicomponent reactive transport in groundwater: 2. Metal mobility in aquifers impacted by acidic mine tailings discharge. Water Resour Res 1994;30:3149–58. [82] Weir GJ, White SP. Surface deposition from fluid flow in a porous medium. Transport Por Media 1996;25:79–96.

T. Graf, R. Therrien / Advances in Water Resources 30 (2007) 742–771 [83] White SP, Mroczek EK. Permeability changes during the evolution of a geothermal field due to the dissolution and precipitation of quartz. Transport Por Media 1998;33:81–101. [84] Yang J, Edwards RN. Predicted groundwater circulation in fractured and unfractured anisotropic porous media driven by nuclear fuel waste heat generation. Can J Earth Sci 2000;37:1301–8. [85] Yeh GT, Tripathi VS. A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resour Res 1989;25(1):93–108.

771

[86] Yoshida J, Nagashima H, Nagasaka M. Numerical experiment on double diffusive currents. In: Brandt A, Fernando HJS, editors. Double-diffusive convection. Washington (DC): American Geophysical Union; 1995. p. 69–79. [87] Zheng C, Bennett GD. Applied contaminant transport modeling. New York: John Wiley; 2002. 621 pp. [88] Zysset A, Stauffer F, Dracos T. Modeling of chemically reactive groundwater transport. Water Resour Res 1994;30(7): 2217–28.