2D reactive transport modeling of the interaction between a marl and a CO2-rich sulfate solution under supercritical CO2 conditions

2D reactive transport modeling; Dávila et al. (2016) 1 2D reactive transport modeling of the interaction between a marl and a 1 CO2-rich sulfate solution under supercritical CO2 conditions 2 Gabriela Dávila a,b,c*, Linda Luquot a,b, Josep M. Soler a,b and Jordi Cama a,b 3 a Institute of Environmental Assessment and Water Research (IDAEA), CSIC. 4 Jordi Girona 18-26, 08034 Barcelona, Catalonia, Spain. 5 b Associated Unit: Hydrogeology Group (UPC-CSIC). 6 c Department of Civil and Environmental Engineering, Universitat Politècnica de Catalunya (UPC). 7 Jordi Girona 1-3, 08034 Barcelona, Catalonia, Spain. 8 *Corresponding author: Gabriela Dávila ([email protected]) 9 Jordi Girona 18-26, Barcelona 08034. Catalonia, Spain 10 Tel #:+34 934006100 11 Fax #: +34932045904 12 Abstract 13 The circulation of CO2-rich solutions through fractured marl cores (caprock) under 14 different flow rates and supercritical CO2 conditions (PTotal = 150 bar, pCO2 = 61 bar 15 and T = 60 °C) led to mineral changes caused mainly by calcite dissolution and to a 16 lesser extent by aluminosilicate dissolution, and by gypsum precipitation adjacent to the 17 fracture walls. Another significant result was the formation of the altered and highly 18 porous zone (Dávila et al., 2016a). Dissolution structures ranged from face to uniform 19 dissolution and wormhole formation depending mainly on the flow rate. 20 2D reactive transport models were used to interpret the results of the percolation 21 experiments (except at 60 mL h-1). They reproduced the variation in the outflow 22 composition with time and the observed width of the altered zone along the fractures. A 23 good match was achieved by using initial Deff values in the rock matrix that ranged from 24 1 × 10-13 m2 s-1 to 3 × 10-13 m2 s-1 under slow flow rates. The Deff value was higher by a 25 factor of 20 (6 × 10-12 m2 s-1) under fast flow. Moreover, a slight variation in the calcite 26 reactive surface areas contributed to the fit of the model to the experimental data. 27 The modeling results reproduced major dissolution of calcite and gypsum precipitation, 28 and minor dissolution of clinochlore. Calcite dissolution was boosted by increasing the 29 flow rate and gypsum precipitation increased at intermediate flow rate (1 mL h-1). Minor 30    2D reactive transport modeling; Dávila et al. (2016) 2 precipitation of dolomite, kaolinite, boehmite and three zeolites (mesolite, stilbite and 31 smectite) along the altered zone occurred. The magnitude of these reactions was 32 consistent with the measured increase in porosity over the altered zone. 33 Keywords: CO2 sequestration, numerical modeling, leakage, marl caprock, calcite 34 dissolution, gypsum precipitation. 35 1.Introduction 36 Leakage of injected CO2 in deep reservoirs may occur through preferential pathways 37 such as faults and fractures (Rutqvist and Tsang, 2002) and through the rock cement 38 interface (Dávila et al., 2016b). Dávila et al. (2016a) showed that the permeability of the 39 Hontomín marl wt.% calcite, < 25 wt.% clay and 10 wt.% 40 quartz was at least 6 orders of magnitude smaller than fracture permeability (kf) in 41 fractured cores. Fluids will therefore flow preferentially through fractures, where the 42 interaction between CO2-rich fluids and the rock matrix can bring about changes in 43 physical and chemical properties caused by dissolution and precipitation (Singurindy 44 and Berkowitz, 2005; Edlmann et al., 2013; Kampman et al., 2014; Chen at al., 2014; 45 Dávila et al., 2016a). The interpretation of the processes is not straightforward given 46 that caprocks are usually composed of a large number of minerals that react differently. 47 Reactive transport modeling has been performed to assess the impact of the 48 interaction between CO2-rich brines and caprocks at the laboratory (Gherardi et al., 49 2007; Credoz et al., 2009; Hao et al., 2013; Tian et al., 2014) and field scales (Gaus et 50 al., 2005). 51 Gherardi et al. (2007) modeled fluid-rock interactions in intact and fractured 52 carbonate-rich shale caprocks using 1D and 2D models, respectively. The dominant role 53 of calcite outweighed the effects of the other mineralogical changes involving 54 precipitation/dissolution of Al-silicate minerals. Minor mineral transformations induced 55 by the advance of the CO2-induced acidic front in the caprock included clay dissolution 56 (illite, chlorite and muscovite) and precipitation reactions (Na-smectite). Other effects 57 entailed the formation of new minerals such as dawsonite, siderite, and ankerite. 58 Credoz et al. (2009) performed and simulated batch experiments with clayey and 59 clay limestone crushed caprocks at pCO2 = 150 bar and T = 80 °C for 90 days. The 60 reactivity of the two rock samples was similar under these experimental conditions. 61 2D reactive transport modeling; Dávila et al. (2016) 3 Dissolution of calcite and dolomite predominated together with minor or negligible 62 precipitation of kaolinite and montmorillonite due to the destabilization of Ca-63 montmorillonite and illite. At a large scale, Tian et al. (2014) performed 1D numerical 64 simulations of a 30 m long vertical column with 1 m2 of cross-sectional area to study 65 the effect of the mineral composition of clay-rich shale and mudstone caprocks at 101 66 bar and 47 °C. The mineral composition of the caprock controlled the trend of porosity 67 change. The buffering capacity of the clay-rich shale was higher than that of the 68 mudstone. The porosity was influenced more by the mineral assemblage of the 69 mudstone than by the clay-rich shale. This study suggested that the mudstone was more 70 suitable for use as a caprock. 71 Hao et al. (2013) used a 3D continuum-scale reactive transport model to 72 simulate core flood experiments for marly dolostone and vuggy limestone reacting with 73 brines equilibrated with pCO2 = 5, 10, 20, 30 bar and T = 60 °C that were performed by 74 Smith et al. (2013). The authors showed the role of physical heterogeneities on the 75 carbonate rocks in the development of the dissolution fronts. Stable dissolution fronts 76 developed in the marly dolostone were fitted using an empirical exponent (n = 3) in the 77 Kozeny Carman equation for porous media (K = Ko·( / o)n, where K and are 78 permeability and porosity). For more impermeable and heterogeneous cores (vuggy 79 limestone) n = 6 8. 80 At field scale, Gaus et al. (2005) evaluated the long term geochemical behavior 81 of silty clay caprock at the Sleipner injection site by performing 1D reactive transport 82 modeling combining reaction kinetics and diffusive transport under CO2 supercritical 83 conditions (PTotal = 101.3 bar; pCO2 = 100 bar and T = 37 °C). If the caprock is assumed 84 to be a homogeneous medium, feldspar alteration would be the dominant long term 85 reaction despite initial carbonate dissolution. These reactions could cause a slight 86 decrease in porosity and a subsequent decrease in permeability. 87 In our study, 2D reactive transport modeling was performed to quantify the 88 dissolution and precipitation processes that occurred in the laboratory percolation 89 experiments in which CO2-rich solutions circulated through fractured marl cores (Dávila 90 et al., 2016a). The changes in the fractures and in the porosity of the rock matrix 91 induced by dissolution of calcite, clinochlore, albite and gypsum (in S-free solutions) 92 2D reactive transport modeling; Dávila et al. (2016) 4 and precipitation of gypsum (in S-rich gypsum-equilibrated solutions) and silicate 93 minerals led to variations in permeability. 94 2.Summary of the experimental results 95 Dávila et al. (2016a) showed that the circulation of CO2-rich solutions with different 96 sulfate contents through fractured cores under supercritical conditions (PTotal = 150 bar, 97 pCO2 = 61 bar and T = 60 °C) and different flow rates (0.2, 1 and 60 mL h-1) resulted in 98 a variety of dissolution patterns. In all cases, dissolution of calcite fostered the 99 formation of a highly porous reacted zone made up of slowly dissolving grains of 100 quartz, illite and pyrite and possible secondary minerals (e.g. gypsum). The dissolution 101 and precipitation reactions controlled the variations in the final pore volume associated 102 with the reacted cores. In all experiments, an increase in the flow rate led to an increase 103 in the calcite dissolution rate and in the final pore volume associated with the reacted 104 core. In the S-rich solution, the volume of dissolved calcite was always larger than the 105 volume of precipitated gypsum. The dissolution of calcite grains and, to a lesser extent, 106 the dissolution of clinochlore and albite, created a high porosity zone ( ranged from 48 107 to 67 %) along the fracture walls. This porosity zone was composed of non-dissolved 108 quartz, pyrite and illite grains, precipitated gypsum (in S-rich solutions) and silicate 109 phases (e.g., kaolinite or amorphous SiO2) and Al-bearing phases (e.g., boehmite). 110 As regards the dissolution patterns, in the S-free solution experiments face dissolution 111 occurred at slow flow rates (peclet number (Pe) of 4 and 21). Wormhole formation was 112 observed at the highest flow rate (Pe = 1267). However, in S-rich solution experiments, 113 face dissolution occurred only at the slowest flow rate. At an intermediate flow rate, 114 both uniform dissolution and wormhole formation took place. A combination of 115 wormhole formation and uniform dissolution occurred at the highest flow rate. Local 116 heterogeneities controlled mineral dissolution along the fracture, which could lead to 117 unexpected dissolution patterns. 118 As regards fracture permeability, under slow flow rates (low Pe) and in S-free solution 119 experiments, kf did not significantly change since face dissolution did not cause fracture 120 aperture to vary. Moreover, detached grains along the fracture could have led to the 121 obstruction of the fluid flow, preventing kf from increasing. In S-rich solutions, a 122 marked decrease in kf was attributed to gypsum precipitation. Under fast flow rate (high 123 2D reactive transport modeling; Dávila et al. (2016) 5 Pe), kf slightly increased with slight variations in aperture in the S-free solution 124 experiment. The increase in kf in the S-rich experiment was attributed to the absence of 125 gypsum precipitates. This increase in kf suggested the existence of flow through the 126 altered rock matrix under these fast flow conditions. 127 3.Reactive transport modeling 128 The simulations were performed using CrunchFlow code (Steefel et al., 2015). The 129 initial mineral composition of the rock was calculated from the initial mineral 130 composition obtained from the XRD-Rietveld analysis. 131 3.1 Description of the CrunchFlow reactive transport code 132 The reactive transport modeling was performed using the CrunchFlow code (Steefel et 133 al., 2015), which solves numerically the mass balance of solutes expressed as 134 j jj j RCC t C q D (j (1) 135 where is porosity, Cj is the concentration of component j (mol m-3), q is the Darcy 136 velocity (m3 m-2 s-1), Rj is the total reaction rate affecting component j (mol m-3 rock s-1) 137 and D is the combined dispersion-diffusion coefficient (m2 s-1). 138 The total reaction rate Rj is given by 139 m m jmj RR (2) 140 where Rm is the rate of precipitation (Rm > 0) or dissolution (Rm < 0) of mineral m in mol 141 m-3 rock s-1 , and jm is the number of the moles of j in mineral m. 142 The reaction rate laws used in the calculations are expressed as 143 1 2 ,1 m m eq i n i terms n HTmmm K IAP aakAR iH (3) 144 where Am is the mineral surface area in m2mineral m-3bulk, km,T is the reaction rate constant 145 in mol m-2 s-1, H n H a is the term describing the effect of pH on the rate, i n i a is the term 146 describing a catalytic/inhibitory effect by another species on the rate, IAP is the ionic 147 2D reactive transport modeling; Dávila et al. (2016) 6 activity product of the solution with respect to the mineral, Keq is the equilibrium 148 constant for the dissolution reaction (ionic activity product at equilibrium) and m2 and 149 m1 are the parameters affecting the dependence of the rate on solution saturation state. 150 The rate constant at temperature T (K) is calculated from 151 TTR Ea kk mTm 11 exp 25 25,, (4) 152 where km,25 is the rate constant at 25 oC, Ea is the apparent activation energy of the 153 overall reaction (J mol-1) and R is the gas constant (J mol-1 K-1). 154 Change in mineral surface area (Am in m2mineral m-3bulk) owing to dissolution is 155 3 2 )( 3 2 )( imi m initial mAA (5) 156 Change due to precipitation is 157 3 2 )(i initial mAA (6) 158 where (i)m is the initial volume fraction of the mineral m and (i) is the initial porosity of 159 the medium. This formulation ensures that as the volume fraction of a mineral reaches 160 zero, so does its surface area. Moreover, in the case of both dissolving and precipitating 161 minerals, the term ( / (i))2/3 requires the surface area of a mineral in contact with fluid 162 to reach zero when the porosity of the medium reaches zero. This formulation is used 163 primarily for primary minerals (i.e., minerals with initial volume fractions > 0). In the 164 case of secondary minerals which precipitate, the value of the initial bulk surface area 165 specified is used provided that precipitation occurs. If this phase subsequently dissolves, 166 the above formulation is used with an arbitrary initial volume fraction of 0.01. 167 3.2 Numerical discretization 168 The 3D cylindrical core samples used for the laboratory experiments (Fig. 1a) were 169 converted to a 2D rectangular symmetry. First, half of the circular section of the core 170 was transformed into a rectangle. The area of this rectangle was equal to the area of the 171 2D reactive transport modeling; Dávila et al. (2016) 7 half circle (area = x·l = ½· ·(l/2)2; where l is the diameter of the core and x is the width 172 of the rectangle). A 2D longitudinal section of the resulting prism was used as the 173 calculation domain (Fig. 1b and c). This domain was composed of a fracture (high 174 permeability zone) and rock matrix (low permeability zone; Fig. 1d). The domain was 175 divided into 36 and 30 nodes from the inlet to the outlet (L) and from the center of the 176 fracture to 3.5 mm in the model, respectively (Table 1). Fracture corresponded only to 177 the first node in the x direction. The model considered advection and dispersion only 178 along the fracture. Solute transport in the rock matrix was caused only by diffusion. 179 3.3 Parameters 180 3.3.1.Rock composition 181 The rock considered in the models was the marl described in Dávila et al. (2016a). It is 182 composed of calcite (71.2 wt.%), quartz (9.7 wt.%), illite (7.1 wt.%), albite (6.5 wt.%), 183 gypsum (2.0 wt.%), clinochlore (2.8 wt.%), anhydrite (0.5 wt.%), and pyrite (0.2 wt.%). 184 The initial porosity of the rock matrix ( (i)) was estimated by image segmentation 185 processing applied to ESEM images, which led to the separation of rock and voids. The 186 initial mineral volume fractions were calculated from the estimated porosity (6.7 %; 187 Table 2). An analysis of the selection of the secondary minerals was performed based 188 on equilibrium batch modeling (Gaus, 2010; Soler et al., 2011; Tian et al., 2014). Initial 189 estimates of the reactive surface areas of the primary minerals were calculated by 190 assuming spheres with radii estimated from the ESEM images (Table 2). The secondary 191 minerals were allowed to start precipitating when the solution reached supersaturation. 192 The initial reactive surface areas for all the secondary minerals were assumed to be the 193 same and sufficiently high (1.0 × 104 m2mineral m-3rock) to allow fast precipitation (local 194 equilibrium). Tutolo et al. (2015) performed simulations and observed that when 195 solutions were supersaturated with respect to Al-bearing minerals (boehmite, gibbsite 196 and diaspore), only boehmite was taken into account since its precipitation was 197 observed under high temperature and pressure conditions. 198 3.3.2.Solution composition 199 The two solutions used in the models were the ones used by Dávila et al., 2016a (Table 200 3). In both solutions, the initial concentration of CO2 was calculated to be 6.52 × 10-1 201 mol kgw-1 according to the Duan and Sun (2003) model with an imposed pCO2 of 61 202 2D reactive transport modeling; Dávila et al. (2016) 8 bar and bearing in mind that PTotal = 150 bar, T = 60 °C and I = 0.6 M. Also, the initial 203 concentration of O2 was estimated to be 3.06 × 10-4 mol L-1 (atmospheric and CO2-204 bottle concentrations). The S-rich solution was undersaturated with respect to calcite 205 (SICal = -3.20) and near equilibrium with respect to gypsum (SIGp = -0.02) at 60 °C. The 206 S-free injected solution was undersaturated with respect to both calcite (SICal = -3.50) 207 and gypsum (SIGp = -7.54) at 60 °C. The calculated initial pH values were 3.26 and 3.29 208 for the S-rich and S-free solutions, respectively. 209 Since the initial pore water composition of these samples was unknown and disturbed 210 during sample preparation, it was assumed to be in equilibrium with calcite and gypsum 211 for the S-rich solution experiments at room T and atmospheric pCO2 (concentration of 212 CO2 = 4.14 × 10-4 mol kgw-1), yielding a pH 7.7 (Table 3). A very small concentration 213 (10-6 mol kgw-1) was assumed for Na+, K+, Al3+, Cl-, Br-, Fe2+ and SiO2(aq). 214 3.3.3.Flow and transport properties 215 Darcy velocity (q), initial effective diffusion coefficient (Deff(i)), and longitudinal 216 dispersivity used in the simulations are shown in Table 2. The flow field used in the 217 reactive transport assumed constant Darcy velocity in the fracture and no flow in the 218 rock matrix. The effective diffusion coefficient (Deff) at 60 °C in the longitudinal and 219 transversal directions was calculated as o m eff DD , where m is the cementation 220 exponent (m = 2.5; Revil and Cathles, 1999) and Do is the diffusion coefficient in water. 221 3.3.4.Thermodynamic and kinetic data 222 One hundred and seven aqueous species were considered in the simulations. The 223 equilibrium constants were taken from the EQ3/6 thermodynamic database (Wolery et 224 al., 1990) included in CrunchFlow and are listed in Table A1 (Appendix 1). The activity 225 coefficients were calculated using the extended Debye-Hückel formulation (b-dot 226 model; Helgeson, 1969) with parameters from the same database. Twenty solid phases 227 (eight primary minerals (calcite, quartz, illite, albite, gypsum, clinochlore, anhydrite and 228 pyrite) and twelve secondary minerals (kaolinite, SiO2(am), dolomite, mesolite, stilbite, 229 smectite, mordenite, scolecite, analcime, wairakite, laumontite, gismondine and 230 boehmite) were considered in the calculations. Equilibrium constants for the mineral 231 dissolution reactions (CrunchFlow, EQ3/6) are given in Table A2 (Appendix 1). The 232 gypsum equilibrium constant (log KGp) value at 60 °C used in this study were lower 233 2D reactive transport modeling; Dávila et al. (2016) 9 than that initially implemented in the database and were log KGp(60°C) = 4.7383 as 234 reported by Nordstrom (2013) and Garcia et al. (2014). 235 Kinetic rate laws for the primary and secondary minerals, rate parameters and 236 activation energies are listed in Table 4. The mineral dissolution and precipitation rates 237 were taken from the literature for each mineral (Palandri and Kharaka, 2004; Xu et al., 238 2012; Bandstra et al., 2008; Bibi et al., 2011; Chou and Wollast, 1985; Hellmann et al., 239 2010; Hamer et al., 2003; Zhang et al., 2015; Domènech et al., 2002 ; Cama et al., 240 2000). The parallel rate laws for minerals describe the explicit dependence of the rates 241 on pH. Rate constants at a temperature different from 25 °C were calculated according 242 to Eq. (4). 243 4.Results and discussion 244 4.1 Dissolution and precipitation reactions 245 4.1.1S-rich injected solution 246 Adjustment of the values of the mineral reactive surface areas (Am) and the initial 247 effective diffusion coefficient (Deff) were used to match the variation in aqueous Ca, S, 248 Mg, K, Si and Fe concentrations with time (Table 2). 249 Fig. 2 shows the match between the experimental and modeled variations in the output 250 concentrations over time in the S-rich injected solution experiment at 1 mL h-1 (exp. 4). 251 The adjusted surface area of calcite was decreased by one order of magnitude from the 252 calculated geometric surface area value to match the output Ca concentration (Table 2), 253 which was higher than the input one owing to calcite dissolution (Fig. 2). The increase 254 in Ca concentration for the first 8 h was associated with the initially large reactive 255 surface area of calcite, resulting in significant calcite dissolution and in a decrease in 256 calcite content in the rock matrix close to the fracture. With time, the reaction was 257 controlled by diffusion through the rock matrix before slowing down. 258 The experimental and simulated output S concentrations were always smaller than the 259 input ones (Fig. 2). The deficit of S was attributed to gypsum precipitation given the 260 excessive Ca concentration and the positive values of the gypsum saturation index (0 < 261 SI the core during the experiment (Fig. S1; supplementary data). 262 2D reactive transport modeling; Dávila et al. (2016) 16 5.From the laboratory to the field scale 449 The experimental and modeling results of the percolation experiments shown in the 450 current paper and in the previous one by Dávila et al. (2016a) indicate that both the 451 formation of dissolution patterns and the variation in fracture permeability are highly 452 dependent on the flow rate which, being the fracture volume practically the same in all 453 the experiments, determines the residence time of the fluid circulating through the 454 fractures. An increase in the residence time (slow flow rate) provokes significant 455 dissolution of calcite only at the fracture inlet, forming face dissolution patterns, and 456 allows for gypsum precipitation, yielding either little variation or a reduction in fracture 457 permeability. Oppositely, shorter residence times (fast flow rate) do not allot for large 458 changes in solution composition, resulting in the formation of wormholes or uniform 459 dissolution and an increase in fracture permeability. 460 The residence time in the experiments ranged from 0.1 to 12 s (Table 2), which is 461 extremely short compared to likely residence times in caprock fractures during CO2 462 injection (e.g., from years, Bachu et al., 1994). Hence, according to our short-term 463 laboratory-scale results, it can be concluded that the slower flow rates in the repository 464 system will lead to limited calcite dissolution and effective gypsum precipitation, which 465 would in turn promote little change or even a decrease in fracture permeability. This 466 important role of flow velocity has already been pointed out by Brunet et al. (2016) for 467 Portland cement - CO2 - brine interactions. 468 At this point, however, we would like to emphasize that the 2D reactive transport 469 modeling applied to quantify the ongoing processes in the short fractured percolation 470 experiments under pCO2 and temperature similar to those at the field site cannot be 471 directly used to predict the long-term behavior in caprock fractures at the field scale. 472 For such modeling, parameters like (1) the dimensions (aperture and length) and 473 morphology of the field-scale fractures, (2) the reactive mineral surface areas, and (3) 474 the fluid flow regime should be known. Without this information, a sensitivity analysis 475 considering the possible variability of these parameters should be performed. Such a 476 modeling exercise is clearly beyond the objectives of this study. 477 2D reactive transport modeling; Dávila et al. (2016) 17 6.Summary and conclusions 478 Two dimensional reactive transport models were employed to interpret the results from 479 percolation experiments with single fracture marl cores during the injection of CO2-rich 480 solutions under supercritical conditions (PTotal = 150 bar, pCO2 = 61 bar and T = 60 °C). 481 In the model, flow circulated through the fracture, and solute transport in the rock 482 matrix was caused only by diffusion. Calculated solution concentration of the outlet 483 compares fairly well with the variation in the measured concentrations under slow flow 484 rates whereas under fast flow a poor match was obtained. The simulations reproduced 485 the dimensions of the dissolution patterns observed in ESEM and XMT images, except 486 for the S-free experiment run under high flow rate (60 mL h-1) in which a wormhole 487 formed (Dávila et al., 2016a) in a porous altered zone. These results considerably 488 improve our understanding of the different dissolution and precipitation processes 489 observed during the percolation experiments. 490 A successful match (except for the fast flow rate) between the experimental and 491 calculated output concentrations was achieved by using Deff values between 1 and 3 × 492 10-13 m2 s-1 and a single value for calcite irrespective of the type of solution with the 493 same reactive surface area (1.9 ± 1.5 × 104 m2mineral m-3 rock). In the fast flow rate 494 experiments (60 mL h-1) to provide more calcite reactivity, the Deff value showed a 495 slight increase (6 × 10-12 m2 s-1) along with the calcite reactive surface area (2.5 ± 0.1 × 496 105 m2mineral m-3 rock). The surface area values of the other minerals did not change. 497 Different compositions of the injected solutions produced diverse effects on the mineral 498 dissolution and precipitation reactions in the Hontomín caprock. The main reaction that 499 occurred was dissolution of calcite in both types of injected solution. For the same 500 experimental time and under the same flow rate, the volume of dissolved calcite was 501 always larger in the S-free solution experiments than in the S-rich ones. In S-rich 502 solution experiments, gypsum precipitated more at 1 mL h-1. Dissolution of clinochlore, 503 albite, illite and pyrite, and precipitation of dolomite, illite, kaolinite, mesolite, smectite, 504 stilbite and boehmite also occurred. Precipitation of secondary Al and Si rich minerals 505 took place to a lesser extent as was expected from the mass balance calculations, but 506 this was not observed in the ESEM and XMT images. 507 The calcite dissolution rate at the contact between the fracture and the rock matrix, 508 where the acid solution reacted with calcite under an advective flow regime, was faster 509 2D reactive transport modeling; Dávila et al. (2016) 18 than the rates calculated as the reaction moved towards the rock matrix. This occurred 510 because the accessibility of the solution to the calcite surface diminished during the 511 formation of the altered zone where solute transport was controlled by diffusion. This 512 phenomenon did not occur under the fast-flow rate experiment. 513 In both types of solution, the calculated porosity in the reacted zone was higher at the 514 fracture wall contact, decreasing with distance normal to fracture. An increase in the 515 flow rate caused a rise in porosity. The increase in porosity was higher at the inlet of the 516 core than at the outlet. This difference is diminished after speeding up the flow rate. The 517 rise in porosity from 6.7 70 % was similar to the volume of dissolved 518 60 %). 519 In contrast to the porosity increase, fracture permeability under slow flow rates tended 520 to decrease in the S-rich experiments and remained fairly constant in S-free 521 experiments. In both types of solution, the calculated porosity was greater under the fast 522 flow rate (60 mL h-1) than under slow flow rates and kf was considerably higher in the S-523 free injected solution. A new model accounting for the wormhole formation is 524 warranted to fully interpret the results from the experiments run under fast flow. 525 Acknowledgements 526 This work was financed by projects CGL2010-20984-C02-01 and CGL2014-54831-C3-527 1-R of the Spanish Government, project 2014SGR (Grup de Recerca SGR) 1456 of the 528 Catalan Government 529 Framework Programme FP7/2007-2013 under grant agreement number 282900). 530 Support was provided for GD by a JAEara la 531 LL by a Juan de la Cierva postdoctoral grant 532 (MINECO, Spain). We would like to thank Natàlia Moreno and Jordi Bellés (IDAEA), 533 and Eva Pelegrí and Maite Romero (SCT-Barcelona University) for analytical 534 assistance. We also wish to thank the anonymous reviewers for their constructive 535 comments that have improved the quality of the paper. 536 Appendix 1 537 Table A1 Equilibrium constants (log Keq) for the homogeneous reactions considered in 538 the reactive transport model. Reactions are written as the dissolution of 1 mol of the 539 2D reactive transport modeling; Dávila et al. (2016) 19 species in the table and in terms of Ca2+, Mg 2+, HCO3-, H+, SO42-, Na+, K+, Al3+, Cl-, 540 Br-, Fe2+, SiO2(aq) and O2(aq). 541 Table A2 Mineral and gas equilibrium constants (log Keq) considered in the reactive 542 transport model. Reactions are written as the dissolution of 1 mol of the species Ca2+, 543 Mg 2+, HCO3-, H+, SO42-, Na+, K+, Al3+, Cl-, Br-, Fe2+, SiO2(aq) and O2(aq). 544 Appendix 2 545 Figure B1 Simulated pH variation with respect to time in the S-rich (left) and in the S-546 free (right) injected solutions at the different flow rates (0.2, 1 and 60 mL h-1). 547 548 2D reactive transport modeling; Dávila et al. (2016) 20 References 549 Bachu, S., Gunter, W.D., Perkins, E.H., 1994. Aquifer disposal of CO2: Hydrodynamic 550 and mineral trapping. Energy Convers. Mgmt. 35, 269-279. 551 Bandstra, J.Z., Buss, H.L., Campen, R.K., Liermann, L.J., Moore, J., Hausrath, E.M., 552 Navarre-Sitchler, A.K., Jang, J-H., Brantley S.L., 2008. Appendix: Compilation of 553 Mineral Dissolution Rates. In S. L. Brantley, J. D. Kubicki, & A. F. White (Eds.), 554 Kinetics of Water Rock Interactions, 737 823. Springer. 555 Bibi, I., Singh, B., Silvester, E., 2011. Dissolution of illite in saline acidic solutions at 556 25°C. Geochim. Cosmochim. Acta 75, 3237 3249. 557 Brunet, J.P.L., Li, L., Karpyn, Z.K., Huerta, N.J., 2016. Fracture opening or self-558 sealing: Critical residence time as a unifying parameter for cement CO2brine 559 interactions. Int. J. Greenh Gas Con. 47, 25 37. 560 Cama, J., Ganor, J., Ayora, C., Lasaga, C.A., 2000. Smectite dissolution kinetics at 80 561 °C and pH 8.8. Geochim. Cosmochim. Acta 64, 2701 2717. 562 Chen, L., Kang, Q., Viswanathan, H.S., Tao, W.Q., 2014. Pore-scale study of 563 dissolution-induced changes in hydrologic properties of rocks with binary 564 minerals. Water Resour. Res. 50, 9343 9365. 565 Chou, L., Wollast, R., 1985. Steady-State kinetics and dissolution mechanisms of albite. 566 Am. J. Sci. 285, 963 993. 567 Credoz, A., Bildstein, O., Jullien, M., Raynal, J., Pétronin, J-C, Lillo, M., Pozo, C., 568 Geniaut, G., 2009. Experimental and modeling study of geochemical reactivity 569 between clayey caprocks and CO2 in geological storage conditions. Energy 570 Procedia 1, 3445 3452. 571 Dávila, G., Cama, J., Galíd, S., Luquot, L., Soler, J.M., 2016b. Efficiency of magnesium 572 hydroxide as engineering seal in the geological sequestration of CO2. Int. J. Greenh 573 Gas Con. 48, 171 185. 574 2D reactive transport modeling; Dávila et al. (2016) 21 Dávila, G., Luquot, L., Soler, J.M., Cama, J., 2016a. Interaction between a fractured 575 marl caprock and CO2-rich sulfate solution under supercritical CO2 conditions. Int. 576 J. Greenh Gas Con. 48, 105 119. 577 Domènech, C., Ayora, C., De Pablo, J., 2002. Oxidative dissolution of pyritic sludge 578 from the Aznalcóllar mine (SW Spain). Chem. Geol. 190, 339 353. 579 Duan, Z., Sun, R., 2003. An improved model calculating CO2 solubility in pure water 580 and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chem. 581 Geol. 193, 257 271. 582 Edlmann, K., Haszeldine, S., McDermott, C., 2013. Experimental investigation into the 583 sealing capability of naturally fractured shale caprocks to supercritical carbon 584 dioxide flow. Environ. Earth Sci. 70, 3393 3409. 585 Garcia-Rios, M., Cama, J., Luquot, L., Soler J.M., 2014. Interaction between CO2-rich 586 sulfate solutions and carbonate reservoir rocks from atmospheric to supercritical 587 CO2 conditions: experiments and modeling. Chem. Geol. 383, 107 122. 588 Gaus, I., 2010. Role and impact of CO2rock interactions during CO2 storage in 589 sedimentary rocks. Int. J. Greenh. Gas Control 4, 73 89. 590 Gaus, I., Azaroual, M., Czernichowski-Lauriol, I., 2005. Reactive transport modelling 591 of the impact of CO2 injection on the clayey caprock at Sleipner (North Sea). 592 Chem. Geol. 217, 319 337. 593 Gherardi, F., Xu, T.F., Pruess, K., 2007. Numerical modeling of self-limiting and self-594 enhancing caprock alteration induced by CO2 storage in a depleted gas reservoir. 595 Chem. Geol. 244, 103 129. 596 Hamer, M., Graham, R.C., Amrhein, C., Bozhilov, K.N., 2003. Dissolution of 597 Ripidolite (Mg, Fe-Chlorite) in Organic and Inorganic Acid Solutions. Soil Sci. 598 Soc. Am. J. 67, 654 661. 599 Hao, Y., Smith, M.M., Sholokhova, Y., Carroll, S.A., 2013. CO2-induced dissolution of 600 low permeability carbonates. Part II: Numerical modeling of experiments. Adv. 601 Water Resour. 62, 388 408. 602 2D reactive transport modeling; Dávila et al. (2016) 22 Helgeson, H.C., 1969. Thermodynamics model of hydrothermal systems at elevated 603 temperatures and pressures. Am. J. Sci. 267, 729 804. 604 Hellmann, R., Daval, D., Tisserand, D., 2010. The dependence of albite feldspar 605 dissolution kinetics on fluid saturation state at acid and basic pH: Progress towards 606 a universal relation. C. R. Geosci. 342, 676 684. 607 Kampman, N., Bickle, M., Wigley, M., Dubacq, B., 2014. Fluid flow and CO2fluid608 mineral interactions during CO2-storage in sedimentary basins. Chem. Geol. 369, 609 22 50. 610 Ketzer, J.M., Iglesias, R., Einloft, S., Dullius, J., Ligabue, R., De Lima, V., 2009. 611 Water rock CO2 interactions in saline aquifers aimed for carbon dioxide storage: 612 Experimental and numerical modeling studies of the Rio Bonito Formation 613 (Permian), southern Brazil. Appl. Geochem. 24, 760 767. 614 Luquot, L., Andreani, M., Gouze, P., Camps, P., 2012. CO2 percolation experiment 615 through chlorite/zeolite-rich sandstone (Pretty Hill Formation Otway Basin616 Australia). Chem. Geol. 294-295, 75 88. 617 Noiriel, C., Luquot, L., Madé, B., Raimbault, L., Gouze, P., van der Lee, J., 2009. 618 Changes in reactive surface area during limestone dissolution: An experimental 619 and modelling study. Chem. Geol. 265, 160 170. 620 Nordstrom, D.K., 2013. Improving internal consistency of standard state 621 thermodynamic data for sulfate ion, portlandite, gypsum, barite, celestine, and 622 associated ions. Procedia Earth Planet. Sci. 7, 624 627. 623 Palandri, J.L., Kharaka, Y.K., 2004. A compilation of rate parameters of water-mineral 624 interaction kinetics for application to geochemical modeling. science for a 625 changing world. 626 Revil, A., Cathles III, L.M., 1999. Permeability of shaly sands. Water Resour. Res. 35, 627 651 662. 628 Rutqvist, J., Tsang, C.F., 2002. A study of caprock hydromechanical changes associated 629 with CO2-injection into a brine formation. Environ. Geol. 42, 296-305. 630 2D reactive transport modeling; Dávila et al. (2016) 23 Singurindy, O., Berkowitz, B., 2005. The role of fractures on coupled dissolution and 631 precipitation patterns in carbonate rocks. Adv. Water Resour. 28, 507 521. 632 Smith, M.M, Sholokhova, Y., Hao, Y., Carroll, S.A., 2013. CO2-induced dissolution of 633 low permeability carbonates. Part I: Characterization and experiments. Adv. Water 634 Resour. 62, 370 387. 635 Soler, J.M., Vuorio, M., Hautojärvi, A., 2011. Reactive transport modeling of the 636 interaction between water and a cementitious grout in a fractured rock. Application 637 to ONKALO (Finland). Appl. Geochem. 26, 1115 1129. 638 Steefel, C.I., Appelo, C.A.J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., 639 Lagneau, V., Lichtner, P.C., Mayer, K.U., Meeussen, J.C.L., Molins, S., Moulton, 640 641 transport codes for subsurface environmental simulation. Computat. Geosci. 19, 642 445-478. 643 Tian, H., Xu, T., Wang, F., Patil, V.V., Sun, Y., Yue, G., 2014. A numerical study of 644 mineral alteration and self-sealing efficiency of a caprock for CO2 geological 645 storage. Acta Geotech. 9, 87 100. 646 Tutolo, B.M., Luhmann, A.J., Kong, X.Z., Saar, M.O., Seyfried, W.E., 2015. CO2 647 sequestration in feldspar-rich sandstone: Coupled evolution of fluid chemistry, 648 mineral reaction rates, and hydrogeochemical properties. Geochim. Cosmochim. 649 Acta 160, 132 154. 650 Wolery, T.J., Jackson, K.J., Bourcier, W.L., Bruton, C.J., Viani, B.E., Knauss, K.G., 651 Delany, J.M., 1990. Current status of the EQ3/6 software package for geochemical 652 modeling, in: Melchior, C., Bassett, R.L. (Eds.). Chemical Modeling of Aqueous 653 Systems II. ACS Symposium 416, 104 116. 654 Xu, J., Fan, C., Teng, H.H., 2012. Calcite dissolution kinetics in view of Gibbs free 655 energy, dislocation density, and pCO2. Chem. Geol. 322 323, 11 18. 656 Yu, Z., Liu, L., Yang, S., Li, S., Yang, Y., 2012. An experimental study of CO2 657 solution rock interaction at in situ pressure temperature reservoir conditions. 658 Chem. Geol. 326-327, 88 101. 659 2D reactive transport modeling; Dávila et al. (2016) 24 Zhang, S., Yang, L., De Paolo, D.J., Steefel, C.I., 2015. Chemical affinity and pH 660 effects on chlorite dissolution kinetics under geological CO2 sequestration related 661 conditions. Chem. Geol. 396, 208 217. 662 663 2D reactive transport modeling; Dávila et al. (2016) 25 Figure captions 664 Figure 1 Schemes showing: a) simplified experimental setup (more details in Dávila et 665 al., 2016a), b) fractured core sample and flow direction, and the spatial discretization 666 corresponding to the fractured marl cores: c) cylindrical coordinates, d) rectangular 667 coordinates and e) mesh distribution along the core sample. 668 Figure 2 Variation of the output concentrations with time under pCO2 of 61 bar and 60 669 °C in S-rich injected solution at 1 mL h-1 (exp. 4) for Ca, S, Mg, K, Fe and Si. Solid 670 symbols and solid line represent the experimental and calculated variations, 671 respectively. The dotted line represents the input solution concentration. 672 Figure 3 Simulated pH variations of the outlet solution with respect to time (left) and 673 with distance normal to fracture and at various positions along the fracture (right) Exp. 4 674 at Q = 1 mL h-1 and S-rich injected solution. 675 Figure 4 Variations in the output species concentrations with time under pCO2 of 61 676 bar and 60 °C during S-rich injected solution experiments at Q = 0.2 mL h-1 (exp. 2; 677 left) and at Q = 60 mL h-1 (exp. 7; right). Solid symbols and lines represent the 678 experimental and calculated (model A, B and C) variations, respectively. The dotted 679 line represents the input solution concentration. 680 Figure 5 Variations in the output concentrations with time under pCO2 of 61 bar and 60 681 °C in S-free injected solution experiments at 0.2, 1 and 60 mL h-1 (exps. 1, 3 and 6): Ca, 682 S, Fe and Si. 683 Figure 6 Variations of the simulated dissolution and precipitation rates of the primary 684 minerals (mol L-1 s-1) with respect to distance normal to fracture at different times at the 685 outlet of the core sample, for the 1 mL h-1 experiment (S-rich): calcite (Cal), gypsum 686 (Gp), clinochlore (Cln), albite (Ab), quartz (Qtz), pyrite (Py), anhydrite (Anh) and illite 687 (Ilt). 688 Figure 7 Variations of the simulated precipitation rates of the secondary minerals (mol 689 L-1 s-1) with respect to distance normal to fracture at different times in a 1 mL h-1 690 experiment (S-rich): dolomite (Dol), mesolite (Ms), smectite (Smc), kaolinite (Kln), 691 stilbite (Stl) and boehmite (Bhm). 692 Table 5 primary phases secondary phases minerals Cal Ilt Gp Cln Ab Kln Dol Ms Slb Bhm VTotal-diss [mm3] VTotalppt [mm3] S-rich solution exp 2; Q = 0.2 mL h - 1 Vmodel [mm 3 ] - 2.90 - 0.03 0.41 - 1.49 -0.01 0.01 0.78 0.19 1.35 0.13 -4.43 2.86 Vbalance [mm 3 ] - 2.80 nc 0.39 - 3.29 -(01.50) 2.134.34 nc nc nc nc -(6.097.59) 2.514.72 Vmodel/V balance 1.04 n c 1.05 0.45 exp 4; Q = 1 mL h - 1 Vmodel [mm 3 ] - 18.84 - 0.16 7.05 - 3.70 -0.18 0.002 2.01 0.38 2.37 0.31 -22.88 12.13 Vbalance [mm 3 ] - 20.34 nc 7.17 - 4.51 -(02.06) 2.695.73 nc nc nc nc -(24.8526.91) 9.8712.90 Vmodel/V balance 0.93 nc 0.98 0.82 exp 7; Q = 60 mL h - 1 Vmodel [mm 3 ] - 7.83 - 0.17 4.73 -1.2 1 -0.15 0.001 0.46 - 0.04 0.11 -9.36 6.13 Vbalance [mm 3 ] - 9.09 nc 13.98 - 1.17 -(00.53) 0.591.38 nc nc nc nc -(10.2510.79) 14.5715.36 Vmodel/V balance 0.86 nc 0.34 1.0 3 S-free solution exp 1; Q = 0.2 mL h - 1 Vmodel [mm 3 ] -3. 07 -0. 03 -0.1 5 -1. 64 -0.01 0.002 0. 83 0.12 1.21 0.15 -4.90 2.31 Vbalance [mm 3 ] - 3.07 nc - 0.14 - 2.29 -(01.05) 1.493.03 nc nc nc nc -(5.366.41) 1.633.17 Vmodel/V balance 1.00 nc 1.07 0.72 exp 3; Q = 1 mL h - 1 Vmodel [mm 3 ] -7.3 0 - 0.03 -0.8 3 -1.9 3 -0.02 0.002 0.8 7 0 .13 1.33 0.17 -10.11 2.51 Vbalance [mm 3 ] - 7.30 nc - 0.62 - 1.53 -(00.70) 0.961.99 nc nc nc nc -(8.839.53) 1.582.61 Vmodel/V balance 1.00 nc 1.34 1.26 exp 6; Q = 60 mL h - 1 Vmodel [mm 3 ] -44.3 1 - 0.1 5 -2. 52 -1.8 9 -0.15 0.001 0.5 5 0. 03 0. 90 0.16 -49.02 1.63 Vbalance [mm 3 ] - 45.82 nc nc - 1.98 -(00.90) (0.982.31) nc nc nc nc -(47.8048.70) 19.0020.33 Vmodel/V balance 0.97 nc nc 0.95 VTotal-diss is the total volume of dissolved mineral and VTotal-ppt is the total volume of precipitated mineral. V = Vfinal - Vinitial, where V > 0 and V < 0 indicate mineral precipitation and mineral dissolution, respectively. Values in parentheses indicate the V range considering the two hypothetical calculations. nc indicates that the volume could not be calculated from the mass balance (Ilt, Dol, Ms and Stl).         