Reactive transport modeling of the interaction between a high-pH plume and a fractured marl: the case of Wellenberg

Applied Geochemistry 18 (2003) 1555–1571 www.elsevier.com/locate/apgeochem

Reactive transport modeling of the interaction between a high-pH plume and a fractured marl: the case of Wellenberg Josep M. Soler1 Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland Received 7 June 2002; accepted 10 February 2003 Editorial handling by R.L. Bassett

Abstract In the context of the proposed low- and intermediate-level radioactive waste repository at Wellenberg (Switzerland), calculations simulating the interaction between hyperalkaline solutions and a fractured marl, at 25
C, have been performed. The aim of these calculations is to evaluate the possible effects of mineral dissolution and precipitation on porosity and permeability changes in such a fractured marl, and their impact on repository performance. Solute transport and chemical reaction are considered in both a high-permeability zone (fracture), where advection is important, and the wall rock, where diffusion is the dominant transport mechanism. The mineral reactions are promoted by the interaction between hyperalkaline solutions derived from the degradation of cement (a major component of the engineered barrier system in the repository) and the host rock. Both diffusive/dispersive and advective solute transport are taken into account in the calculations. Mineral reactions are described by kinetic rate laws. The fluid flow system under consideration is a two-dimensional porous medium (marl, 1% porosity), with a high-permeability zone simulating a fracture (10% porosity) crossing the domain. The dimensions of the domain are 6 m per 1 m, and the fracture width is 10 cm. The fluid flow field is updated during the course of the simulations. Permeabilities are updated according to Kozeny’s equation. The composition of the solutions entering the domain is derived from modeling studies of the degradation of cement under the conditions at the proposed underground repository at Wellenberg. Two different cases have been considered in the calculations. These 2 cases are representative of 2 different stages in the process of degradation of cement (pH 13.5 and pH 12.5). In both cases, the flow velocity in the fracture diminishes with time, due to a decrease in porosity. This decrease in porosity is caused by the precipitation of calcite (replacement of dolomite by calcite) and other secondary minerals (brucite, sepiolite, analcime, natrolite, tobermorite). However, the decrease in porosity and flow velocity is much more pronounced in the lower pH case. The extent of the zone of mineral alteration along the fracture is also much more limited in the lower pH case. The reduction of porosity in the fractures would be highly beneficial for repository performance, since it would mean that the solutions coming from the repository and potentially carrying radionuclides in solution would have to flow through low-conductive rock before they would be able to get to higher-conductive features. The biggest uncertainty in the reaction rates used in the calculations arises from the surface areas of the primary minerals. Additional calculations making use of smaller surface areas have also been performed. The results show that the smaller surface areas (and therefore smaller reaction rates) result in a smaller reactivity of the system and smaller porosity changes. # 2003 Elsevier Ltd. All rights reserved.

1

E-mail address: [email protected] (J.M. Soler). Current address: Institut de Cie`ncies de la Terra ‘‘Jaume Almera’’ (CSIC), Lluı´s Sole´ i Sabarı´s s/n, 08028 Barcelona, Spain.

0883-2927/03/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0883-2927(03)00048-9

1556

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

1. Introduction A project for a repository for low- and intermediatelevel radioactive waste at Wellenberg (Canton Nidwalden, Switzerland) is currently under consideration. Cement is a major component of the engineered barrier system in such a repository, and therefore, the interaction between the hyperalkaline solutions derived from the degradation of cement and the host rock (marl) may change the physical and chemical properties of the rock in the vicinity of the repository. The objective of this study is to provide some insight into the consequences of such interaction (changes in mineralogy, porosity, permeability), based on reactive transport calculations. Unlike some previous investigations, this study accounts for the important feedback between changes in porosity induced by mineral reaction and changes in fluid flow. This is purely a numerical modeling exercise, where solute transport and chemical reaction, driven by the flux of the hyperalkaline solutions, are considered in both a high-permeability zone (fracture), where advection is important, and the wall rock, where diffusion is the dominant transport mechanism. The calculations are based on the available mineralogical and structural information at the site, available thermodynamic and kinetic data for the system of interest (extrapolation of reaction rates to high-pH conditions has been required) and available information and estimates of the transport parameters. Important degrees of uncertainty are associated with some of the parameters used in the calculations. Also, the results cannot be compared to a real system (laboratory or field experiment) and therefore lack confirmation. However, the overall conclusions are consistent with other experimental and modeling studies of the interaction between hyperalkaline solutions and different types of rocks, i.e., the progressive sealing of porosity (Steefel and Lichtner, 1994; Lichtner and Eikenberg, 1995; Soler, 1998, 2000; Adler, 2001; Soler and Ma¨der, 2002; Savage et al., 2002). The results should be taken as an indication of the possible evolution of the system given the current knowledge and conceptual model of the relevant chemical and transport processes in the vicinity of the repository.

2. The model 2.1. Reactive transport calculations: The GIMRT program The GIMRT program (Steefel and Yabusaki, 1996) is part of a numerical software package for simulating multicomponent reactive transport in porous media. GIMRT (Global Implicit Multicomponent Reactive Transport) numerically solves the reaction-transport differential equations in either one or two dimensions,

making use of a one-step or global implicit approach, i.e., solving the transport and reaction terms simultaneously (Steefel and Lasaga, 1994). The reaction-transport equations, which describe the conservation of solute mass, are given by (Steefel, 1992; Steefel and Lasaga, 1994)  @
f MH2 O Ci þ r @t

  u~f MH2 O Ci  Dr f MH2 O Ci ¼ Ri ði ¼ 1; 2; :::Ntot Þ

ð1Þ

where Ci is the molal concentration of a species in solution, MH2 O is the mass fraction of H2 O,  is porosity, f is the fluid density, u~ is the Darcy velocity, D is the combined dispersion-diffusion tensor, Ri is the total reaction rate of species i in solution (in units of moles per unit volume rock per unit time), and Ntot is the total number of aqueous species. Porosity, , is embedded in both the flux terms and the reactions terms, as described below. The Darcy velocity is related to the true velocity of the fluid, v~, by u~ ¼ ~v

ð2Þ

and D is defined as the sum of the mechanical or kine matic dispersion, D , and the molecular diffusion coefficient in water, D0 , divided by the formation resistivity factor, F

D¼D þ

D0 F

The kinematic dispersion is written as   ui uj Dij ¼ T u~ þ ðL  T Þ   u~

ð3Þ

ð4Þ

where L and T are the longitudinal and transverse dispersivities, respectively (Bear, 1979; de Marsily, 1986). The model assumes that the principal direction of flow is aligned with the grid, i.e., Dij is a diagonal matrix. The definition of F used in the model is based on Archie’s Law and is given by F ¼ m

ð5Þ

where m (the ‘‘cementation exponent’’) has reported values between 1.3 and 5.4, and values around 2 for deeply buried compacted sediments (Dullien, 1979; Horseman et al., 1996; Ullman and Aller, 1982). The reaction term Ri results from the sum of all the individual mineral-water reactions which affect the concentration of the ith species Nm X Ri ¼  im rm ð6Þ m¼1

1557

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

where rm is the net rate of precipitation (rm > 0) or dissolution (rm < 0) of mineral per unit volume of rock, im is the number of moles of species i per mole of mineral m, and Nm is the number of minerals present in the rock. In order to describe how the reaction rate term, Ri , is treated in GIMRT, it is useful at this moment to introduce the concepts of primary and secondary species. The total number of species in solution (Ntot ), can be distributed between Nc primary species and Nx secondary species. The number of primary species (Nc ), which is the number of independent chemical components in the system, will be given by the total number of species (Ntot ) minus the number of independent chemical equilibria relating them. Following the approach previously described by Reed (1982), Lichtner (1985), Kirkner and Reeves (1988), and Steefel (1992), the different equilibrium relationships can be written as the production or destruction of 1 mol of secondary species, and take the form Nc X Ai () ij Aj ði ¼ 1; 2; :::Nx Þ

ð7Þ

j¼1

where Aj and Ai are the chemical formulas of the primary and secondary species, respectively, and ij is the number of moles of primary species j in 1 mol of secondary species i. The chemical equilibria provide an algebraic link between the primary and secondary species via the law of mass action for each reaction. The concentration of a secondary species i is given by Xi ¼

Ki1 i1

Nc Y


j Cj

ij

Cj þ

Nx X ij Xi i¼1

defines the total concentration of species j, Uj (Reed, 1982; Lichtner, 1985; Kirkner and Reeves, 1988). Also, in this formulation of the transport-reaction equation, it is assumed that the diffusion coefficients are the same for all the aqueous species (both primary and secondary). The mineral reaction rate laws implemented in GIMRT are of the form ! NY c þNx p rm ¼ sgnðlog½Qm =Km ÞAm km ai fðDGÞ ð11Þ i¼1

where sgn means sign, Qm is the ion activity product, Km is the equilibrium constant of the mineral dissolution reaction (ionic activity product at equilibrium), Am is the mineral surface area per unit volume of rock, km is the growth or dissolution rate constant (in units of moles of mineral per unit surface area per time), ai is the activity of an inhibiting or catalyzing species, raised to an empirically determined power p, and fðDGÞ is a function describing the dependence of the rate on solution saturation state. This function guarantees that the rate is zero at equilibrium; fðDG ¼ 0Þ ¼ 0

ði ¼ 1; 2; :::Nx Þ

ð8Þ The ion activity product is defined as

j¼1

where the Ki are the equilibrium constants of reaction (8), and j and i are the activity coefficients of the primary and secondary species, respectively. The activity coefficients can either be arbitrarily chosen to be all equal to unity, or calculated according to the extended Debye–Hu¨ckel formulation, given by logi ¼ 

The term

: AZi2 I 1=2 þ bI : 1 þ Bai I 1=2

ð9Þ

Z is the ionic charge, I is: the ionic strength, and the : parameters A, B, ai , and b are from the EQ3 database (Wolery et al., 1990). The resulting transport-reaction equations are of the form " !# Nx X @ f MH2 O Cj þ ij Xi þr @t i¼1 " ! ( !)# Nx Nx X X  u~f MH2 O Cj þ ij Xi Dr f MH2 O Cj þ ij Xi i¼1

i¼1

¼ Rj ðj ¼ 1; 2; :::Nc Þ ð10Þ

Qm ¼

Nc Y  aj mj

ð12Þ

j¼1

where aj are the activities of the species making up the mineral m, and mj are the stoichiometric coefficients. The effect of solution saturation state is given by the last term of Eq. (11), which can be expressed in terms of the Gibbs Free Energy of the reaction ðDGÞ  n  M n

   Qm DG   1 ð13Þ fðDGÞ ¼ exp M  1 ¼    Km RT M and n are two positive numbers which are normally determined experimentally. If no experimental information is available, it is common practice to make both M and n equal to unity, which reflects the dependence of the net rate of an elementary reaction on solution saturation state, as derived from Transition State Theory (Lasaga, 1998). In this study the original code was modified so the program could also accept a form of fðDGÞ given by   M   Q m M DG  fðDGÞ ¼   ¼ ln ð14Þ Km  RT

1558

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

ð15Þ

The two-dimensional flow domain has a rectangular shape (Fig. 1). The numerical discretization of the domain in the model is based on a cell-centered finite difference scheme. Hydraulic heads are constant along the left- and right-hand-side boundaries, defining an overall hydraulic gradient from left to right. The top and bottom boundaries are no-flow boundaries. The array of Darcy velocities is calculated by solving the equation of conservation of fluid mass, assuming a constant density of the fluid. The initial flow field is calculated by solving the equation of conservation of fluid mass at steady state, which is given by

ð16Þ

r u~ ¼ 0

which has been found to apply to gibbsite and kaolinite precipitation rates (Nagy et al., 1991; Nagy and Lasaga, 1992). Also, the program was modified to allow for different dissolution and precipitation rate laws [Eq. (11)] for one same mineral. In the modified version of GIMRT used in this study, and unlike the original version, mineral grains are assumed to be spherical, with their whole surface area in contact with fluid (‘‘floating’’ grains). Mineral volume fractions and surface areas are calculated according to 4 m ¼ Rm3 Nm 3 Am ¼ 4 Rm2 Nm

where Rm is the grain radius, and Nm is the number of grains per volume of rock (constant for each mineral). Porosity, which is updated along the calculation, is given by ¼1

Nn X

@ @t

ð21Þ

ð17Þ

m¼1

Mineral volume fractions are updated by making use of the calculated reaction rates and the duration of the time increment for each time step along the simulation (amounts of mineral dissolved or precipitated at each time step). Grain radii and mineral surface areas are updated according to Eqs. (15) and (16). Notice that the spherical grain geometry is only used for the calculation of mineral surface areas and is not taken into account in the calculation of porosity. The discretization of the reaction-transport differential equations [Eq. (10)] is done following the integrated finite difference method (Patankar, 1980; de Marsily, 1986), using a fully implicit approach for the time derivatives. The resulting algebraic non-linear equations are solved using Newton’s method. 2.2. Flow calculations A two-dimensional flow program has been added to the GIMRT reactive transport code to calculate the initial fluid flow field and its evolution during the simulations. The flow field used in the reactive transport simulations consists of an array of Darcy velocities (~ u), which are given by u~ ¼ Krh

At later times, the changes of porosity and permeability of the system are taken into account, and the flow field is updated using r u~ ¼ 

m

ð20Þ

ð18Þ

where K is the hydraulic conductivity of the medium, and h is the hydraulic head. The porous medium is assumed to be isotropic (K is a scalar), with K given by kf g K¼ ð19Þ  where k is the permeability, g is the gravitational acceleration, and  is the viscosity of the fluid.

Both Eqs. (20) and (21) are solved for the hydraulic heads by means of a cell-centered finite difference discretization scheme. The system of algebraic equations resulting from the discretization is solved using the Successive Over Relaxation Technique, and once the hydraulic heads are known, Darcy velocities are calculated from Eq. (18). The functional dependence of permeability on porosity, which is used for the update of the flow field, has been assumed to follow Kozeny’s equation (see for instance Bear, 1972), which has the form k ¼ c0

3 S2

ð22Þ

where S is the total mineral surface area per volume of rock, and c0 is a constant. The update of the flow field is performed by numerically solving Eq. (21) when the porosity in the domain has changed by a given amount. At each time step, the relative porosity change ðjnew  old j=old Þ at each point in the flow domain is calculated, and the maximum change is recorded. When the summation over successive time steps of the maximum relative porosity changes is equal or greater than a specified amount (0.1, equivalent to 10%), the flow field is updated.

3. Parameters The parameters used in the calculations are intended to simulate the conditions in the host rock of the planned repository for low- and intermediate-level radioactive waste at Wellenberg, about 18 km SE of Luzern (Nagra, 1994, 1997). These host rocks are the marls of the Palfris Formation and Vitznau Marls, which were deposited in the Early Cretaceous (Valanginian). The adjacent Globigerina Marls and Schimberg Shales, deposited in the Tertiary, display properties similar to

1559

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

Fig. 1. Schematic diagram showing the geometry and boundary conditions of the flow domain used in the model. Left- and right-hand-side boundaries are constant head boundaries. Top and bottom boundaries are no-flow boundaries. Lx=6 m. Ly=1 m. Fracture width is 10 cm.

the Early Cretaceous marls and have also been considered suitable for hosting a repository. The main water-conducting features in the host rock are cataclastic shear zones, which display total thicknesses ranging from decimeters up to a few meters. These cataclastic shear zones contain fractures filled with fault gouge (highly porous ground rock), with thicknesses of the order of millimeters. The porosity of the fault gouge ranges from 5% to 20%. A zone of strongly deformed rock with thicknesses of the order of several centimeters and porosities of the order of 3% can be found on both sides of the fault gouges. The undeformed marls display porosities of the order of 1–2%. A simplified structural model has been used in the calculations. The domain of the simulations is a 2-dimensional section of a porous medium (marl, 1% porosity), crossed by a high-permeability zone simulating a fracture (10% porosity). The fracture crosses the whole domain, from the left- to the right-hand-side boundary, which is equivalent to the assumption of a worst case scenario, with a high-conductive zone starting right at the repository. 3.1. Dimensions

3.2. Rock composition The initial mineralogical composition and porosity of both wall rock and high-permeability zone (fracture) are given in Tables 1 and 2, respectively, and are based on Table 1 Initial mineralogical composition of the wall rocka

Calcite Dolomite Quartz Muscovite Porosity

Ly ¼ 1 m The grid spacing is 0.1 m in the x direction, and 0.05 m in the y direction. The high-permeability zone (fracture) crosses the domain from left to right parallel to the x axis, and is 10 cm thick (2 nodes along the y direction).

Surface area (m2/m3 rock)

Number of grains per m3 rock

50 6 15 28 1

1495 180 450 839 –

1.19108 1.43107 3.58107 6.68107 –

a Surface areas calculated according to a spherical grain geometry.

Table 2 Initial mineralogical composition of the high-permeability zonea

The dimensions of the domain are (see Fig. 1): Lx ¼ 6 m

Volume %

Calcite Dolomite Quartz Muscovite Porosity a

Volume %

Surface area (m2/m3 rock)

Number of grains per m3 rock

45 5 14 26 10

1345 150 420 780 –

1.07108 1.19107 3.34107 6.21107 –

Surface areas calculated according to a spherical grain geometry.

1560

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

the reported composition of the Valanginian marls at Wellenberg (Nagra, 1994, 1997). Muscovite is used as a surrogate for illite, mixed layer illite/smectite, and chlorite. An initial grain radius of 1 mm has been assumed for all the primary minerals. Surface areas and number of grains per volume of rock are then calculated from Eqs. (15) and (16). Regarding the secondary minerals, they are allowed to start precipitating when the solution reaches supersaturation ðDG > 0Þ with respect to that mineral. The initial surface area available for precipitation of secondary minerals is given by a fixed number of grains (nuclei) per volume of rock, and an initial grain radius equal to 0.01 mm. The initial number of grains for all the secondary minerals was 1011 grains/m3rock. The value of the initial surface for secondary minerals does not have any major influence on the simulation results at the time scales of interest (thousands to tens of thousands of years). However, this is not the case for the primary minerals, and the evolution of the system under consideration is dependent on this parameter. Further discussion is presented in Section 3.5 (reaction rates). 3.3. Solution composition The composition of the solutions entering the domain is derived from a modeling study of the degradation of cement under the conditions at the proposed underground repository at Wellenberg (Neall, 1994). These are highly alkaline solutions (pH 12.5–13.5), at equilibrium with calcite, and undersaturated with respect to the other primary minerals in the marl. Neall (1994) used a chemical equilibrium model, based on the cement degradation model developed by Berner (1987, 1990), to simulate the degradation of different types of cement in contact with Wellenberg groundwaters. Six different types of cement (the types to be used in the construction of the repository) and 2 types of groundwater were considered by Neall (1994). The 2 different groundwaters correspond to the two groundwaters observed in the Wellenberg area: a shallower and more dilute NaHCO3-type water, and a deeper and more saline NaCl-type water (Nagra, 1993). The compositions of these 2 groundwaters have also been used as initial conditions (the initial water in the pores of the rock) in the simulations. The results from the cement degradation model, supported by experimental evidence, show that there are 3 distinct phases in the degradation process. In the first stage, the solution composition is dominated by the dissolution of the alkali hydroxides in the cement, giving an initial pH of about 13.5. In the second stage, the pH is about 12.5 and the composition is controlled by the dissolution of the Ca(OH)2 component of the CSH (Calcium–Silicate–Hydrate) gel. The third stage is characterized by the dissolution of phases other than the CSH

gel (katoite, ettringite, secondary gibbsite) and results in significant changes in solution chemistry, with pH going down to normal groundwater values. The solutions resulting from the degradation of the different types of cement are similar to each other in composition. Only the solution compositions corresponding to the degradation of average Swiss Portland Cement have been used in the calculations. A solution representative of the first stage of the degradation of cement (pH about 13.5) and another one representative of the second stage (pH about 12.5) have been used as solutions entering the fractured marl. As a result, 4 different cases have been considered: 1. 2. 3. 4.

pH pH pH pH

13.44, 12.50, 13.48, 12.56,

NaCl-type groundwater NaCl-type groundwater NaHCO3-type groundwater NaHCO3-type groundwater

Since the results for cases 1 and 2 (NaCl-type groundwater) are very similar to the results for cases 3 and 4 (NaHCO3-type groundwater), only the results for the first 2 cases will be shown. Neall (1994) estimated the duration of the first phase of cement degradation (pH13.5) to be between 2104 and 2106 a, based on the porosity of the cement and the average hydraulic conductivity and hydraulic gradient in the host rock. The estimated duration of the second phase (pH12.5) ranged from 2106 to 2108 a. Typical time spans of performance assessment calculations, which include the periods of peak release of all the different radionuclides from repositories, are in the range of 106–107 a. 3.4. Equilibrium constants Nineteen minerals (Nm ¼ 19) and 41 species in solution (Ntot ¼ 41) have been taken into account in the calculations. All the chemical equilibria in solution are listed in Table 3, and are written according to Eq. (7) (destruction of 1 mol of secondary species). The equilibrium constants for all the mineral reactions are given in Table 4. The values of all the equilibrium constants are from the thermodynamic database included in GIMRT, which is based on the EQ3/6 database (Wolery et al., 1990). Activity coefficients are calculated according to Eq. (9). The activity of water is taken to be unity. 3.5. Reaction rates Reaction rate parameters for the different mineral reactions are given in Table 5. Reaction rates for the primary minerals (calcite, dolomite, muscovite, quartz), and for secondary gibbsite and kaolinite, are based on experimental work published in the literature. Since

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571 Table 3 List of chemical equilibria in solution (speciation reactions). The reactions are written according to Eq. (8), i.e., as the destruction of 1 mol of secondary speciesa Equilibria

log K

 Alþ3 () AlðOHÞ 4 4OH  AlOHþ2 () AlðOHÞ 4 3OH   AlðOHÞþ ( ) AlðOHÞ  2OH 2 4  AlðOHÞ3 ðaqÞ () AlðOHÞ 4 OH  H3 SiO ( ) SiO ð aq Þ þ OH þ H2 O 2 4  H2 SiO2 4 () SiO2 ðaqÞ þ 2OH  NaHCO3 ðaqÞ () Naþ þ CO2 3  OH þ 2 NaCO 3 () Na þ CO3 þ 

33.83 24.85 15.94 8.01 4.18 5.08 þ H2 O 3.51 0.51 NaOHðaqÞ () Na þ OH 0.18 NaHSiO3 ðaqÞ () Naþ þ SiO2 ðaqÞ þ OH 5.69 þ 2 0.82 NaSO 4 () Na þ SO4 KOHðaqÞ () Kþ þ OH 0.46 þ 2 KSO 0.88 4 () K þ SO4 CaH2 SiO4 ðaqÞ () Caþ2 þ SiO2 ðaqÞ þ 2OH 9.43 þ2 þ SiO2 ðaqÞ þ OH þ H2 O 5.20 CaH3 SiOþ 4 () Ca CaðH3 SiO4 Þ2 () Caþ2 þ 2SiO2 ðaqÞ þ 2OH þ 2H2 O 12.94 þ2  CaHCOþ þ CO2 2.62 3  OH þ H2 O 3 () Ca þ2 2 CaCO3 ðaqÞ () Ca þ CO3 3.33 1.15 CaOHþ () Ca2þ þ OH CaSO4 ðaqÞ () Caþ2 þ SO2 2.11 4 MgCO3 ðaqÞ () Mgþ2 þ CO2 2.98 10.51 MgH2 SiO4 ðaqÞ () Mgþ2 þ SiO2 ðaqÞ þ 2OH þ2 MgH3 SiOþ þ SiO2 ðaqÞ þ OH þ H2 O 5.45 4 () Mg þ þ2  MgHCO3 () Mg þ CO2 2.63 3  OH þ H2 O MgOHþ () Mgþ2 þ OH 2.21 2.41 MgSO4 ðaqÞ () Mgþ2 þ SO2 4 2  HCO 3.67 3 () CO3  OH þ H2 O  CO2 ðaqÞ () CO2 11.32 3  2OH þ H2 O  H2 SO4 ðaqÞ () SO2 29.01 4  2OH þ 2H2 O  2  12.02 HSO4 () SO4  OH þ H2 O Hþ ()  OH þ H2 O 14.00 a All equilibrium constants are from the database included in GIMRT, which is based on the EQ3/6 database (Wolery et al., 1990).

most experimental rate data has been obtained under pH conditions ranging from slightly acidic to slightly alkaline, extrapolation of the pH dependencies of rates under alkaline conditions to pH values about 11–12 has been necessary. Relatively fast rate constants have been used for the rest of the secondary minerals. These rate constants have been set to an arbitrary value of 109 mol/m2/s. The results at the time scales of interest (thousands to tens of thousands of years), and for the given set of kinetic parameters for the primary minerals, resemble the local equilibrium solution (simulations run with rate constants 2 orders of magnitude smaller for all secondary minerals show extremely similar results). The largest source of uncertainty for the rates of primary

1561

minerals arises from the surface area term. The degree of uncertainty in this term is of the order of a few orders of magnitude, which is larger than the uncertainty for the reaction rate constants (uncertainties of up to about an order of magnitude; White and Brantley, 1995). Surface areas could conceivable be 3 orders of magnitude larger than the ones specified in Table 1 if specific surface areas were of the order of 10 m2/g (a typical value for fine-grained granular mineral assemblages) and the whole surface area were in contact with solution. If not all the mineral surface area were in contact with solution, available reactive surface area could be much smaller. Such uncertainty has clearly a big effect on the model results that will be presented. The trend is to have larger changes of porosity (driven by faster reaction rates) the larger the surface area. However, the direction of the change (e.g., decrease in porosity and permeability at or near the fracture inlet, as will be shown) remains the same, given the fact that the rates for the secondary minerals are fast enough (conditions close to local equilibrium for the secondary minerals). Given these conditions, the general trends shown in the results should apply, although spatial and temporal scales will depend on the value of the surface areas. Examples of this dependency will be shown. No explicit dependence of the rates on pH has been used in the simulations (no catalysis by means of a apHþ term in the rate law), although the rates are assumed to apply to the prevailing alkaline conditions in the marl during the course of the simulations, as discussed above. The rate laws used for all the minerals are a simplified version of Eq. (11), with no catalytic/inhibitory effects taken into account, and are given by rm ¼ sgnðlog½Qm =Km ÞAm km fðDGÞ

ð23Þ

3.6. Transport properties The initial hydraulic conductivities of the wall rock and high-permeability zone are given in Table 6. Constant hydraulic heads on the left- and right-hand-side boundaries (h2=12.4 m; h1=10 m) provide a constant overall hydraulic gradient of 0.4 [(h2h1)/Lx=0.4]. Both the hydraulic conductivities and the value of the hydraulic gradient are based on the reported conditions at Wellenberg (Nagra, 1994, 1997). Fluid flow is initially parallel to the x axis (Fig. 1), with Darcy velocities in the wall rock and in the fracture equal to 3.311013 m/s (1.04105 m/a) and 4.01010 m/s (1.26102 m/a), respectively. The molecular diffusion coefficients, cementation exponents, and dispersivities used in the simulations are given in Table 6. They are based on common values reported in the literature (see for instance de Marsily, 1986; Domenico and Schwartz, 1990; Ullman and Aller, 1982).

1562

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

Table 4 Equilibrium constants for the different mineral reactionsa Mineral

Reaction

log K

Calcite Dolomite Muscovite Quartz Analcime Anhydrite Brucite Ettringite Foshagite Gibbsite Gypsum Hillebrandite Kaolinite Katoite Laumontite Natrolite Portlandite Sepiolite Tobermorite

CaCO3 þ Hþ () Caþ2 þ HCO 3 CaMgðCO3 Þ2 þ2Hþ () Ca2þ þ Mg2þ þ 2HCO 3 KAl3 Si3 O10 ðOHÞ2 þ10Hþ () Kþ þ 3Al3þ þ 3SiO2 ðaqÞ þ 6H2 O SiO2 () SiO2 ðaqÞ Na0:96 Al0:96 Si2:04 O6 H2 O þ 3:84Hþ () 0:96Al3þ þ 0:96Naþ þ 2:04SiO2 ðaqÞ þ 2:92H2 O CaSO4 () Caþ2 þ SO2 4 MgðOHÞ2 þ2Hþ () Mgþ2 þ 2H2 O þ2 Ca6 Al2 ðSO4 Þ3 ðOHÞ12 26H2 O þ 12Hþ () 2Alþ3 þ 3SO2 þ 38H2 O 4 þ 6Ca Ca4 Si3 O9 ðOHÞ2 0:5H2 O þ 8Hþ () 3SiO2 ðaqÞ þ 4Caþ2 þ 5:5H2 O AlðOHÞ3 þ3Hþ () Alþ3 þ 3H2 O CaSO4 2H2 O () Caþ2 þ SO2 4 þ 2H2 O Ca2 SiO3 ðOHÞ2 0:17H2 O þ 4Hþ () SiO2 ðaqÞ þ 2Caþ2 þ 3:17H2 O Al2 Si2 O5 ðOHÞ4 þ6Hþ () 2Alþ3 þ 2SiO2 ðaqÞ þ 5H2 O Ca3 Al2 ðOHÞ12 þ12Hþ () 2Alþ3 þ 3Caþ2 þ 12H2 O CaAl2 Si4 O12 4H2 O þ 8Hþ () Caþ2 þ 2Alþ3 þ 4SiO2 ðaqÞ þ 8H2 O Na2 Al2 Si3 O10 2H2 O þ 8Hþ () 2Al3þ þ 2Naþ þ 3SiO2 ðaqÞ þ 6H2 O CaðOHÞ2 þ2Hþ () Caþ2 þ 2H2 O Mg4 Si6 O15 ðOHÞ2 6H2 O þ 8Hþ () 4Mgþ2 þ 6SiO2 ðaqÞ þ 11H2 O Ca5 Si6 O16 ðOHÞ2 9:5H2 O þ 10Hþ () 5Caþ2 þ 6SiO2 ðaqÞ þ 15:5H2 O

1.85 2.51 11.02 4.00 5.32 4.31 16.30 60.81 65.92 6.80 4.48 36.82 5.10 77.23 11.96 16.81 22.56 30.44 63.84

a The first 4 reactions correspond to the primary minerals that make up the marl. The rest correspond to the secondary minerals that have been taken into account in the simulations. The reactions are written here the way they are written in the thermodynamic database included in GIMRT. During the simulations, the program rewrites the reactions in terms of primary species, and recalculates the values of the equilibrium constants accordingly.

Table 5 Parameters for the mineral reaction rate lawsa Mineral

km (mol/m2/s) M

Calcite 1.0107 Dolomite 2.0109 Muscovite (prec.) 3.01013 Muscovite (diss.) Quartz Gibbsite (prec.) Gibbsite (diss.) Kaolinite (prec.)

1.01012 1.01011 4.11012 1.01011 9.01014

Kaolinite (diss.) Other phases

3.01013 1.0109

n

References

0.5 1.0 1.0

2.0 Chou et al. (1989); Morse and Mackenzie (1990); Sjo¨berg (1978) 1.0 Chou et al. (1989); Morse and Mackenzie (1990) – Kalinowski and Schweda (1996); Knauss and Wolery (1989); Lin and Clemency (1981); Nickel (1973); Stumm et al. (1987); Nagy et al. (1991) 0.85 1.0 1.0 1.0 Dove (1994); Rimstidt and Barnes (1980) 1.1 – Bloom and Erich (1987); Mogollo´n et al. (1996); Nagy and Lasaga (1992) 7.0 9.0 1.0 – Carroll and Walther (1990); Carroll-Webb and Walther (1988); Ganor et al. (1995); Nagy et al. (1991); Wieland and Stumm (1992) 0.85 1.0 1.0 1.0 –

a For muscovite, gibbsite, and kaolinite, the precipitation and dissolution rate laws are different. For the rest of minerals, the same rate law applies to both dissolution and precipitation. Also, when both M and n are given, the dependence of the rate on solution saturation [ f ðDGÞ] state follows Eq. (15). If only M is given, the dependence follows Eq. (16).

Table 6 Initial transport properties assumed in the calculations Hydr. cond. (fracture) Hydr. cond. (wall rock) Molec. diffusion coeff. Cementation exponent Longitud. dispersivity Transverse dispersivity

K (fracture) K (wall rock) D0 m L T

109 m/s 8.271013 m/s 109 m2/s 2 0.1 m 0.01 m

4. Results and discusion Results for the 2 different cases considered (see Section 3.3) are shown below. Although the calculations are fully 2-dimensional, emphasis is placed on the evolution of system along the fracture, where the amount of mineral alteration is significant. Mineral alteration is only very minor in the wall rock, although it is more important in the high pH case (the zone of mineral

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571 Table 7 Composition (total molalities and pH) of the solution entering the domain (incoming) and the initial fluid in the pores of the marla

pH Al Si Na K Ca Mg Carbonate Sulfate Cl a

Incoming

Initial

13.44 1.14103 3.91104 5.87101 3.75101 1.01103 1.30108 2.53104 1.16103 4.19101

7.50 6.50109 1.06104 4.61101 7.10104 5.05103 2.65104 3.39103 1.07105 4.69101

Chlorine is used for charge-balance and is non-reactive.

Fig. 2. Mineral content vs. distance along the fracture at t=5000 a. Case 1.

1563

alteration extends up to about 10 cm laterally from the fracture near the fracture inlet after 5000 a of reaction) than in the low pH cases. 4.1. Case 1: pH 13.44, NaCl-type groundwater The composition of the solution entering the system and the initial composition of the solution in the pores of the rock is shown in Table 7. The incoming solution is at equilibrium with respect to calcite and undersaturated with respect to the other primary minerals in the marl. The initial solution is at equilibrium with the marl. Fig. 2 shows mineral content along the fracture after 5000 a of circulation of the high-pH solution. Dolomite has been completely dissolved and partially replaced by calcite in the first 4 m along the fracture. There has also been minor muscovite and quartz dissolution. Regarding the Mg-containing secondary phases (Fig. 2b), there is a brucite-rich zone close to the fracture inlet, followed by a sepiolite-rich zone. Natrolite and tobermorite are also present in the first half of the domain and there is an analcime-containing zone in the last 1.5 m of the domain. The changes in mineral content have caused a slight increase of porosity in the first 2 m along the fracture and a slight decrease afterwards (Fig. 2a). These changes in porosity and the associated changes in permeability have translated into a significant reduction of the Darcy velocity. Water still flows preferentially along the fracture (Fig. 3), but the Darcy velocity is now a factor of 7 smaller than the initial velocity.

Fig. 3. Velocity field for case 1, at t=5000 a. Fluid flows preferentially through the fracture. Flow outside the fracture is negligible (length of the arrows is proportional to the magnitude of the Darcy velocity). The dimensions of the domain are 6 m per 1 m (labels on the axes indicate only fraction of total length).

1564

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

Fig. 4 shows mineral reaction rates along the fracture. Negative rates mean dissolution and positive rates mean precipitation. The most visible feature in Fig. 4a is the region of dolomite dissolution and simultaneous precipitation of calcite, which shows a double peak. The rates for both reactions are practically identical. The peak in dolomite dissolution and calcite precipitation at x=4 m coincides also with a peak for brucite precipitation (Fig. 4b), although there is also minor precipitation of sepiolite and natrolite. The overall reaction in this

region of the fracture corresponding to these rate maxima can be approximately written as CaMgðCO3 Þ2 þ2OH ! CaCO3 þ MgðOHÞ2 þCO2 3 dolomite

calcite

brucite

In the second peak (x=4.5 m) there is no precipitation of brucite, but significant precipitation of natrolite instead. Silica is provided by the dissolution of analcime (replacement of analcime by natrolite).

Fig. 4. Reaction rates along the fracture for (a) primary and (b) secondary minerals. Case 1. Negative rates mean dissolution; positive rates mean precipitation. The smaller figures are the same plots with larger vertical scales, so the magnitude of the large dissolution and precipitation peaks can be observed.

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

The advance of the calcite-dolomite reaction front is shown in Fig. 5. The rate of advance of the front diminishes with time, as the flow velocity diminishes. The advance of the double reaction front can also be observed in the evolution of pH with time (Fig. 6). The location of the reaction fronts at each time is marked by inflections in the pH curves. Notice that the first reaction front advances more quickly than the second (the 2 inflections are closer together at later times). In Fig. 4b it can also be observed that there is a first small peak of brucite precipitation at x=2 m, corresponding to the dissolution of preexisting sepiolite. Between the two zones of brucite precipitation there is

1565

replacement (simultaneous dissolution/precipitation) of brucite by sepiolite. Regarding the other primary minerals, quartz and muscovite also dissolve along the profile. Regarding other secondary minerals, there is precipitation of natrolite and tobermorite in the first half of the domain. The precipitation of tobermorite, which takes Ca from solution, causes minor dissolution of calcite close to the fracture inlet. Total concentrations in solution for the different species are plotted in Fig. 7. Carbonate concentrations are very high and Ca concentrations are very low, due to the fact that only half the carbonate released into solution by the dissolution of dolomite is consumed by the precipitation of calcite (see reaction above), combined with the fast rate for the overall reaction (i.e., fast release of carbonate into solution). If an additional

Fig. 5. Location of the calcite–dolomite reaction front vs. time for cases 1 (input pH 13.44) and 2 (input pH 12.50).

Fig. 6. pH vs. distance along the fracture at times equal to 100, 1000 and 5000 a. Case 1.

Fig. 7. Concentrations (molalities) vs. distance along the fracture at t=5000 a. Case 1. Sulfate and Cl concentrations (not shown for clarity) are the same as in the incoming solution (see Table 7) and constant along the profile. Magnesium concentration is of the order of 107 m.

1566

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

Table 8 Composition (total molalities and pH) of the solution entering the domain (incoming) and the initial fluid in the pores of the marla

pH Al Si Na K Ca Mg Carbonate Sulfate Cl a

Incoming

Initial

12.50 1.32104 1.40105 4.61101 7.10104 3.51102 1.40107 9.64106 1.30105 4.70101

7.50 6.50109 1.06104 4.61101 7.10104 5.05103 2.65104 3.39103 1.07105 4.69101

Chlorine is used for charge-balance and is non-reactive.

Fig. 8. Mineral content vs. distance along the fracture at t=20 000 a. Case 2. The horizontal scale for the second plot (secondary minerals) goes down only to x=1 m (the extent of the zone of mineral alteration).

carbonate-bearing phase were to control carbonate concentrations in solution, the tendency would be to increase the dissolution of dolomite, increasing the release of Mg into solution and the precipitation of secondary Mg-bearing phases. 4.2. Case 2: pH 12.50, NaCl-type groundwater Table 8 shows the composition of the solution entering the system and the initial composition of the solution in the pores of the rock. The incoming solution is at equilibrium with respect to calcite and undersaturated with respect to the other primary minerals in the marl. The initial solution is at equilibrium with the marl. Fig. 8 shows mineral content along the fracture after 20 000 a of circulation of the high-pH solution. Notice that the spatial extent of the zone of mineral alteration (about 0.8 m along the fracture) is much smaller than in case 1 (input pH about 13.5), due to the pronounced reduction of porosity at the fracture inlet, which causes an extreme reduction of fluid flow through the system. The reduction of porosity at the fracture inlet causes also a channeling of the flow from the wall rock into the fracture (Fig. 9). The Darcy velocity along the fracture is about a factor of 75 smaller than the initial velocity. Fig. 8a shows that the dolomite in the fracture has been completely dissolved and replaced by calcite in the first half meter of the domain. There has also been minor reaction of quartz (dissolution) and muscovite (dissolution and precipitation). Fig. 8b shows the presence of secondary tobermorite, natrolite, sepiolite and minor brucite and analcime.

Fig. 9. Velocity field for case 2, at t=20 000 a. Fluid channels into the fracture at the inlet side of the domain. The length of the arrows is proportional to the magnitude of the Darcy velocity. The dimensions of the domain are 6 m per 1 m (labels on the axes indicate only fraction of total length).

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

Reaction rates along the fracture are shown in Fig. 10. In Fig. 10a it can be observed that there is a zone of quartz dissolution next to the fracture inlet, followed immediately by a zone where dolomite dissolves and calcite precipitates. The magnitude of the calcite precipitation rate is twice the dolomite dissolution rate. The overall reaction can be approximately written as CaMgðCO3 Þ2 þCaþ2 ! 2CaCO3 þ Mgþ2 dolomite

calcite

Brucite had precipitated in earlier stages of the process. However, this brucite is later redissolved and the Mg is taken up by the precipitation of sepiolite (Fig. 8b). The zone of brucite dissolution also coincides with the dissolution of previously precipitated tobermorite, and it is followed by a zone where analcime replaces previouslyprecipitated natrolite. The sealing of the fracture and the reduction of fluid flow through the system is also reflected by the decrease

1567

in the rate of advance of the calcite–dolomite reaction front (Fig. 5), the evolution of pH (Fig. 11) and the composition of the solution along the fracture (Fig. 12). Total concentrations in solution for the different species after 20 000 a of reaction are plotted in Fig. 12. Carbonate concentrations in solution are much smaller than in the previous case, since all the carbonate released into solution by the dissolution of dolomite is consumed by the precipitation of calcite (see reaction above). 4.3. Additional calculations with smaller surface areas Two sets of additional calculations have been included in the study in order to have some insight into the effect of the initial surface areas for primary minerals on the model results. These additional calculations correspond to the 2 cases described above, but with initial surface areas for primary minerals 4.6 times smaller than the ones previously used. 4.3.1. Case 1 with smaller surface areas All the parameters for this simulation (except surface areas for primary minerals) are the same as for case 1 (see Section 4.1). Fig. 13 is a plot of mineral content along the fracture after 5000 a of circulation of the high-pH solution. Comparing with Fig. 2, it is easily observed that dolomite has been completely dissolved along the fracture in the domain and partially replaced by calcite. Quartz and muscovite have only experienced minor dissolution and the only secondary phases present are brucite and

Fig. 10. Reaction rates along the fracture for (a) primary and (b) secondary minerals. Case 2. Negative rates mean dissolution; positive rates mean precipitation. The horizontal scale goes down only to x=1 m (the extent of the zone of mineral alteration).

Fig. 11. pH vs. distance along the fracture at times equal to 1000, 5000, 10 000 and 20 000 a. Case 2. The horizontal scale goes down only to x=1 m (the extent of the zone of mineral alteration).

1568

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

Fig. 12. Concentrations (molalities) vs. distance along the fracture at t=20 000 a. Case 2. Sulfate and Cl concentrations (not shown for clarity) are the same as in the incoming solution (see Table 8) and constant along the profile. Sodium concentration is 0.44 m and also constant along the profile. Fig. 13. Mineral content vs. distance along the fracture at t=5000 a. Case 1 with smaller surface areas.

natrolite. The changes in porosity are only very small, which causes the flow velocity to be faster than in case 1. Darcy velocity along the fracture is about 1.17102 m/ a (it was about 1.7103 m/a in case 1 after 5000 a; initial velocity was 1.26102 m/a). Because of the fast flow velocity, the calcite–dolomite reaction front has advanced faster than in case 1 (it has actually moved out of the domain after about 1500 a). The main consequence of the smaller surface areas for the primary minerals is a smaller change in porosity, which means that flow velocity along the fracture is faster and less favorable for the performance of a repository than the case with larger surface areas. However, Darcy velocity has still slightly decreased compared to the initial velocity.

4.3.2. Case 2 with smaller surface areas All the parameters for this simulation (except surface areas for primary minerals) are the same as for case 2 (see Section 4.2). Fig. 14 shows mineral content along the fracture after 20 000 a of circulation of the high-pH solution. Comparing with case 2 (Fig. 8), it can be observed that both profiles are rather similar, with a small difference in the porosity near the fracture inlet (smaller amount of secondary minerals). Porosity at the fracture inlet is slightly larger for the case with smaller surface areas, which causes a slightly larger Darcy velocity along the fracture (5.3104 m/a, compared to 1.7104 m/a in case 2). As in the previous case, the consequence of the smaller surface areas for the primary minerals is a smaller

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

Fig. 14. Mineral content vs. distance along the fracture at t=20000 a. Case 2 with smaller surface areas. The horizontal scale for the second plot (secondary minerals) goes down only to x=1.5 m (the extent of the zone of mineral alteration).

change in porosity, which means that flow velocity along the fracture is faster and less favorable for the performance of a repository than the case with larger surface areas. The difference is, however, smaller than in the previous case.

5. Conclusions In the reference frame of the proposed low- and intermediate-level radioactive waste repository at Wellenberg, the alteration of a fractured marl due to the circulation of hyperalkaline solutions derived from the

1569

degradation of cement has been studied by means of 2-dimensional reactive transport simulations. A modified version of the GIMRT software package (Steefel and Yabusaki, 1996) has been used for the reactive transport calculations. Both diffusive/dispersive and advective solute transport are taken into account. Mineral reactions are described by kinetic rate laws. The changes in porosity caused by mineral dissolution and precipitation are used to update the fluid flow field during the calculations. Flow fields are calculated using a separate finite difference code. Four different cases have been considered. These 4 cases are representative of the 2 different types of groundwater (NaHCO3-type and NaCl-type) present at Wellenberg and 2 different stages in the process of degradation of cement (pH 12.5 and pH 13.5). Only the results for the NaCl-type groundwater have been shown. The results of the calculations show only small differences between the NaHCO3 system and the NaCl system. However, the results with the 2 different pH values are quite different. In both cases, the flow velocity in the fracture diminishes with time, due to a decrease in porosity. However, the decrease in porosity and flow velocity is much more pronounced in the lower pH case. Also, the extent of the zone of mineral alteration along the fracture is much more limited in the lower pH case (1–2 m along the fracture after 20 000 a). For the higher pH case, mineral alteration affects the whole length of the flow domain along the fracture (6 m) after only 5000 a. One of the main mineral reactions in these systems is the dissolution of the dolomite initially present in the rock (fracture) and its replacement (partial or total) by calcite. A markedly advancing reaction front is defined by this reaction in the higher pH system. The Mg released by the dissolution of dolomite is taken up by the precipitation of brucite and sepiolite (a Mg-phyllosilicate). Brucite dominates near the fracture inlet and at earlier stages of the process. Sepiolite is the dominant Mg-containing phase further from the fracture inlet and also at later stages (brucite tends to dissolve and be replaced by sepiolite with time in the lower pH cases). Other secondary minerals that precipitate during the simulations are analcime and natrolite (zeolites) and tobermorite (a CSH phase). The amount of reaction affecting quartz (dissolution) and muscovite (dissolution and precipitation) is relatively minor. Since the biggest source of uncertainty for the reaction rates arises from the value of the surface areas used for the primary minerals, additional calculations making use of smaller values of these areas have also been performed. The results show that the smaller surface areas result in an overall smaller reactivity of the system and smaller porosity changes. These smaller porosity changes mean than porosities near the fracture inlet are

1570

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571

larger than in the calculations with larger surface areas. The larger porosities translate into larger Darcy velocities along the fracture (less favorable for the performance of a repository), showing that the coupling between chemical reaction, which drives the changes in porosity, and fluid flow is very important. In summary, the results of the simulations suggest that the net result of the interaction of the high-pH solutions with the fractured marl is a reduction of porosity and water velocity along the fractures. This result would actually be highly beneficial for the performance of a repository hosted by the Valanginian marls at Wellenberg. The magnitude and scale (spatial and temporal) of this evolution is open to a large degree of uncertainty. However, the results of this modeling study are consistent in general terms (precipitation and nature of secondary phases, decrease in porosity/permeability) with the results from experimental work performed on comparable rocks (Adler, 2001) and with the results from other modeling studies (Steefel and Lichtner, 1994; Lichtner and Eikenberg, 1995; Soler, 1998, 2000; Soler and Ma¨der, 2002; Savage et al., 2002).

Acknowledgements The author wishes to thank the insightful comments on the manuscript by Andreas Jakob and Urs Berner and the detailed reviews by James G. Brown and an anonymous reviewer. Partial financial support by the Swiss National Cooperative for the Disposal of Radioactive Waste (Nagra) is also gratefully acknowledged. References Adler, M., 2001. Interaction of Claystone and Hyperalkaline Solutions at 30
C: A Combined Experimental and Modeling Study. PhD Diss., Univ. of Bern. Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York. Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York. Berner, U., 1987. Modelling porewater chemistry in hydrated cement. Mat. Res. Soc. Symp. Proc. 84, 319. Berner, U., 1990. A Thermodynamic Description of the Evolution of Pore Water Chemistry and Uranium Speciation During the Degradation of Cement. Paul Scherrer Institut Bericht No. 62 (also published as Nagra Technical Report NTB 90-12). Bloom, P.R., Erich, M.S., 1987. Effect of solution composition on the rate and mechanism of gibbsite dissolution in acid solutions. Soil Sci. Soc. Am. J. 51, 1131–1136. Carroll, S.A., Walther, J.V., 1990. Kaolinite dissolution at 25, 60 and 80
C. Am. J. Sci. 290, 797–810. Carroll-Webb, S.A., Walther, J.V., 1988. A surface complex reaction model for the pH-dependence of corundum and kaolinite dissolution rates. Geochim. Cosmochim. Acta 52, 2609–2623.

Chou, L., Garrels, R.M., Wollast, R., 1989. Comparative study of the dissolution kinetics and mechanisms of carbonates in aqueous solutions. Chem. Geol. 78, 269–282. Domenico, P.A., Schwartz, F.W., 1990. Physical and Chemical Hydrogeology. John Wiley & Sons, New York. Dove, P.M., 1994. The dissolution kinetics of quartz in sodium chloride solutions at 25 to 300
C. Am. J. Sci. 294, 665–712. Dullien, F.A.L., 1979. Porous Media. Academic Press, New York. Ganor, J., Mogollo´n, J.L., Lasaga, A.C., 1995. The effect of pH on kaolinite dissolution rates and on activation energy. Geochim. Cosmochim. Acta 59, 1037–1052. Horseman, S.T., Higgo, J.J.W, Alexander, J., Harrington, J.F., 1996. Water, Gas and Solute Movement through Argillaceous Media. Report CC-96/1. NEA-OECD, Paris. Kalinowski, B.E., Schweda, P., 1996. Kinetics of muscovite, phlogopite, and biotite dissolution and alteration at pH 1-4, room temperature. Geochim. Cosmochim. Acta 60, 367–385. Kirkner, D.J., Reeves, H., 1988. Multicomponent mass transport with homogeneous and heterogeneous chemical reactions: effect of the chemistry on the choice of numerical algorithm. 1. Theory. Water Resour. Res. 24, 1719–1729. Knauss, K.G., Wolery, T.J., 1989. Muscovite dissolution kinetics as a function of pH and time at 70
C. Geochim. Cosmochim. Acta 53, 1493–1502. Lasaga, A.C., 1998. Kinetic Theory in the Earth Sciences. Princeton University Press, Princeton. Lichtner, P.C., 1985. Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems. Geochim. Cosmochim. Acta 49, 779–800. Lichtner P.C., Eikenberg J., 1995. Propagation of a Hyperalkaline Plume into the Geological Barrier Surounding a Radioactive Waste Repository. Paul Scherrer Institut Bericht No. 95-1. Wu¨relingen. Lin, F.C., Clemency, C.V., 1981. The kinetics of dissolution of muscovites at 25
C and 1 atm CO2 partial pressure. Geochim. Cosmochim. Acta 45, 571–576. de Marsily, G., 1986. Quantitative Hydrogeology. Academic Press, San Diego. Mogollo´n, J.L., Ganor, J., Soler, J.M., Lasaga, A.C., 1996. Column experiments and the full dissolution rate of gibbsite. Am. J. Sci. 296, 729–765. Morse, J.W., Mackenzie, F.T., 1990. Geochemistry of Sedimentary Carbonates. Developments in Geology 48. Elsevier, Amsterdam. Nagra, 1993. Untersuchungen zur Standorteignung im Hinblick auf die Endlagerung Schwach- und Mittelaktiver Abfa¨lle—Geologische Grundlagen und Datensatz zur Beurteilung der Langzeitsicherheit des Endlagers fu¨r Schwach und Mittelaktive Abfa¨lle am Standort Wellenberg (Gemeinde Wolfenschiessen, NW). Nagra Technischer Bericht NTB 93-28. Nagra, Wettingen. Nagra, 1994. Endlager fu¨r Schwach- und Mittelaktive Abfa¨lle (Endlager SMA)- Bericht zur Langzeitsicherheit des Endlagers SMA am Standort Wellenberg (Gemeinde Wolfenschiessen, NW). Nagra Technischer Bericht NTB 94-06. Nagra, Wettingen. Nagra, 1997. Geosynthese Wellenberg 1996. Ergebnisse der Untersuchungsphasen I und II. Nagra Technischer Bericht NTB 96-01. Nagra, Wettingen. Nagy, K.L., Blum, A.E., Lasaga, A.C., 1991. Dissolution and

J.M. Soler / Applied Geochemistry 18 (2003) 1555–1571 precipitation kinetics of kaolinite at 80
C and pH 3: the dependence on solution saturation state. Am. J. Sci. 291, 649–686. Nagy, K.L., Lasaga, A.C., 1992. Dissolution and precipitation kinetics of gibbsite at 80
C and pH 3: the dependence on solution saturation state. Geochim. Cosmochim. Acta 56, 3093–3111. Neall, F.B., 1994. Modelling of the Near-Field Chemistry of the SMA Repository at the Wellenberg Site. Paul Scherrer Institut Bericht 94-18 (also published as Nagra Tech. Rep. NTB 94-03). Nickel, E., 1973. Experimental dissolution of light and heavy minerals in comparison with weathering and intrastratal solution. Contrib. Sedimentol 1, 1–68. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publications, New York. Reed, M.H., 1982. Calculation of multicomponent chemical equilibria and reaction process in systems involving minerals, gases and aqueous phases. Geochim. Cosmochim. Acta 46, 513–528. Rimstidt, J.D., Barnes, H.L., 1980. The kinetics of silica-water reactions. Geochim. Cosmochim. Acta 44, 1683–1699. Savage, D., Noy, D., Mihara, M., 2002. Modelling the interaction of bentonite with hyperalkaline fluids. Appl. Geochem. 17, 207–223. Sjo¨berg, E.L., 1978. Kinetics and mechanism of calcite dissolution in aqueous solutions at low temperatures. Stockholm Contrib. Geol 32, 1–96. Soler, J.M., 1998. Reactive transport modelling of the interaction between a high pH plume and a fractured marl. Mineral. Mag. 62A, 1427–1428. Soler, J.M., 2000. One-dimensional reactive transport modelling of the interaction between a high-pH plume and a fractured granodiorite. The GTS-HPF project. Goldschmidt 2000, September 3rd–8th, 2000, Oxford, UK. J. Conf. Abstracts 5, 947.

1571

Soler, J.M., Ma¨der, U.K., 2002. The GTS-HPF experiment: reaction-induced permeability changes. Goldschmidt 2002, August 18th–23rd, 2002, Davos, Switzerland. Geochim. Cosmochim. Acta 66, A726. Steefel, C.I., 1992. Coupled Fluid Flow and Chemical Reaction: Model Development and Application to Water-Rock Interaction. PhD Diss., Yale Univ., New Haven. Steefel, C.I., Lasaga, A.C., 1994. A coupled model for transport of multiple chemical species and kinetic precipitation/ dissolution reactions with applications to reactive flow in single phase hydrothermal systems. Am. J. Sci. 294, 529– 592. Steefel, C.I., Lichtner, P.C., 1994. Diffusion and reaction in rock matrix bordering a hyperalkaline fluid-filled fracture. Geochim. Cosmochim. Acta 58, 3595–3612. Steefel, C.I., Yabusaki, S.B., 1996. OS3D/GIMRT, Software for Multicomponent–Multidimensional Reactive Transport: User’s Manual and Programmer’s Guide. PNL-11166. Pacific Northwest National Laboratory, Richland. Stumm, W., Wehrli, B., Wieland, E., 1987. Surface complexation and its impact on geochemical kinetics. Croatica Chem. Acta 60, 429–456. Ullman, W.J., Aller, R.C., 1982. Diffusion coefficients in nearshore marine sediments. Limnol. Oceanog. 27, 552–556. White, A.F., Brantley, S.L. (Eds.), 1995. Chemical Weathering Rates of Silicate Minerals. Reviews in Mineralogy, Vol. 31. Mineralogical Society of America, Washington, DC. Wieland, E., Stumm, W., 1992. Dissolution kinetics of kaolinite in acidic aqueous solutions at 25
C. Geochim. Cosmochim. Acta 56, 3339–3355. Wolery, T.J., Jackson, K.J., Bourcier, W.L., Bruton, C.J., Viani, B.E., Knauss, K.G., Delany, J.M., 1990. Current status of the EQ3/6 software package for geochemical modeling. In: Melchior, C., Bassett, R.L. (Eds.), Chemical Modeling of Aqueous Systems II. ACS Symposium Series No. 416. ACS, pp. 104–116.